chapter xvi, example 3. two parametrizations of some data c from * a practical guide to splines * by c. de boor calls splint(bsplvb,banfac/slv),bsplpp(bsplvb*),banslv*,ppvalu(interv) c parameter k=4,kpkm1=7, n=8,npk=12, npiece=6,npoint=21 c integer i,icount,iflag,kp1,l c real bcoef(n),break(npiece),ds,q(n,kpkm1),s(n),scrtch(k,k) c * ,ss,t(npk),x(n),xcoef(k,npiece),xx(npoint),y(n) c * ,ycoef(k,npiece),yy(npoint) integer i,icount,iflag,k,kpkm1,kp1,l,n,npoint data k,kpkm1,n,npoint /4,7,8,21/ real bcoef(8),break(6),ds,q(8,7),s(8),scrtch(4,4) * ,ss,t(12),x(8),xcoef(4,6),xx(21),y(8) * ,ycoef(4,6),yy(21) data x /0.,.1,.2,.3,.301,.4,.5,.6/ c *** compute y-component and set 'natural' parametrization do 1 i=1,n y(i) = (x(i)-.3)**2 1 s(i) = x(i) print 601 601 format(26h 'natural' parametrization/6x,1hx,11x,1hy) icount = 1 c *** convert data abscissae to knots. note that second and second c last data abscissae are not knots. 5 do 6 i=1,k t(i) = s(1) 6 t(n+i) = s(n) kp1 = k+1 do 7 i=kp1,n 7 t(i) = s(i+2-k) c *** interpolate to x-component call splint(s,x,t,n,k,q,bcoef,iflag) call bsplpp(t,bcoef,n,k,scrtch,break,xcoef,l) c *** interpolate to y-component. since data abscissae and knots are c the same for both components, we only need to use backsubstitution do 10 i=1,n 10 bcoef(i) = y(i) call banslv(q,kpkm1,n,k-1,k-1,bcoef) call bsplpp(t,bcoef,n,k,scrtch,break,ycoef,l) c *** evaluate curve at some points near the potential trouble spot, c the fourth and fifth data points. ss = s(3) ds = (s(6)-s(3))/float(npoint-1) do 20 i=1,npoint xx(i) = ppvalu(break,xcoef,l,k,ss,0) yy(i) = ppvalu(break,ycoef,l,k,ss,0) 20 ss = ss + ds print 620,(xx(i),yy(i),i=1,npoint) 620 format(2f12.7) if (icount .ge. 2) stop c *** now repeat the whole process with uniform parametrization icount = icount + 1 do 30 i=1,n 30 s(i) = float(i) print 630 630 format(/26h 'uniform' parametrization/6x,1hx,11x,1hy) go to 5 end