Contributed by: Doug Rickman, Global Hydrology and Climate Center, MSFC, NASA
This is aimed at intermediate REXX users.
Splines are very useful mathematical tools. These two routines take a set of irregularly spaced x,y points and return a set of y values at evenly spaced intervals. The code is based on the algorithms SPLINE and SPLINT in Numerical Recipes, FORTRAN version. To use store the x,y data in the paired arrays x. and y. (OK, if one is picky x. and y. are not arrays, but compound variables.!) The range of x values for which computed y values are to be obtained is from InitialLine2P to LastLine2P. The computed values are in the array Elevation..
n=X.0
/* Must have slopes at the beginning and ending of the line. */
YP1=(y.2-y.1)/(x.2-x.1)
nm1=n-1
YPN=(y.n-y.nm1)/(x.n-x.nm1)
call spline
do xe = InitialLine2P to LastLine2P
klo=1
khi=n
do while khi-klo>1
k=trunc((khi+klo)/2)
if X.k>Xe then khi=k
else klo=k
end /* while khi-klo>1 then do ... */
h=X.khi-X.klo
a=(X.khi-Xe)/h
b=(Xe-X.klo)/h
ye=a*Y.klo + b*Y.khi + ( (a*a*a -a) * Y2.klo + ( b*b*b -b)*Y2.khi )*h*h /6
Elevation.xe=strip(format(ye,8,1))
end xe
return 1
/* Given X.,Y.,n,YP1,YP2, return the array Y2. */
/* Based on the algorithm SPLINE in Numerical Recipes, FORTRAN version. */
/* X.,Y.=input table of observed x,y values. n=number of observed points. */
/* YP1 and YP2 are estimated slopes at first and last points. */
SPLINE:
if YP1>99e+30 then do
y2.1=0
u.1=0
end
else do
y2.1=-0.5
u.1=(3/(x.2-x.1))*((y.2-y.1)/(x.2-x.1)-YP1)
end
do i=2 to n-1
im1=i-1
ip1=i+1
sig=(x.i-x.im1)/(x.ip1-x.im1)
p=sig*Y2.im1+2
Y2.i=(sig-1)/p
u.i=(6*((y.iP1-y.i)/(x.ip1-x.i)-(y.i-y.im1) / (x.i-x.im1)) ,
/ (x.ip1-x.im1) - sig*u.im1)/p
end i
nm1=n-1
if (YPN>99e+30) then do
qn=0
un=0
end
else do
qn=0.5
un=(3/(x.n-x.nm1))*(YPN-(y.n-y.nm1)/(x.n-x.nm1))
end
Y2.n=(un-qn*u.nm1)/(qn*Y2.nm1-1)
do k=n-1 to 1 by -1
kp1=k+1
Y2.k=Y2.k*Y2.kp1+u.k
end k
return 1
end /* do */