1 subroutine bispev(tx,nx,ty,ny,c,kx,ky,x,mx,y,my,z,wrk,lwrk, 2 * iwrk,kwrk,ier) 3c subroutine bispev evaluates on a grid (x(i),y(j)),i=1,...,mx; j=1,... 4c ,my a bivariate spline s(x,y) of degrees kx and ky, given in the 5c b-spline representation. 6c 7c calling sequence: 8c call bispev(tx,nx,ty,ny,c,kx,ky,x,mx,y,my,z,wrk,lwrk, 9c * iwrk,kwrk,ier) 10c 11c input parameters: 12c tx : real array, length nx, which contains the position of the 13c knots in the x-direction. 14c nx : integer, giving the total number of knots in the x-direction 15c ty : real array, length ny, which contains the position of the 16c knots in the y-direction. 17c ny : integer, giving the total number of knots in the y-direction 18c c : real array, length (nx-kx-1)*(ny-ky-1), which contains the 19c b-spline coefficients. 20c kx,ky : integer values, giving the degrees of the spline. 21c x : real array of dimension (mx). 22c before entry x(i) must be set to the x co-ordinate of the 23c i-th grid point along the x-axis. 24c tx(kx+1)<=x(i-1)<=x(i)<=tx(nx-kx), i=2,...,mx. 25c mx : on entry mx must specify the number of grid points along 26c the x-axis. mx >=1. 27c y : real array of dimension (my). 28c before entry y(j) must be set to the y co-ordinate of the 29c j-th grid point along the y-axis. 30c ty(ky+1)<=y(j-1)<=y(j)<=ty(ny-ky), j=2,...,my. 31c my : on entry my must specify the number of grid points along 32c the y-axis. my >=1. 33c wrk : real array of dimension lwrk. used as workspace. 34c lwrk : integer, specifying the dimension of wrk. 35c lwrk >= mx*(kx+1)+my*(ky+1) 36c iwrk : integer array of dimension kwrk. used as workspace. 37c kwrk : integer, specifying the dimension of iwrk. kwrk >= mx+my. 38c 39c output parameters: 40c z : real array of dimension (mx*my). 41c on successful exit z(my*(i-1)+j) contains the value of s(x,y) 42c at the point (x(i),y(j)),i=1,...,mx;j=1,...,my. 43c ier : integer error flag 44c ier=0 : normal return 45c ier=10: invalid input data (see restrictions) 46c 47c restrictions: 48c mx >=1, my >=1, lwrk>=mx*(kx+1)+my*(ky+1), kwrk>=mx+my 49c tx(kx+1) <= x(i-1) <= x(i) <= tx(nx-kx), i=2,...,mx 50c ty(ky+1) <= y(j-1) <= y(j) <= ty(ny-ky), j=2,...,my 51c 52c other subroutines required: 53c fpbisp,fpbspl 54c 55c references : 56c de boor c : on calculating with b-splines, j. approximation theory 57c 6 (1972) 50-62. 58c cox m.g. : the numerical evaluation of b-splines, j. inst. maths 59c applics 10 (1972) 134-149. 60c dierckx p. : curve and surface fitting with splines, monographs on 61c numerical analysis, oxford university press, 1993. 62c 63c author : 64c p.dierckx 65c dept. computer science, k.u.leuven 66c celestijnenlaan 200a, b-3001 heverlee, belgium. 67c e-mail : Paul.Dierckx@cs.kuleuven.ac.be 68c 69c latest update : march 1987 70c 71c ..scalar arguments.. 72 integer nx,ny,kx,ky,mx,my,lwrk,kwrk,ier 73c ..array arguments.. 74 integer iwrk(kwrk) 75 real*8 tx(nx),ty(ny),c((nx-kx-1)*(ny-ky-1)),x(mx),y(my),z(mx*my), 76 * wrk(lwrk) 77c ..local scalars.. 78 integer i,iw,lwest 79c .. 80c before starting computations a data check is made. if the input data 81c are invalid control is immediately repassed to the calling program. 82 ier = 10 83 lwest = (kx+1)*mx+(ky+1)*my 84 if(lwrk.lt.lwest) go to 100 85 if(kwrk.lt.(mx+my)) go to 100 86 if (mx.lt.1) go to 100 87 if (mx.eq.1) go to 30 88 go to 10 89 10 do 20 i=2,mx 90 if(x(i).lt.x(i-1)) go to 100 91 20 continue 92 30 if (my.lt.1) go to 100 93 if (my.eq.1) go to 60 94 go to 40 95 40 do 50 i=2,my 96 if(y(i).lt.y(i-1)) go to 100 97 50 continue 98 60 ier = 0 99 iw = mx*(kx+1)+1 100 call fpbisp(tx,nx,ty,ny,c,kx,ky,x,mx,y,my,z,wrk(1),wrk(iw), 101 * iwrk(1),iwrk(mx+1)) 102 100 return 103 end 104