' program demo_sfit.bas March 30, 2017

' this program demonstrates the procedure for calling
' the subroutine surfit.bas which calculates fitting
' coefficients for 3-dimensional surfaces of the form
' z = f(x,y) from x, y and z data points input by the
' user. It also illustrates how to interpolate an x-y
' data pair provided by the user.

' this program demonstrates the use of surfit to find
' fitting coefficients for the surface defined by

' z = -cos(3*pi*x/10)*cos(pi*y/10)

' for -5 <= x <= +5 and -5 <= y <= +5

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

option default float

option base 1

' dimension data and coefficient arrays

dim xdata(121), ydata(121), zdata(121), coef(11, 11)

print " "
print "program demo_sfit"
print "-----------------"
print " "
print "this program demonstrates the procedure for calling"
print "the subroutine surfit.bas which calculates fitting"
print "coefficients for 3-dimensional surfaces of the form"
print "z = f(x,y) from x, y and z data points input by the"
print "user."
print " "
print "this program demonstrates the use of surfit to find"
print "fitting coefficients for the surface defined by"
print " "
print " z = -cos(3*pi*x/10)*cos(pi*y/10)"
print " "
print " for -5 <= x <= +5 and -5 <= y <= +5"

' generate data points for the surface defined by
' z = -cos(3*pi*x/10)*cos(pi*y/10)
' for -5 <= x <= +5 and -5 <= y <= +5

ndata% = 0

for i% = -5 to 5 step 1

for j% = -5 to 5 step 1

ndata% = ndata% + 1

x = i%

xdata(ndata%) = i%

y = j%

ydata(ndata%) = j%

z = -cos(3.0 * pi * x / 10.0) * cos(pi * y / 10.0)

zdata(ndata%) = z

next j%

next i%

do

print " "

print "please input the degree of the surface fit (degree <= 10)"

input ndeg%

loop until (ndeg% >= 0 and ndeg% <= 10)

print " "
print "please wait, computing coefficients";
print " "

surfit(ndeg%, ndata%, xdata(), ydata(), zdata(), coef())

' print fitting coefficients

for i% = 1 to ndeg% + 1

print " "
print "fitting coefficients"
print " "

i1% = i% + 1

for j% = 1 to i%

print "coef("; i1% - j% - 1; ","; j% - 1; ") =";

print coef(j%, i1% - j%)

next j%

next i%

do

print " "
print "would you like to interpolate an x, y data point"
print "(y = yes, n = no)"

input fit$

loop until (fit$ = "y" or fit$ = "n")

if (fit$ = "y") then

do
print
print "please input the x data point (-5 <= x <= +5)"
input x

loop until (abs(x) <= 5.0)

do

print
print "please input the y data point (-5 <= y <= +5)"
input y

loop until (abs(y) <= 5.0)

' perform interpolation

sfunction(ndeg%, coef(), x, y, z)

' print results of interpolation

print " "
print "x data point = ", x
print " "
print "y data point = ", y
print " "
print "z data point = ", z
print " "

' analytic z-value

za = -cos(3.0 * pi * x / 10.0) * cos(pi * y / 10.0)

print "za = ", za

print "z_analytic - z_fit = ", za - z
print " "

end if

end

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

sub surfit (ndeg%, ndata%, xdata(), ydata(), zdata(), coef())

' surface fit subroutine

' input

' ndeg% = degree of surface fit (ndeg% <= 10)
' ndata% = number of x, y and z data points
' xdata() = vector of x data values (ndata% rows)
' ydata() = vector of y data values (ndata% rows)
' zdata() = vector of z data values (ndata% rows)

' output

' coef() = array of surface fit coefficients (11 rows by 11 columns)

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

local ndeg1% = ndeg% + 1

local ndeg2% = ndeg% + 2

local nrow% = ndeg1% * ndeg2% / 2

local ncol% = nrow% + 1

local nalf% = 2 * ndeg%

local nalf1% = nalf% + 1

local nalf2% = nalf% + 2

local xa(ndeg1% * ndeg2% / 2)

local zadd(ndeg1%, ndeg1%)

local hsmg(2 * ndeg1%, 2 * ndeg1%)

local xmatrix(ndeg1% * ndeg2% / 2, ndeg1% * ndeg2% / 2 + 1)

for i% = 1 to nalf1%

for j% = 1 to nalf2% - i%

hsmg(i%, j%) = 0.0

next j%

next i%

for i% = 1 to ndeg1%

for j% = 1 to ndeg2% - i%

zadd(i%, j%) = 0.0

next j%

next i%

z1 = 0.0

for i% = 1 to ndata%

z1 = z1 + zdata(i%)

next i%

if (ndeg% <> 0) then

for i% = 1 to ndata%

a = 1.0

xa(1) = 1.0

x = xdata(i%)

y = ydata(i%)

z = zdata(i%)

for j% = 2 to nalf1%

a = a * x

xa(j%) = a

if (j% <= ndeg1%) then zadd(1, j%) = zadd(1, j%) + a * z

hsmg(1, j%) = hsmg(1, j%) + a

next j%

for l% = 2 to nalf1%

for m% = 1 to nalf2% - l%

xa(m%) = xa(m%) * y

hsmg(l%, m%) = hsmg(l%, m%) + xa(m%)

next m%

if (l% <= ndeg1%) then

for n% = 1 to ndeg2% - l%

zadd(l%, n%) = zadd(l%, n%) + xa(n%) * z

next n%

end if

next l%

next i%

end if

hsmg(1, 1) = ndata%

zadd(1, 1) = z1

m% = 0

for i% = 1 to ndeg1%

for j% = 1 to ndeg2% - i%

n% = 0

m% = m% + 1

xmatrix(m%, ncol%) = zadd(i%, j%)

for k% = i% to ndeg% + i%

for l% = j% to ndeg% + i% + j% - k%

n% = n% + 1

xmatrix(m%, n%) = hsmg(k%, l%)

next l%

next k%

next j%

next i%

if (ndeg% = 0) then

xa(1) = xmatrix(1, 2) / xmatrix(1, 1)

return

end if

nr1% = nrow% - 1

for k% = 1 to nr1%

nb% = 0

amax = xmatrix(k%, k%)

i1% = k% + 1

for ki% = i1% to nrow%

if (abs(xmatrix(ki%, k%)) > abs(amax)) then

amax = xmatrix(ki%, k%)

nb% = ki%

end if

next ki%

if (nb% <> 0) then

for nl% = k% to ncol%

xsave = xmatrix(nb%, nl%)

xmatrix(nb%, nl%) = xmatrix(k%, nl%)

xmatrix(k%, nl%) = xsave

next nl%

end if

const1 = xmatrix(k%, k%)

for i% = i1% to nrow%

xa(i%) = -xmatrix(i%, k%) / const1

next i%

for j% = i1% to ncol%

const1 = xmatrix(k%, j%)

for i% = i1% to nrow%

xmatrix(i%, j%) = xmatrix(i%, j%) + const1 * xa(i%)

next i%

next j%

next k%

xa(nrow%) = xmatrix(nrow%, ncol%) / xmatrix(nrow%, nrow%)

for i% = 1 to nr1%

k% = nrow% - i%

sum = 0.0

const1 = xmatrix(k%, k%)

for j% = k% + 1 to i% + k%

sum = sum + xmatrix(k%, j%) * xa(j%)

next j%

xa(k%) = (xmatrix(k%, ncol%) - sum) / const1

next i%

' load coef array

k% = 0

for i% = 1 to ndeg1%

for j% = 1 to ndeg2% - i%

k% = k% + 1

coef(i%, j%) = xa(k%)

next j%

next i%

end sub

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

sub sfunction(ndeg%, coef(), x, y, z)

' interpolate surface function subroutine

' z = f(x,y)

' input

' ndeg% = degree of surface fit
' coef() = array of surface fit coefficents
' (up to 11 rows by 11 columns)
' x = x data value
' y = y data value

' output

' z = interpolated z data value of surface at x, y

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

local a, y1

local ndeg1% = ndeg% + 1

local ndeg2% = ndeg% + 2

local work(ndeg1% * ndeg2% / 2)

if (ndeg% = 0) then

z = coef(1, 1)

exit sub

end if

a = 1.0

y1 = 0.0

work(1) = 1.0

for i% = 2 to ndeg1%

a = a * x

work(i%) = a

y1 = y1 + a * coef(1, i%)

next i%

for i% = 2 to ndeg1%

for j% = 1 to ndeg2% - i%

work(j%) = work(j%) * y

y1 = y1 + work(j%) * coef(i%, j%)

next j%

next i%

z = y1 + coef(1, 1)

end sub