' demos_sfit.bas April 1, 2017

' this program demonstrates how to interact
' with the cubic spline subroutines which
' "fit" and interpolate one-dimensional
' tabular data of the form y = f(x)

' MMBASIC eXtreme version

'''''''''''''''''''''''''

option default float

option base 1

dim x(168), y(168), xk(168), yk(168)

dim xfit(168), yfit(168)

dim c1(168), c2(168), c3(168)

' read tabular data

readdata

' define number of x and y data pairs

ndata% = 168

' sample every fifth data value and load
' "working" data arrays for spline operations

nsp% = 5

spdata(x(), y(), ndata%, nsp%, xk(), yk(), npts%)

' generate cubic spline coefficient arrays

spcoef(xk(), yk(), npts%, c1(), c2(), c3())

' evaluate the cubic spline "fit" at every third
' data point and display the interpolated data

nspj% = 3

i% = 0

for jj% = 1 to ndata%

if ((jj% mod nspj%) <> 0) then

' null

else

ii% = ii% + 1

xval = x(jj%)

spval(xval, xk(), yk(), npts%, c1(), c2(), c3(), sx)

xfit(ii%) = xval

yfit(ii%) = sx

print ii%, xfit(ii%), yfit(ii%)

end if

next jj%

end

''''''''''''''''''''''''''''''''''''''''''''''''''''
''''''''''''''''''''''''''''''''''''''''''''''''''''

sub spdata (x(), y(), ndata%, nsp%, xk(), yk(), npts%)

' sample and load data arrays for spline operations subroutine

' input

' x() = x data array
' y() = y data array
' ndata% = number of data points in data arrays
' nsp% = data points sampling interval

' output

' xk() = sampled x data array
' yk() = sampled y data array
' npts% = number of data points in sampled arrays

''''''''''''''''''''''''''''''''''''''''''''''''''

xk(1) = x(1)

yk(1) = y(1)

npts% = 1

for i% = 1 to ndata%

if ((i% mod nsp%) <> 0)) then

' null

else

npts% = npts% + 1

xk(npts%) = x(i%)

yk(npts%) = y(i%)

j% = i%

end if

next i%

if (j% = ndata%) then

' null

else

npts% = npts% + 1

xk(npts%) = x(ndata%)

yk(npts%) = y(ndata%)

end if

end sub

'''''''''''''''''''''''''''''''''''''''''
'''''''''''''''''''''''''''''''''''''''''

sub spcoef (x(), y(), n%, c1(), c2(), c3())

' calculate spline coefficients subroutine

' input

' x() = x data array
' y() = y data array
' n% = number of data points in arrays

' output

' c1(), c2(), c3() = array of spline coefficients

''''''''''''''''''''''''''''''''''''''''''''''''''

if (n% < 2) then

print "not enough data points; n must be >= 2"

return

end if

if (n% = 2) then

c1(1) = (y(2) - y(1)) / (x(2) - x(1))

c2(1) = 0.0

c3(1) = 0.0

c1(2) = c1(1)

c2(2) = 0.0

c3(2) = 0.0

return

end if

nm1% = n% - 1

c3(1) = x(2) - x(1)

c2(2) = (y(2) - y(1)) / c3(1)

for i% = 2 to nm1%

c3(i%) = x(i% + 1) - x(i%)

c1(i%) = 2 * (c3(i% - 1) + c3(i%))

c2(i% + 1) = (y(i% + 1) - y(i%)) / c3(i%)

c2(i%) = c2(i% + 1) - c2(i%)

next i%

c1(1) = -c3(1)

c1(n%) = -c3(n% - 1)

if (n% <> 3) then

c2(1) = c2(3) / (x(4) - x(2)) - c2(2) / (x(3) - x(1))

c2(n%) = c2(n% - 1) / (x(n%) - x(n% - 2)) - c2(n% - 2) / (x(n% - 1) - x(n% - 3))

c2(1) = c2(1) * c3(1) ^ 2 / (x(4) - x(1))

c2(n%) = -c2(n%) * c3(n% - 1) ^ 2 / (x(n%) - x(n% - 3))

else

c2(1) = 0.0

c2(n%) = 0.0

end if

for i% = 2 to n%

t = c3(i% - 1) / c1(i% - 1)

c1(i%) = c1(i%) - t * c3(i% - 1)

c2(i%) = c2(i%) - t * c2(i% - 1)

next i%

c2(n%) = c2(n%) / c1(n%)

for j% = 1 to nm1%

i% = n% - j%

c2(i%) = (c2(i%) - c3(i%) * c2(i% + 1)) / c1(i%)

next j%

c1(n%) = (y(n%) - y(nm1%)) / c3(nm1%) + c3(nm1%) * (c2(nm1%) + 2 * c2(n%))

for i% = 1 to nm1%

c1(i%) = (y(i% + 1) - y(i%)) / c3(i%) - c3(i%) * (c2(i% + 1) + 2 * c2(i%))

c3(i%) = (c2(i% + 1) - c2(i%)) / c3(i%)

c2(i%) = 3.0 * c2(i%)

next i%

c2(n%) = 3.0 * c2(n%)

c3(n%) = c3(n% - 1)

end sub

'''''''''''''''''''''''''''''''''''''''''''''''
'''''''''''''''''''''''''''''''''''''''''''''''

sub spval (u, x(), y(), n%, c1(), c2(), c3(), sx)

' evaluate spline fit subroutine

' input

' u = x value to be "fit"
' x() = x data array
' y() = y data array
' n% = number of data points in arrays
' c1(), c2(), c3() = spline fit coefficient arrays

' output

' sx = spline "fit" at u value

'''''''''''''''''''''''''''''''

i% = 1

if ((n% < 2) or (u > x(n%)) or (u < x(1))) then

print " "

print "not enough data points or x value outside range"

return

end if

if (i% >= n%) then

i% = 1

end if

if (u <= x(i% + 1)) then

' null

else

i% = 1

j% = n% + 1

do

k% = int((i% + j%) / 2)

if (u < x(k%)) then

j% = k%

end if

if (u >= x(k%)) then

i% = k%

end if

if (j% > (i% + 1)) then

' null

else

exit do

end if

loop

end if

dx = u - x(i%)

c1i = c1(i%)

c2i = c2(i%)

c3i = c3(i%)

prod = dx * (c1i + dx * (c2i + dx * c3i))

sx = y(i%) + prod

end sub

''''''''''
''''''''''

sub readdata

' read tabular data subroutine

for i% = 1 to 168

read y(i%), x(i%)

next i%

data 80.4, 1994.9170
data 79.6, 1995.0004
data 78.8, 1995.0836
data 78.2, 1995.1670
data 77.7, 1995.2504
data 77.3, 1995.3336
data 76.7, 1995.4170
data 76.0, 1995.5004
data 75.3, 1995.5836
data 74.5, 1995.6670
data 73.7, 1995.7504
data 73.2, 1995.8336
data 72.9, 1995.9170
data 72.5, 1996.0004
data 72.3, 1996.0836
data 72.2, 1996.1670
data 72.1, 1996.2504
data 72.1, 1996.3336
data 72.3, 1996.4170
data 72.6, 1996.5004
data 72.9, 1996.5836
data 73.5, 1996.6670
data 74.1, 1996.7504
data 75.1, 1996.8336
data 76.0, 1996.9170
data 77.0, 1997.0004
data 78.0, 1997.0836
data 79.3, 1997.1670
data 80.6, 1997.2504
data 82.0, 1997.3336
data 83.7, 1997.4170
data 85.6, 1997.5004
data 87.4, 1997.5836
data 89.3, 1997.6670
data 91.4, 1997.7504
data 93.7, 1997.8336
data 95.9, 1997.9170
data 98.2, 1998.0004
data 100.6, 1998.0836
data 103.1, 1998.1670
data 105.9, 1998.2504
data 108.9, 1998.3336
data 111.8, 1998.4170
data 114.8, 1998.5004
data 118.0, 1998.5836
data 120.9, 1998.6670
data 123.7, 1998.7504
data 126.5, 1998.8336
data 129.2, 1998.9170
data 132.0, 1999.0004
data 134.6, 1999.0836
data 137.2, 1999.1670
data 139.6, 1999.2504
data 142.0, 1999.3336
data 144.3, 1999.4170
data 146.4, 1999.5004
data 148.4, 1999.5836
data 150.2, 1999.6670
data 151.9, 1999.7504
data 153.4, 1999.8336
data 154.7, 1999.9170
data 155.8, 2000.0004
data 156.6, 2000.0836
data 157.3, 2000.1670
data 157.6, 2000.2504
data 157.8, 2000.3336
data 155.9, 2000.4170
data 154.2, 2000.5004
data 152.8, 2000.5836
data 151.5, 2000.6670
data 150.0, 2000.7504
data 148.4, 2000.8336
data 146.4, 2000.9170
data 144.7, 2001.0004
data 143.1, 2001.0836
data 141.4, 2001.1670
data 139.5, 2001.2504
data 137.9, 2001.3336
data 137.0, 2001.4170
data 136.0, 2001.5004
data 134.2, 2001.5836
data 132.5, 2001.6670
data 131.3, 2001.7504
data 130.1, 2001.8336
data 128.7, 2001.9170
data 127.1, 2002.0004
data 125.3, 2002.0836
data 123.5, 2002.1670
data 122.0, 2002.2504
data 120.6, 2002.3336
data 118.7, 2002.4170
data 116.5, 2002.5004
data 114.5, 2002.5836
data 112.9, 2002.6670
data 111.6, 2002.7504
data 110.0, 2002.8336
data 108.1, 2002.9170
data 106.6, 2003.0004
data 105.1, 2003.0836
data 103.5, 2003.1670
data 102.1, 2003.2504
data 100.9, 2003.3336
data 99.7, 2003.4170
data 98.6, 2003.5004
data 97.6, 2003.5836
data 96.7, 2003.6670
data 95.8, 2003.7504
data 94.6, 2003.8336
data 93.5, 2003.9170
data 92.3, 2004.0004
data 91.4, 2004.0836
data 90.6, 2004.1670
data 89.7, 2004.2504
data 88.6, 2004.3336
data 87.8, 2004.4170
data 87.0, 2004.5004
data 86.0, 2004.5836
data 84.8, 2004.6670
data 83.6, 2004.7504
data 82.5, 2004.8336
data 81.9, 2004.9170
data 81.3, 2005.0004
data 80.4, 2005.0836
data 79.6, 2005.1670
data 79.0, 2005.2504
data 78.4, 2005.3336
data 77.8, 2005.4170
data 77.3, 2005.5004
data 77.0, 2005.5836
data 76.6, 2005.6670
data 76.1, 2005.7504
data 75.6, 2005.8336
data 75.1, 2005.9170
data 74.6, 2006.0004
data 74.2, 2006.0836
data 74.0, 2006.1670
data 73.7, 2006.2504
data 73.4, 2006.3336
data 73.1, 2006.4170
data 72.6, 2006.5004
data 72.3, 2006.5836
data 71.9, 2006.6670
data 71.7, 2006.7504
data 71.6, 2006.8336
data 71.5, 2006.9170
data 71.5, 2007.0004
data 71.7, 2007.0836
data 71.8, 2007.1670
data 72.1, 2007.2504
data 72.5, 2007.3336
data 73.0, 2007.4170
data 73.4, 2007.5004
data 74.1, 2007.5836
data 74.9, 2007.6670
data 76.0, 2007.7504
data 77.2, 2007.8336
data 78.3, 2007.9170
data 79.6, 2008.0004
data 80.9, 2008.0836
data 82.4, 2008.1670
data 84.1, 2008.2504
data 86.2, 2008.3336
data 88.2, 2008.4170
data 90.3, 2008.5004
data 92.4, 2008.5836
data 94.5, 2008.6670
data 96.7, 2008.7504
data 98.9, 2008.8336

end sub