REM Cubic Spline Interpolation

dim x(100)
dim y(100)
dim d(100)
dim f(100)
dim g(100)
dim h(100)

cls

print
print "Cubic Spline Interpolation"
print
print "Data must be input from the smallest x value to the largest..."
print
input "Number of data points to enter?  ";n
for r=1 to n
    [again]
    print
    print "x for point ";r;" = ";
    input " ";x1
    print "y for point ";r;" = ";
    input " ";y1
    x(r)=x1
    y(r)=y1
    if r<>1 and x(r)<x(r-1) then
        print "x values must be input from smallest to largest - try again"
        goto [again]
        end
    end if
next r

rem compute derivative
gosub [deriv]

rem compute cubic spline coefficients
gosub [csc]

rem this is the output routine
while q$<>"q"
    gosub [submenu]
    if q$="s" then
        gosub [single]
    else
        gosub [increment]
    end if
wend

[submenu]
cls
j=0
print
print "Do you want a single point (s),"
print "a list based on an increment (i),"
print "or quit (q)? "
print
input "s,i,q  ";q$
if q$="s" then
    cls
    print
    print "Enter x values one by one."
    print "Returns to submenu if value of x is outside"
    print "the range ";x(1);" to ";x(n)
    gosub [single]
end if
if q$="i" then
    gosub [increment]
end if
if q$="q" then end
gosub [submenu]
return

[single]
print
input " x=   ";j
if j<x(1) or j>x(n) then
    timer 3000,[jumpa]
    print "out of range"
    wait
    [jumpa]
    gosub [submenu]
end if
for k=2 to n
    if j<x(k) then
        k=k-1
        gosub [poly]
        print " y(x)= ";v
        gosub [single]
    end if
next k
return

[increment]
cls
print
input "Enter increment: ";u
print
j=x(1)
k=1
[more]
gosub [poly]
gosub [poly]
print j;
print "  ";
print v
j=j+u
if j<x(k+1) then
    goto [more]
end if
k=k+1
if k<>n then goto [more]
[more1]
print "press 'c' to continue..."
input " ";key$
if key$<>"c" then goto [more1]
return

[deriv]
rem compute the derivatives
j=2
gosub [coeff]
d(1)=2*a*x(1)+b
d(2)=2*a*x(2)+b
for j=3 to n-1
    gosub [coeff]
    d(j)=2*a*x(j)+b
next j
d(n)=2*a*x(n)+b
return

[coeff]
rem coefficents of a parabola
a=(y(j-1)-y(j))/(x(j-1)-x(j))
a=a-(y(j)-y(j+1))/(x(j)-x(j+1))
a=a/(x(j-1)-x(j+1))
b=(y(j-1)-y(j))/(x(j-1)-x(j))
b=b-a*(x(j-1)+x(j))
return

[csc]
rem cubic spline coefficients
for j=1 to n-1
    f(j)=d(j)*(x(j+1)-x(j))
    g(j)=3*y(j+1)-d(j+1)*(x(j+1)-x(j))
    g(j)=g(j)-3*y(j)-2*f(j)
    h(j)=y(j+1)-y(j)-f(j)-g(j)
next j
gosub [poly]
return

[poly]
rem evaluate polynomial
w=(j-x(k))/(x(k+1)-x(k))
v=y(k)+f(k)*w+g(k)*w*w+h(k)*w*w*w
return