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