chapter xiv, additional smoothing example, supplied by C. Reinsch calls bsplpp(bsplvb), ppvalu(interv), smooth(setupq,chol1d) c perturbed values of sin , with the opportunity of scaling both C x and y arbitrarily (which does bring ruin to the version of C smooth in the original pgs. integer i,j,npoint,nm1 parameter (npoint=91, nm1 = npoint-1) real skalx, skaly, skalesy(6) data skalesy /1.e-12, 1., 100., 300., 400., 500. / real a(npoint,4),dely,dy(npoint),freq, rands(91) * ,scrtch(npoint,7),sfp,sfpchk,x(npoint),xx,y(npoint) C the following random data were generated by matlab data rands / *-0.5621, -0.9059, 0.3577, 0.3586, 0.8694, -0.2330, 0.0388, *0.6619, -0.9309, -0.8931, 0.0594, 0.3423, -0.9846, -0.2332, *-0.8663, -0.1650, 0.3735, 0.1780, 0.8609, 0.6923, 0.0539, *-0.8161, 0.3078, -0.1680, 0.4024, 0.8206, 0.5244, -0.4751, *-0.9051, 0.4722, -0.3435, 0.2653, 0.5128, 0.9821, -0.2693, *-0.5059, 0.9651, 0.4453, 0.5067, 0.3030, -0.8546, 0.2633, *0.7694, -0.4546, -0.1272, 0.5330, -0.0445, -0.5245, -0.4502, *-0.2815, -0.6670, -0.0270, 0.7953, 0.8184, -0.8789, 0.8093, *0.0090, 0.0326, -0.3619, 0.9733, -0.0120, -0.4677, -0.8185, *0.8955, -0.8525, 0.0014, -0.2317, -0.4458, 0.8276, 0.0595, *-0.0711, 0.8820, -0.8998, 0.5230, 0.5404, 0.6556, -0.7493, *-0.9683, 0.3769, 0.7365, 0.2591, 0.4724, 0.4508, 0.9989, *0.7771, -0.5336, -0.3874, -0.2980, 0.0265, 0.1822, 0.6920/ freq = 3.14159265/9. skalx = 1. do 100 j=1,6 skaly = skalesy(j) dely = skaly*.05/sqrt(3.) do 10 i=1,npoint xx = float(i-1)*.1 x(i) = skalx*xx y(i) = skaly*(sin(xx*freq) + .05*rands(i)) 10 dy(i) = dely sfp = smooth ( x, y, dy, npoint, 80., scrtch, a ) c also compute sfp from scratch here sfpchk = 0. do 90 i=1,npoint sfpchk = sfpchk + ((y(i) - a(i,1))/dy(i))**2 90 continue print *, 'scales used: (',skalx,', ', skaly, '), sfp = ', sfp, sfpchk 100 continue c also try interpolation, to be sure: sfp = smooth ( x, y, dy, npoint, 0., scrtch, a ) sfpchk = 0. c the following check is trivial here since a(:,1) <-- y in this case. do 110 i=1,npoint sfpchk = sfpchk + ((y(i) - a(i,1))/dy(i))**2 110 continue print *, 'same scales as last, but S = 0, gives sfp = ', sfp, sfpchk stop end