subroutine parder(tx,nx,ty,ny,c,kx,ky,nux,nuy,x,mx,y,my,z, * wrk,lwrk,iwrk,kwrk,ier) c subroutine parder evaluates on a grid (x(i),y(j)),i=1,...,mx; j=1,... c ,my the partial derivative ( order nux,nuy) of a bivariate spline c s(x,y) of degrees kx and ky, given in the b-spline representation. c c calling sequence: c call parder(tx,nx,ty,ny,c,kx,ky,nux,nuy,x,mx,y,my,z,wrk,lwrk, c * iwrk,kwrk,ier) c c input parameters: c tx : real array, length nx, which contains the position of the c knots in the x-direction. c nx : integer, giving the total number of knots in the x-direction c ty : real array, length ny, which contains the position of the c knots in the y-direction. c ny : integer, giving the total number of knots in the y-direction c c : real array, length (nx-kx-1)*(ny-ky-1), which contains the c b-spline coefficients. c kx,ky : integer values, giving the degrees of the spline. c nux : integer values, specifying the order of the partial c nuy derivative. 0<=nux=1. c y : real array of dimension (my). c before entry y(j) must be set to the y co-ordinate of the c j-th grid point along the y-axis. c ty(ky+1)<=y(j-1)<=y(j)<=ty(ny-ky), j=2,...,my. c my : on entry my must specify the number of grid points along c the y-axis. my >=1. c wrk : real array of dimension lwrk. used as workspace. c lwrk : integer, specifying the dimension of wrk. c lwrk >= mx*(kx+1-nux)+my*(ky+1-nuy)+(nx-kx-1)*(ny-ky-1) c iwrk : integer array of dimension kwrk. used as workspace. c kwrk : integer, specifying the dimension of iwrk. kwrk >= mx+my. c c output parameters: c z : real array of dimension (mx*my). c on succesful exit z(my*(i-1)+j) contains the value of the c specified partial derivative of s(x,y) at the point c (x(i),y(j)),i=1,...,mx;j=1,...,my. c ier : integer error flag c ier=0 : normal return c ier=10: invalid input data (see restrictions) c c restrictions: c mx >=1, my >=1, 0 <= nux < kx, 0 <= nuy < ky, kwrk>=mx+my c lwrk>=mx*(kx+1-nux)+my*(ky+1-nuy)+(nx-kx-1)*(ny-ky-1), c tx(kx+1) <= x(i-1) <= x(i) <= tx(nx-kx), i=2,...,mx c ty(ky+1) <= y(j-1) <= y(j) <= ty(ny-ky), j=2,...,my c c other subroutines required: c fpbisp,fpbspl c c references : c de boor c : on calculating with b-splines, j. approximation theory c 6 (1972) 50-62. c dierckx p. : curve and surface fitting with splines, monographs on c numerical analysis, oxford university press, 1993. c c author : c p.dierckx c dept. computer science, k.u.leuven c celestijnenlaan 200a, b-3001 heverlee, belgium. c e-mail : Paul.Dierckx@cs.kuleuven.ac.be c c latest update : march 1989 c c ..scalar arguments.. integer nx,ny,kx,ky,nux,nuy,mx,my,lwrk,kwrk,ier c ..array arguments.. integer iwrk(kwrk) real*8 tx(nx),ty(ny),c((nx-kx-1)*(ny-ky-1)),x(mx),y(my),z(mx*my), * wrk(lwrk) c ..local scalars.. integer i,iwx,iwy,j,kkx,kky,kx1,ky1,lx,ly,lwest,l1,l2,m,m0,m1, * nc,nkx1,nky1,nxx,nyy real*8 ak,fac c .. c before starting computations a data check is made. if the input data c are invalid control is immediately repassed to the calling program. ier = 10 kx1 = kx+1 ky1 = ky+1 nkx1 = nx-kx1 nky1 = ny-ky1 nc = nkx1*nky1 if(nux.lt.0 .or. nux.ge.kx) go to 400 if(nuy.lt.0 .or. nuy.ge.ky) go to 400 lwest = nc +(kx1-nux)*mx+(ky1-nuy)*my if(lwrk.lt.lwest) go to 400 if(kwrk.lt.(mx+my)) go to 400 if (mx.lt.1) go to 400 if (mx.eq.1) go to 30 go to 10 10 do 20 i=2,mx if(x(i).lt.x(i-1)) go to 400 20 continue 30 if (my.lt.1) go to 400 if (my.eq.1) go to 60 go to 40 40 do 50 i=2,my if(y(i).lt.y(i-1)) go to 400 50 continue 60 ier = 0 nxx = nkx1 nyy = nky1 kkx = kx kky = ky c the partial derivative of order (nux,nuy) of a bivariate spline of c degrees kx,ky is a bivariate spline of degrees kx-nux,ky-nuy. c we calculate the b-spline coefficients of this spline do 70 i=1,nc wrk(i) = c(i) 70 continue if(nux.eq.0) go to 200 lx = 1 do 100 j=1,nux ak = kkx nxx = nxx-1 l1 = lx m0 = 1 do 90 i=1,nxx l1 = l1+1 l2 = l1+kkx fac = tx(l2)-tx(l1) if(fac.le.0.) go to 90 do 80 m=1,nyy m1 = m0+nyy wrk(m0) = (wrk(m1)-wrk(m0))*ak/fac m0 = m0+1 80 continue 90 continue lx = lx+1 kkx = kkx-1 100 continue 200 if(nuy.eq.0) go to 300 ly = 1 do 230 j=1,nuy ak = kky nyy = nyy-1 l1 = ly do 220 i=1,nyy l1 = l1+1 l2 = l1+kky fac = ty(l2)-ty(l1) if(fac.le.0.) go to 220 m0 = i do 210 m=1,nxx m1 = m0+1 wrk(m0) = (wrk(m1)-wrk(m0))*ak/fac m0 = m0+nky1 210 continue 220 continue ly = ly+1 kky = kky-1 230 continue m0 = nyy m1 = nky1 do 250 m=2,nxx do 240 i=1,nyy m0 = m0+1 m1 = m1+1 wrk(m0) = wrk(m1) 240 continue m1 = m1+nuy 250 continue c we partition the working space and evaluate the partial derivative 300 iwx = 1+nxx*nyy iwy = iwx+mx*(kx1-nux) call fpbisp(tx(nux+1),nx-2*nux,ty(nuy+1),ny-2*nuy,wrk,kkx,kky, * x,mx,y,my,z,wrk(iwx),wrk(iwy),iwrk(1),iwrk(mx+1)) 400 return end