1
subroutine spgrid(iopt,ider,mu,u,mv,v,r,r0,r1,s,nuest,nvest,
2
* nu,tu,nv,tv,c,fp,wrk,lwrk,iwrk,kwrk,ier)
3
c given the function values r(i,j) on the latitude-longitude grid
4
c (u(i),v(j)), i=1,...,mu ; j=1,...,mv , spgrid determines a smooth
5
c bicubic spline approximation on the rectangular domain 0<=u<=pi,
6
c vb<=v<=ve (vb = v(1), ve=vb+2*pi).
7
c this approximation s(u,v) will satisfy the properties
9
c (1) s(0,v) = s(0,0) = dr(1)
11
c d s(0,v) d s(0,0) d s(0,pi/2)
12
c (2) -------- = cos(v)* -------- + sin(v)* -----------
15
c = cos(v)*dr(2)+sin(v)*dr(3)
17
c (3) s(pi,v) = s(pi,0) = dr(4)
19
c d s(pi,v) d s(pi,0) d s(pi,pi/2)
20
c (4) -------- = cos(v)* --------- + sin(v)* ------------
23
c = cos(v)*dr(5)+sin(v)*dr(6)
25
c and will be periodic in the variable v, i.e.
29
c (5) --------- = --------- 0 <=u<= pi , j=0,1,2
33
c the number of knots of s(u,v) and their position tu(i),i=1,2,...,nu;
34
c tv(j),j=1,2,...,nv, is chosen automatically by the routine. the
35
c smoothness of s(u,v) is achieved by minimalizing the discontinuity
36
c jumps of the derivatives of the spline at the knots. the amount of
37
c smoothness of s(u,v) is determined by the condition that
38
c fp=sumi=1,mu(sumj=1,mv((r(i,j)-s(u(i),v(j)))**2))+(r0-s(0,v))**2
39
c + (r1-s(pi,v))**2 <= s, with s a given non-negative constant.
40
c the fit s(u,v) is given in its b-spline representation and can be
41
c evaluated by means of routine bispev
44
c call spgrid(iopt,ider,mu,u,mv,v,r,r0,r1,s,nuest,nvest,nu,tu,
45
c * ,nv,tv,c,fp,wrk,lwrk,iwrk,kwrk,ier)
48
c iopt : integer array of dimension 3, specifying different options.
50
c iopt(1):on entry iopt(1) must specify whether a least-squares spline
51
c (iopt(1)=-1) or a smoothing spline (iopt(1)=0 or 1) must be
53
c if iopt(1)=0 the routine will start with an initial set of
54
c knots tu(i)=0,tu(i+4)=pi,i=1,...,4;tv(i)=v(1)+(i-4)*2*pi,
56
c if iopt(1)=1 the routine will continue with the set of knots
57
c found at the last call of the routine.
58
c attention: a call with iopt(1)=1 must always be immediately
59
c preceded by another call with iopt(1) = 1 or iopt(1) = 0.
60
c iopt(2):on entry iopt(2) must specify the requested order of conti-
61
c nuity at the pole u=0.
62
c if iopt(2)=0 only condition (1) must be fulfilled and
63
c if iopt(2)=1 conditions (1)+(2) must be fulfilled.
64
c iopt(3):on entry iopt(3) must specify the requested order of conti-
65
c nuity at the pole u=pi.
66
c if iopt(3)=0 only condition (3) must be fulfilled and
67
c if iopt(3)=1 conditions (3)+(4) must be fulfilled.
68
c ider : integer array of dimension 4, specifying different options.
70
c ider(1):on entry ider(1) must specify whether (ider(1)=0 or 1) or not
71
c (ider(1)=-1) there is a data value r0 at the pole u=0.
72
c if ider(1)=1, r0 will be considered to be the right function
73
c value, and it will be fitted exactly (s(0,v)=r0).
74
c if ider(1)=0, r0 will be considered to be a data value just
75
c like the other data values r(i,j).
76
c ider(2):on entry ider(2) must specify whether (ider(2)=1) or not
77
c (ider(2)=0) the approximation has vanishing derivatives
78
c dr(2) and dr(3) at the pole u=0 (in case iopt(2)=1)
79
c ider(3):on entry ider(3) must specify whether (ider(3)=0 or 1) or not
80
c (ider(3)=-1) there is a data value r1 at the pole u=pi.
81
c if ider(3)=1, r1 will be considered to be the right function
82
c value, and it will be fitted exactly (s(pi,v)=r1).
83
c if ider(3)=0, r1 will be considered to be a data value just
84
c like the other data values r(i,j).
85
c ider(4):on entry ider(4) must specify whether (ider(4)=1) or not
86
c (ider(4)=0) the approximation has vanishing derivatives
87
c dr(5) and dr(6) at the pole u=pi (in case iopt(3)=1)
88
c mu : integer. on entry mu must specify the number of grid points
89
c along the u-axis. unchanged on exit.
90
c mu >= 1, mu >=mumin=4-i0-i1-ider(2)-ider(4) with
91
c i0=min(1,ider(1)+1), i1=min(1,ider(3)+1)
92
c u : real array of dimension at least (mu). before entry, u(i)
93
c must be set to the u-co-ordinate of the i-th grid point
94
c along the u-axis, for i=1,2,...,mu. these values must be
95
c supplied in strictly ascending order. unchanged on exit.
97
c mv : integer. on entry mv must specify the number of grid points
98
c along the v-axis. mv > 3 . unchanged on exit.
99
c v : real array of dimension at least (mv). before entry, v(j)
100
c must be set to the v-co-ordinate of the j-th grid point
101
c along the v-axis, for j=1,2,...,mv. these values must be
102
c supplied in strictly ascending order. unchanged on exit.
103
c -pi <= v(1) < pi , v(mv) < v(1)+2*pi.
104
c r : real array of dimension at least (mu*mv).
105
c before entry, r(mv*(i-1)+j) must be set to the data value at
106
c the grid point (u(i),v(j)) for i=1,...,mu and j=1,...,mv.
108
c r0 : real value. on entry (if ider(1) >=0 ) r0 must specify the
109
c data value at the pole u=0. unchanged on exit.
110
c r1 : real value. on entry (if ider(1) >=0 ) r1 must specify the
111
c data value at the pole u=pi. unchanged on exit.
112
c s : real. on entry (if iopt(1)>=0) s must specify the smoothing
113
c factor. s >=0. unchanged on exit.
114
c for advice on the choice of s see further comments
115
c nuest : integer. unchanged on exit.
116
c nvest : integer. unchanged on exit.
117
c on entry, nuest and nvest must specify an upper bound for the
118
c number of knots required in the u- and v-directions respect.
119
c these numbers will also determine the storage space needed by
120
c the routine. nuest >= 8, nvest >= 8.
121
c in most practical situation nuest = mu/2, nvest=mv/2, will
122
c be sufficient. always large enough are nuest=mu+6+iopt(2)+
123
c iopt(3), nvest = mv+7, the number of knots needed for
124
c interpolation (s=0). see also further comments.
126
c unless ier=10 (in case iopt(1)>=0), nu will contain the total
127
c number of knots with respect to the u-variable, of the spline
128
c approximation returned. if the computation mode iopt(1)=1 is
129
c used, the value of nu should be left unchanged between sub-
130
c sequent calls. in case iopt(1)=-1, the value of nu should be
131
c specified on entry.
132
c tu : real array of dimension at least (nuest).
133
c on succesful exit, this array will contain the knots of the
134
c spline with respect to the u-variable, i.e. the position of
135
c the interior knots tu(5),...,tu(nu-4) as well as the position
136
c of the additional knots tu(1)=...=tu(4)=0 and tu(nu-3)=...=
137
c tu(nu)=pi needed for the b-spline representation.
138
c if the computation mode iopt(1)=1 is used,the values of tu(1)
139
c ...,tu(nu) should be left unchanged between subsequent calls.
140
c if the computation mode iopt(1)=-1 is used, the values tu(5),
141
c ...tu(nu-4) must be supplied by the user, before entry.
142
c see also the restrictions (ier=10).
144
c unless ier=10 (in case iopt(1)>=0), nv will contain the total
145
c number of knots with respect to the v-variable, of the spline
146
c approximation returned. if the computation mode iopt(1)=1 is
147
c used, the value of nv should be left unchanged between sub-
148
c sequent calls. in case iopt(1) = -1, the value of nv should
149
c be specified on entry.
150
c tv : real array of dimension at least (nvest).
151
c on succesful exit, this array will contain the knots of the
152
c spline with respect to the v-variable, i.e. the position of
153
c the interior knots tv(5),...,tv(nv-4) as well as the position
154
c of the additional knots tv(1),...,tv(4) and tv(nv-3),...,
155
c tv(nv) needed for the b-spline representation.
156
c if the computation mode iopt(1)=1 is used,the values of tv(1)
157
c ...,tv(nv) should be left unchanged between subsequent calls.
158
c if the computation mode iopt(1)=-1 is used, the values tv(5),
159
c ...tv(nv-4) must be supplied by the user, before entry.
160
c see also the restrictions (ier=10).
161
c c : real array of dimension at least (nuest-4)*(nvest-4).
162
c on succesful exit, c contains the coefficients of the spline
163
c approximation s(u,v)
164
c fp : real. unless ier=10, fp contains the sum of squared
165
c residuals of the spline approximation returned.
166
c wrk : real array of dimension (lwrk). used as workspace.
167
c if the computation mode iopt(1)=1 is used the values of
168
c wrk(1),..,wrk(12) should be left unchanged between subsequent
170
c lwrk : integer. on entry lwrk must specify the actual dimension of
171
c the array wrk as declared in the calling (sub)program.
172
c lwrk must not be too small.
173
c lwrk >= 12+nuest*(mv+nvest+3)+nvest*24+4*mu+8*mv+q
174
c where q is the larger of (mv+nvest) and nuest.
175
c iwrk : integer array of dimension (kwrk). used as workspace.
176
c if the computation mode iopt(1)=1 is used the values of
177
c iwrk(1),.,iwrk(5) should be left unchanged between subsequent
179
c kwrk : integer. on entry kwrk must specify the actual dimension of
180
c the array iwrk as declared in the calling (sub)program.
181
c kwrk >= 5+mu+mv+nuest+nvest.
182
c ier : integer. unless the routine detects an error, ier contains a
183
c non-positive value on exit, i.e.
184
c ier=0 : normal return. the spline returned has a residual sum of
185
c squares fp such that abs(fp-s)/s <= tol with tol a relat-
186
c ive tolerance set to 0.001 by the program.
187
c ier=-1 : normal return. the spline returned is an interpolating
189
c ier=-2 : normal return. the spline returned is the least-squares
190
c constrained polynomial. in this extreme case fp gives the
191
c upper bound for the smoothing factor s.
192
c ier=1 : error. the required storage space exceeds the available
193
c storage space, as specified by the parameters nuest and
195
c probably causes : nuest or nvest too small. if these param-
196
c eters are already large, it may also indicate that s is
198
c the approximation returned is the least-squares spline
199
c according to the current set of knots. the parameter fp
200
c gives the corresponding sum of squared residuals (fp>s).
201
c ier=2 : error. a theoretically impossible result was found during
202
c the iteration proces for finding a smoothing spline with
203
c fp = s. probably causes : s too small.
204
c there is an approximation returned but the corresponding
205
c sum of squared residuals does not satisfy the condition
207
c ier=3 : error. the maximal number of iterations maxit (set to 20
208
c by the program) allowed for finding a smoothing spline
209
c with fp=s has been reached. probably causes : s too small
210
c there is an approximation returned but the corresponding
211
c sum of squared residuals does not satisfy the condition
213
c ier=10 : error. on entry, the input data are controlled on validity
214
c the following restrictions must be satisfied.
215
c -1<=iopt(1)<=1, 0<=iopt(2)<=1, 0<=iopt(3)<=1,
216
c -1<=ider(1)<=1, 0<=ider(2)<=1, ider(2)=0 if iopt(2)=0.
217
c -1<=ider(3)<=1, 0<=ider(4)<=1, ider(4)=0 if iopt(3)=0.
218
c mu >= mumin (see above), mv >= 4, nuest >=8, nvest >= 8,
219
c kwrk>=5+mu+mv+nuest+nvest,
220
c lwrk >= 12+nuest*(mv+nvest+3)+nvest*24+4*mu+8*mv+
221
c max(nuest,mv+nvest)
222
c 0< u(i-1)<u(i)< pi,i=2,..,mu,
223
c -pi<=v(1)< pi, v(1)<v(i-1)<v(i)<v(1)+2*pi, i=3,...,mv
224
c if iopt(1)=-1: 8<=nu<=min(nuest,mu+6+iopt(2)+iopt(3))
225
c 0<tu(5)<tu(6)<...<tu(nu-4)< pi
226
c 8<=nv<=min(nvest,mv+7)
227
c v(1)<tv(5)<tv(6)<...<tv(nv-4)<v(1)+2*pi
228
c the schoenberg-whitney conditions, i.e. there must
229
c be subset of grid co-ordinates uu(p) and vv(q) such
230
c that tu(p) < uu(p) < tu(p+4) ,p=1,...,nu-4
231
c (iopt(2)=1 and iopt(3)=1 also count for a uu-value
232
c tv(q) < vv(q) < tv(q+4) ,q=1,...,nv-4
233
c (vv(q) is either a value v(j) or v(j)+2*pi)
234
c if iopt(1)>=0: s>=0
235
c if s=0: nuest>=mu+6+iopt(2)+iopt(3), nvest>=mv+7
236
c if one of these conditions is found to be violated,control
237
c is immediately repassed to the calling program. in that
238
c case there is no approximation returned.
241
c spgrid does not allow individual weighting of the data-values.
242
c so, if these were determined to widely different accuracies, then
243
c perhaps the general data set routine sphere should rather be used
244
c in spite of efficiency.
245
c by means of the parameter s, the user can control the tradeoff
246
c between closeness of fit and smoothness of fit of the approximation.
247
c if s is too large, the spline will be too smooth and signal will be
248
c lost ; if s is too small the spline will pick up too much noise. in
249
c the extreme cases the program will return an interpolating spline if
250
c s=0 and the constrained least-squares polynomial(degrees 3,0)if s is
251
c very large. between these extremes, a properly chosen s will result
252
c in a good compromise between closeness of fit and smoothness of fit.
253
c to decide whether an approximation, corresponding to a certain s is
254
c satisfactory the user is highly recommended to inspect the fits
256
c recommended values for s depend on the accuracy of the data values.
257
c if the user has an idea of the statistical errors on the data, he
258
c can also find a proper estimate for s. for, by assuming that, if he
259
c specifies the right s, spgrid will return a spline s(u,v) which
260
c exactly reproduces the function underlying the data he can evaluate
261
c the sum((r(i,j)-s(u(i),v(j)))**2) to find a good estimate for this s
262
c for example, if he knows that the statistical errors on his r(i,j)-
263
c values is not greater than 0.1, he may expect that a good s should
264
c have a value not larger than mu*mv*(0.1)**2.
265
c if nothing is known about the statistical error in r(i,j), s must
266
c be determined by trial and error, taking account of the comments
267
c above. the best is then to start with a very large value of s (to
268
c determine the least-squares polynomial and the corresponding upper
269
c bound fp0 for s) and then to progressively decrease the value of s
270
c ( say by a factor 10 in the beginning, i.e. s=fp0/10,fp0/100,...
271
c and more carefully as the approximation shows more detail) to
272
c obtain closer fits.
273
c to economize the search for a good s-value the program provides with
274
c different modes of computation. at the first call of the routine, or
275
c whenever he wants to restart with the initial set of knots the user
276
c must set iopt(1)=0.
277
c if iopt(1) = 1 the program will continue with the knots found at
278
c the last call of the routine. this will save a lot of computation
279
c time if spgrid is called repeatedly for different values of s.
280
c the number of knots of the spline returned and their location will
281
c depend on the value of s and on the complexity of the shape of the
282
c function underlying the data. if the computation mode iopt(1) = 1
283
c is used, the knots returned may also depend on the s-values at
284
c previous calls (if these were smaller). therefore, if after a number
285
c of trials with different s-values and iopt(1)=1,the user can finally
286
c accept a fit as satisfactory, it may be worthwhile for him to call
287
c spgrid once more with the chosen value for s but now with iopt(1)=0.
288
c indeed, spgrid may then return an approximation of the same quality
289
c of fit but with fewer knots and therefore better if data reduction
290
c is also an important objective for the user.
291
c the number of knots may also depend on the upper bounds nuest and
292
c nvest. indeed, if at a certain stage in spgrid the number of knots
293
c in one direction (say nu) has reached the value of its upper bound
294
c (nuest), then from that moment on all subsequent knots are added
295
c in the other (v) direction. this may indicate that the value of
296
c nuest is too small. on the other hand, it gives the user the option
297
c of limiting the number of knots the routine locates in any direction
298
c for example, by setting nuest=8 (the lowest allowable value for
299
c nuest), the user can indicate that he wants an approximation which
300
c is a simple cubic polynomial in the variable u.
302
c other subroutines required:
303
c fpspgr,fpchec,fpchep,fpknot,fpopsp,fprati,fpgrsp,fpsysy,fpback,
304
c fpbacp,fpbspl,fpcyt1,fpcyt2,fpdisc,fpgivs,fprota
307
c dierckx p. : fast algorithms for smoothing data over a disc or a
308
c sphere using tensor product splines, in "algorithms
309
c for approximation", ed. j.c.mason and m.g.cox,
310
c clarendon press oxford, 1987, pp. 51-65
311
c dierckx p. : fast algorithms for smoothing data over a disc or a
312
c sphere using tensor product splines, report tw73, dept.
313
c computer science,k.u.leuven, 1985.
314
c dierckx p. : curve and surface fitting with splines, monographs on
315
c numerical analysis, oxford university press, 1993.
319
c dept. computer science, k.u. leuven
320
c celestijnenlaan 200a, b-3001 heverlee, belgium.
321
c e-mail : Paul.Dierckx@cs.kuleuven.ac.be
323
c creation date : july 1985
324
c latest update : march 1989
327
c ..scalar arguments..
329
integer mu,mv,nuest,nvest,nu,nv,lwrk,kwrk,ier
330
c ..array arguments..
331
integer iopt(3),ider(4),iwrk(kwrk)
332
real*8 u(mu),v(mv),r(mu*mv),c((nuest-4)*(nvest-4)),tu(nuest),
333
* tv(nvest),wrk(lwrk)
335
real*8 per,pi,tol,uu,ve,rmax,rmin,one,half,rn,rb,re
336
integer i,i1,i2,j,jwrk,j1,j2,kndu,kndv,knru,knrv,kwest,l,
337
* ldr,lfpu,lfpv,lwest,lww,m,maxit,mumin,muu,nc
338
c ..function references..
341
c ..subroutine references..
342
c fpchec,fpchep,fpspgr
347
pi = datan2(0d0,-one)
350
c we set up the parameters tol and maxit.
353
c before starting computations, a data check is made. if the input data
354
c are invalid, control is immediately repassed to the calling program.
356
if(iopt(1).lt.(-1) .or. iopt(1).gt.1) go to 200
357
if(iopt(2).lt.0 .or. iopt(2).gt.1) go to 200
358
if(iopt(3).lt.0 .or. iopt(3).gt.1) go to 200
359
if(ider(1).lt.(-1) .or. ider(1).gt.1) go to 200
360
if(ider(2).lt.0 .or. ider(2).gt.1) go to 200
361
if(ider(2).eq.1 .and. iopt(2).eq.0) go to 200
362
if(ider(3).lt.(-1) .or. ider(3).gt.1) go to 200
363
if(ider(4).lt.0 .or. ider(4).gt.1) go to 200
364
if(ider(4).eq.1 .and. iopt(3).eq.0) go to 200
366
if(ider(1).ge.0) mumin = mumin-1
367
if(iopt(2).eq.1 .and. ider(2).eq.1) mumin = mumin-1
368
if(ider(3).ge.0) mumin = mumin-1
369
if(iopt(3).eq.1 .and. ider(4).eq.1) mumin = mumin-1
370
if(mumin.eq.0) mumin = 1
371
if(mu.lt.mumin .or. mv.lt.4) go to 200
372
if(nuest.lt.8 .or. nvest.lt.8) go to 200
374
nc = (nuest-4)*(nvest-4)
375
lwest = 12+nuest*(mv+nvest+3)+24*nvest+4*mu+8*mv+
376
* max0(nuest,mv+nvest)
377
kwest = 5+mu+mv+nuest+nvest
378
if(lwrk.lt.lwest .or. kwrk.lt.kwest) go to 200
379
if(u(1).le.0. .or. u(mu).ge.pi) go to 200
382
if(u(i-1).ge.u(i)) go to 200
384
30 if(v(1).lt. (-pi) .or. v(1).ge.pi ) go to 200
385
if(v(mv).ge.v(1)+per) go to 200
387
if(v(i-1).ge.v(i)) go to 200
389
if(iopt(1).gt.0) go to 140
390
c if not given, we compute an estimate for r0.
392
if(ider(1).lt.0) go to 45
400
c if not given, we compute an estimate for r1.
401
55 if(ider(3).lt.0) go to 60
411
c we determine the range of r-values.
415
if(r(i).lt.rmin) rmin = r(i)
416
if(r(i).gt.rmax) rmax = r(i)
428
if(iopt(1).eq.0) go to 140
429
if(nu.lt.8 .or. nu.gt.nuest) go to 200
430
if(nv.lt.11 .or. nv.gt.nvest) go to 200
439
if(iopt(2).eq.0) go to 100
442
if(uu.gt.tu(5)) uu = tu(5)
448
if(iopt(3).eq.0) go to 120
451
if(uu.lt.tu(nu-4)) uu = tu(nu-4)
452
wrk(l) = uu+(pi-uu)*half
456
call fpchec(wrk(13),muu,tu,nu,3,ier)
457
if(ier.ne.0) go to 200
478
call fpchep(wrk(13),mv+1,tv,nv,3,ier)
480
140 if(s.lt.0.) go to 200
481
if(s.eq.0. .and. (nuest.lt.(mu+6+iopt(2)+iopt(3)) .or.
482
* nvest.lt.(mv+7)) ) go to 200
483
c we partition the working space and determine the spline approximation
488
jwrk = lwrk-12-nuest-nvest
493
call fpspgr(iopt,ider,u,mu,v,mv,r,m,rb,re,s,nuest,nvest,tol,maxit,
495
* nc,nu,tu,nv,tv,c,fp,wrk(1),wrk(2),wrk(3),wrk(4),wrk(lfpu),
496
* wrk(lfpv),wrk(ldr),wrk(11),iwrk(1),iwrk(2),iwrk(3),iwrk(4),
497
* iwrk(5),iwrk(knru),iwrk(knrv),iwrk(kndu),iwrk(kndv),wrk(lww),
added
removed
Lib/interpolate/fitpack/spgrid.f
