1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
|
subroutine pardeu(tx,nx,ty,ny,c,kx,ky,nux,nuy,x,y,z,m,
* wrk,lwrk,iwrk,kwrk,ier)
c subroutine pardeu evaluates on a set of points (x(i),y(i)),i=1,...,m
c 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<kx, 0<=nuy<ky.
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 >= 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 (m).
c on succesful exit z(i) contains the value of the
c specified partial derivative of s(x,y) at the point
c (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 lwrk>=m*(kx+1-nux)+m*(ky+1-nuy)+(nx-kx-1)*(ny-ky-1),
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,m,lwrk,kwrk,ier,nux,nuy
c ..array arguments..
integer iwrk(kwrk)
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 i,iwx,iwy,j,kkx,kky,kx1,ky1,lx,ly,lwest,l1,l2,mm,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)*m+(ky1-nuy)*m
if(lwrk.lt.lwest) go to 400
if(kwrk.lt.(m+m)) go to 400
if (m.lt.1) go to 400
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 mm=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 mm=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 mm=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+m*(kx1-nux)
do 390 i=1,m
call fpbisp(tx(nux+1),nx-2*nux,ty(nuy+1),ny-2*nuy,wrk,kkx,kky,
* x(i),1,y(i),1,z(i),wrk(iwx),wrk(iwy),iwrk(1),iwrk(2))
390 continue
400 return
end
|
