subroutine sphere(iopt,m,teta,phi,r,w,s,ntest,npest,eps, * nt,tt,np,tp,c,fp,wrk1,lwrk1,wrk2,lwrk2,iwrk,kwrk,ier) c subroutine sphere determines a smooth bicubic spherical spline c approximation s(teta,phi), 0 <= teta <= pi ; 0 <= phi <= 2*pi c to a given set of data points (teta(i),phi(i),r(i)),i=1,2,...,m. c such a spline has the following specific properties c c (1) s(0,phi) = constant 0 <=phi<= 2*pi. c c (2) s(pi,phi) = constant 0 <=phi<= 2*pi c c j j c d s(teta,0) d s(teta,2*pi) c (3) ----------- = ------------ 0 <=teta<=pi, j=0,1,2 c j j c d phi d phi c c d s(0,phi) d s(0,0) d s(0,pi/2) c (4) ---------- = -------- *cos(phi) + ----------- *sin(phi) c d teta d teta d teta c c d s(pi,phi) d s(pi,0) d s(pi,pi/2) c (5) ----------- = ---------*cos(phi) + ------------*sin(phi) c d teta d teta d teta c c if iopt =-1 sphere calculates a weighted least-squares spherical c spline according to a given set of knots in teta- and phi- direction. c if iopt >=0, the number of knots in each direction and their position c tt(j),j=1,2,...,nt ; tp(j),j=1,2,...,np are chosen automatically by c the routine. the smoothness of s(teta,phi) is then achieved by mini- c malizing the discontinuity jumps of the derivatives of the spline c at the knots. the amount of smoothness of s(teta,phi) is determined c by the condition that fp = sum((w(i)*(r(i)-s(teta(i),phi(i))))**2) c be <= s, with s a given non-negative constant. c the spherical spline is given in the standard b-spline representation c of bicubic splines and can be evaluated by means of subroutine bispev c c calling sequence: c call sphere(iopt,m,teta,phi,r,w,s,ntest,npest,eps, c * nt,tt,np,tp,c,fp,wrk1,lwrk1,wrk2,lwrk2,iwrk,kwrk,ier) c c parameters: c iopt : integer flag. on entry iopt must specify whether a weighted c least-squares spherical spline (iopt=-1) or a smoothing c spherical spline (iopt=0 or 1) must be determined. c if iopt=0 the routine will start with an initial set of knots c tt(i)=0,tt(i+4)=pi,i=1,...,4;tp(i)=0,tp(i+4)=2*pi,i=1,...,4. c if iopt=1 the routine will continue with the set of knots c found at the last call of the routine. c attention: a call with iopt=1 must always be immediately pre- c ceded by another call with iopt=1 or iopt=0. c unchanged on exit. c m : integer. on entry m must specify the number of data points. c m >= 2. unchanged on exit. c teta : real array of dimension at least (m). c phi : real array of dimension at least (m). c r : real array of dimension at least (m). c before entry,teta(i),phi(i),r(i) must be set to the spherical c co-ordinates of the i-th data point, for i=1,...,m.the order c of the data points is immaterial. unchanged on exit. c w : real array of dimension at least (m). before entry, w(i) must c be set to the i-th value in the set of weights. the w(i) must c be strictly positive. unchanged on exit. c s : real. on entry (in case iopt>=0) s must specify the smoothing c factor. s >=0. unchanged on exit. c for advice on the choice of s see further comments c ntest : integer. unchanged on exit. c npest : integer. unchanged on exit. c on entry, ntest and npest must specify an upper bound for the c number of knots required in the teta- and phi-directions. c these numbers will also determine the storage space needed by c the routine. ntest >= 8, npest >= 8. c in most practical situation ntest = npest = 8+sqrt(m/2) will c be sufficient. see also further comments. c eps : real. c on entry, eps must specify a threshold for determining the c effective rank of an over-determined linear system of equat- c ions. 0 < eps < 1. if the number of decimal digits in the c computer representation of a real number is q, then 10**(-q) c is a suitable value for eps in most practical applications. c unchanged on exit. c nt : integer. c unless ier=10 (in case iopt >=0), nt will contain the total c number of knots with respect to the teta-variable, of the c spline approximation returned. if the computation mode iopt=1 c is used, the value of nt should be left unchanged between c subsequent calls. c in case iopt=-1, the value of nt should be specified on entry c tt : real array of dimension at least ntest. c on succesful exit, this array will contain the knots of the c spline with respect to the teta-variable, i.e. the position c of the interior knots tt(5),...,tt(nt-4) as well as the c position of the additional knots tt(1)=...=tt(4)=0 and c tt(nt-3)=...=tt(nt)=pi needed for the b-spline representation c if the computation mode iopt=1 is used, the values of tt(1), c ...,tt(nt) should be left unchanged between subsequent calls. c if the computation mode iopt=-1 is used, the values tt(5), c ...tt(nt-4) must be supplied by the user, before entry. c see also the restrictions (ier=10). c np : integer. c unless ier=10 (in case iopt >=0), np will contain the total c number of knots with respect to the phi-variable, of the c spline approximation returned. if the computation mode iopt=1 c is used, the value of np should be left unchanged between c subsequent calls. c in case iopt=-1, the value of np (>=9) should be specified c on entry. c tp : real array of dimension at least npest. c on succesful exit, this array will contain the knots of the c spline with respect to the phi-variable, i.e. the position of c the interior knots tp(5),...,tp(np-4) as well as the position c of the additional knots tp(1),...,tp(4) and tp(np-3),..., c tp(np) needed for the b-spline representation. c if the computation mode iopt=1 is used, the values of tp(1), c ...,tp(np) should be left unchanged between subsequent calls. c if the computation mode iopt=-1 is used, the values tp(5), c ...tp(np-4) must be supplied by the user, before entry. c see also the restrictions (ier=10). c c : real array of dimension at least (ntest-4)*(npest-4). c on succesful exit, c contains the coefficients of the spline c approximation s(teta,phi). c fp : real. unless ier=10, fp contains the weighted sum of c squared residuals of the spline approximation returned. c wrk1 : real array of dimension (lwrk1). used as workspace. c if the computation mode iopt=1 is used the value of wrk1(1) c should be left unchanged between subsequent calls. c on exit wrk1(2),wrk1(3),...,wrk1(1+ncof) will contain the c values d(i)/max(d(i)),i=1,...,ncof=6+(np-7)*(nt-8) c with d(i) the i-th diagonal element of the reduced triangular c matrix for calculating the b-spline coefficients. it includes c those elements whose square is less than eps,which are treat- c ed as 0 in the case of presumed rank deficiency (ier<-2). c lwrk1 : integer. on entry lwrk1 must specify the actual dimension of c the array wrk1 as declared in the calling (sub)program. c lwrk1 must not be too small. let c u = ntest-7, v = npest-7, then c lwrk1 >= 185+52*v+10*u+14*u*v+8*(u-1)*v**2+8*m c wrk2 : real array of dimension (lwrk2). used as workspace, but c only in the case a rank deficient system is encountered. c lwrk2 : integer. on entry lwrk2 must specify the actual dimension of c the array wrk2 as declared in the calling (sub)program. c lwrk2 > 0 . a save upper bound for lwrk2 = 48+21*v+7*u*v+ c 4*(u-1)*v**2 where u,v are as above. if there are enough data c points, scattered uniformly over the approximation domain c and if the smoothing factor s is not too small, there is a c good chance that this extra workspace is not needed. a lot c of memory might therefore be saved by setting lwrk2=1. c (see also ier > 10) c iwrk : integer array of dimension (kwrk). used as workspace. c kwrk : integer. on entry kwrk must specify the actual dimension of c the array iwrk as declared in the calling (sub)program. c kwrk >= m+(ntest-7)*(npest-7). c ier : integer. unless the routine detects an error, ier contains a c non-positive value on exit, i.e. c ier=0 : normal return. the spline returned has a residual sum of c squares fp such that abs(fp-s)/s <= tol with tol a relat- c ive tolerance set to 0.001 by the program. c ier=-1 : normal return. the spline returned is a spherical c interpolating spline (fp=0). c ier=-2 : normal return. the spline returned is the weighted least- c squares constrained polynomial . in this extreme case c fp gives the upper bound for the smoothing factor s. c ier<-2 : warning. the coefficients of the spline returned have been c computed as the minimal norm least-squares solution of a c (numerically) rank deficient system. (-ier) gives the rank. c especially if the rank deficiency which can be computed as c 6+(nt-8)*(np-7)+ier, is large the results may be inaccurate c they could also seriously depend on the value of eps. c ier=1 : error. the required storage space exceeds the available c storage space, as specified by the parameters ntest and c npest. c probably causes : ntest or npest too small. if these param- c eters are already large, it may also indicate that s is c too small c the approximation returned is the weighted least-squares c spherical spline according to the current set of knots. c the parameter fp gives the corresponding weighted sum of c squared residuals (fp>s). c ier=2 : error. a theoretically impossible result was found during c the iteration proces for finding a smoothing spline with c fp = s. probably causes : s too small or badly chosen eps. c there is an approximation returned but the corresponding c weighted sum of squared residuals does not satisfy the c condition abs(fp-s)/s < tol. c ier=3 : error. the maximal number of iterations maxit (set to 20 c by the program) allowed for finding a smoothing spline c with fp=s has been reached. probably causes : s too small c there is an approximation returned but the corresponding c weighted sum of squared residuals does not satisfy the c condition abs(fp-s)/s < tol. c ier=4 : error. no more knots can be added because the dimension c of the spherical spline 6+(nt-8)*(np-7) already exceeds c the number of data points m. c probably causes : either s or m too small. c the approximation returned is the weighted least-squares c spherical spline according to the current set of knots. c the parameter fp gives the corresponding weighted sum of c squared residuals (fp>s). c ier=5 : error. no more knots can be added because the additional c knot would (quasi) coincide with an old one. c probably causes : s too small or too large a weight to an c inaccurate data point. c the approximation returned is the weighted least-squares c spherical spline according to the current set of knots. c the parameter fp gives the corresponding weighted sum of c squared residuals (fp>s). c ier=10 : error. on entry, the input data are controlled on validity c the following restrictions must be satisfied. c -1<=iopt<=1, m>=2, ntest>=8 ,npest >=8, 00, i=1,...,m c lwrk1 >= 185+52*v+10*u+14*u*v+8*(u-1)*v**2+8*m c kwrk >= m+(ntest-7)*(npest-7) c if iopt=-1: 8<=nt<=ntest , 9<=np<=npest c 0=0: s>=0 c if one of these conditions is found to be violated,control c is immediately repassed to the calling program. in that c case there is no approximation returned. c ier>10 : error. lwrk2 is too small, i.e. there is not enough work- c space for computing the minimal least-squares solution of c a rank deficient system of linear equations. ier gives the c requested value for lwrk2. there is no approximation re- c turned but, having saved the information contained in nt, c np,tt,tp,wrk1, and having adjusted the value of lwrk2 and c the dimension of the array wrk2 accordingly, the user can c continue at the point the program was left, by calling c sphere with iopt=1. c c further comments: c by means of the parameter s, the user can control the tradeoff c between closeness of fit and smoothness of fit of the approximation. c if s is too large, the spline will be too smooth and signal will be c lost ; if s is too small the spline will pick up too much noise. in c the extreme cases the program will return an interpolating spline if c s=0 and the constrained weighted least-squares polynomial if s is c very large. between these extremes, a properly chosen s will result c in a good compromise between closeness of fit and smoothness of fit. c to decide whether an approximation, corresponding to a certain s is c satisfactory the user is highly recommended to inspect the fits c graphically. c recommended values for s depend on the weights w(i). if these are c taken as 1/d(i) with d(i) an estimate of the standard deviation of c r(i), a good s-value should be found in the range (m-sqrt(2*m),m+ c sqrt(2*m)). if nothing is known about the statistical error in r(i) c each w(i) can be set equal to one and s determined by trial and c error, taking account of the comments above. the best is then to c start with a very large value of s ( to determine the least-squares c polynomial and the corresponding upper bound fp0 for s) and then to c progressively decrease the value of s ( say by a factor 10 in the c beginning, i.e. s=fp0/10, fp0/100,...and more carefully as the c approximation shows more detail) to obtain closer fits. c to choose s very small is strongly discouraged. this considerably c increases computation time and memory requirements. it may also c cause rank-deficiency (ier<-2) and endager numerical stability. c to economize the search for a good s-value the program provides with c different modes of computation. at the first call of the routine, or c whenever he wants to restart with the initial set of knots the user c must set iopt=0. c if iopt=1 the program will continue with the set of knots found at c the last call of the routine. this will save a lot of computation c time if sphere is called repeatedly for different values of s. c the number of knots of the spline returned and their location will c depend on the value of s and on the complexity of the shape of the c function underlying the data. if the computation mode iopt=1 c is used, the knots returned may also depend on the s-values at c previous calls (if these were smaller). therefore, if after a number c of trials with different s-values and iopt=1, the user can finally c accept a fit as satisfactory, it may be worthwhile for him to call c sphere once more with the selected value for s but now with iopt=0. c indeed, sphere may then return an approximation of the same quality c of fit but with fewer knots and therefore better if data reduction c is also an important objective for the user. c the number of knots may also depend on the upper bounds ntest and c npest. indeed, if at a certain stage in sphere the number of knots c in one direction (say nt) has reached the value of its upper bound c (ntest), then from that moment on all subsequent knots are added c in the other (phi) direction. this may indicate that the value of c ntest is too small. on the other hand, it gives the user the option c of limiting the number of knots the routine locates in any direction c for example, by setting ntest=8 (the lowest allowable value for c ntest), the user can indicate that he wants an approximation which c is a cubic polynomial in the variable teta. c c other subroutines required: c fpback,fpbspl,fpsphe,fpdisc,fpgivs,fprank,fprati,fprota,fporde, c fprpsp c c references: c dierckx p. : algorithms for smoothing data on the sphere with tensor c product splines, computing 32 (1984) 319-342. c dierckx p. : algorithms for smoothing data on the sphere with tensor c product splines, report tw62, dept. computer science, c k.u.leuven, 1983. 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 creation date : july 1983 c latest update : march 1989 c c .. c ..scalar arguments.. real s,eps,fp integer iopt,m,ntest,npest,nt,np,lwrk1,lwrk2,kwrk,ier c ..array arguments.. real teta(m),phi(m),r(m),w(m),tt(ntest),tp(npest), * c((ntest-4)*(npest-4)),wrk1(lwrk1),wrk2(lwrk2) integer iwrk(kwrk) c ..local scalars.. real tol,pi,pi2,one integer i,ib1,ib3,ki,kn,kwest,la,lbt,lcc,lcs,lro,j * lbp,lco,lf,lff,lfp,lh,lq,lst,lsp,lwest,maxit,ncest,ncc,ntt, * npp,nreg,nrint,ncof,nt4,np4 c ..function references.. real atan c ..subroutine references.. c fpsphe c .. c set constants one = 0.1e+01 c we set up the parameters tol and maxit. maxit = 20 tol = 0.1e-02 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 if(eps.le.0. .or. eps.ge.1.) go to 80 if(iopt.lt.(-1) .or. iopt.gt.1) go to 80 if(m.lt.2) go to 80 if(ntest.lt.8 .or. npest.lt.8) go to 80 nt4 = ntest-4 np4 = npest-4 ncest = nt4*np4 ntt = ntest-7 npp = npest-7 ncc = 6+npp*(ntt-1) nrint = ntt+npp nreg = ntt*npp ncof = 6+3*npp ib1 = 4*npp ib3 = ib1+3 if(ncof.gt.ib1) ib1 = ncof if(ncof.gt.ib3) ib3 = ncof lwest = 185+52*npp+10*ntt+14*ntt*npp+8*(m+(ntt-1)*npp**2) kwest = m+nreg if(lwrk1.lt.lwest .or. kwrk.lt.kwest) go to 80 if(iopt.gt.0) go to 60 pi = atan(one)*4 pi2 = pi+pi do 20 i=1,m if(w(i).le.0.) go to 80 if(teta(i).lt.0. .or. teta(i).gt.pi) go to 80 if(phi(i) .lt.0. .or. phi(i).gt.pi2) go to 80 20 continue if(iopt.eq.0) go to 60 ntt = nt-8 if(ntt.lt.0 .or. nt.gt.ntest) go to 80 if(ntt.eq.0) go to 40 tt(4) = 0. do 30 i=1,ntt j = i+4 if(tt(j).le.tt(j-1) .or. tt(j).ge.pi) go to 80 30 continue 40 npp = np-8 if(npp.lt.1 .or. np.gt.npest) go to 80 tp(4) = 0. do 50 i=1,npp j = i+4 if(tp(j).le.tp(j-1) .or. tp(j).ge.pi2) go to 80 50 continue go to 70 60 if(s.lt.0.) go to 80 70 ier = 0 c we partition the working space and determine the spline approximation kn = 1 ki = kn+m lq = 2 la = lq+ncc*ib3 lf = la+ncc*ib1 lff = lf+ncc lfp = lff+ncest lco = lfp+nrint lh = lco+nrint lbt = lh+ib3 lbp = lbt+5*ntest lro = lbp+5*npest lcc = lro+npest lcs = lcc+npest lst = lcs+npest lsp = lst+m*4 call fpsphe(iopt,m,teta,phi,r,w,s,ntest,npest,eps,tol,maxit, * ib1,ib3,ncest,ncc,nrint,nreg,nt,tt,np,tp,c,fp,wrk1(1),wrk1(lfp), * wrk1(lco),wrk1(lf),wrk1(lff),wrk1(lro),wrk1(lcc),wrk1(lcs), * wrk1(la),wrk1(lq),wrk1(lbt),wrk1(lbp),wrk1(lst),wrk1(lsp), * wrk1(lh),iwrk(ki),iwrk(kn),wrk2,lwrk2,ier) 80 return end