67 lines
2.3 KiB
Fortran
67 lines
2.3 KiB
Fortran
recursive subroutine bispeu(tx,nx,ty,ny,c,kx,ky,x,y,z,m,wrk,
|
|
* lwrk, ier)
|
|
implicit none
|
|
c subroutine bispeu evaluates on a set of points (x(i),y(i)),i=1,...,m
|
|
c a bivariate spline s(x,y) of degrees kx and ky, given in the
|
|
c b-spline representation.
|
|
c
|
|
c calling sequence:
|
|
c call bispeu(tx,nx,ty,ny,c,kx,ky,x,y,z,m,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 x : real array of dimension (mx).
|
|
c y : real array of dimension (my).
|
|
c m : on entry m must specify the number points. m >= 1.
|
|
c wrk : real array of dimension lwrk. used as workspace.
|
|
c lwrk : integer, specifying the dimension of wrk.
|
|
c lwrk >= kx+ky+2
|
|
c
|
|
c output parameters:
|
|
c z : real array of dimension m.
|
|
c on successful exit z(i) contains the value of s(x,y)
|
|
c at the point (x(i),y(i)), i=1,...,m.
|
|
c ier : integer error flag
|
|
c ier=0 : normal return
|
|
c ier=10: invalid input data (see restrictions)
|
|
c
|
|
c restrictions:
|
|
c m >=1, lwrk>=mx*(kx+1)+my*(ky+1), kwrk>=mx+my
|
|
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 ..scalar arguments..
|
|
integer nx,ny,kx,ky,m,lwrk,ier
|
|
c ..array arguments..
|
|
real*8 tx(nx),ty(ny),c((nx-kx-1)*(ny-ky-1)),x(m),y(m),z(m),
|
|
* wrk(lwrk)
|
|
c ..local scalars..
|
|
integer iwrk(2)
|
|
integer i, lwest
|
|
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
|
|
lwest = kx+ky+2
|
|
if (lwrk.lt.lwest) go to 100
|
|
if (m.lt.1) go to 100
|
|
ier = 0
|
|
do 10 i=1,m
|
|
call fpbisp(tx,nx,ty,ny,c,kx,ky,x(i),1,y(i),1,z(i),wrk(1),
|
|
* wrk(kx+2),iwrk(1),iwrk(2))
|
|
10 continue
|
|
100 return
|
|
end
|