' demo_lint2.bas March 8, 2017

' this program demonstrates the procedure for calling
' the subroutines interp2 which performs a bilinear
' interpolation of tabular data of the form z = f(x,y)
' input by the user.

' note: the data must be ordered as follows:

' x(1) < x(2) < x(3) < .... < x(n)
' y(1) < y(2) < y(3) < .... < y(n)

' this program demonstrates the use of interp2 for the
' interpolation of the function defined by

' z = f(x,y) = sin(x + y)

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

print " "

print "program demo_lint2"

print " "
print "this program demonstrates the procedure for calling the"
print "the subroutine interp2 which performs a bilinear"
print "interpolation of tabular data of the form z = f(x,y)"
print "input by the user."
print " "
print "note: the data must be ordered as follows:"
print " "
print "x(1) < x(2) < x(3) < .... < x(n)"
print "y(1) < y(2) < y(3) < .... < y(n)"
print " "
print "this program demonstrates the use of interp2 for the"
print "interpolation of the function defined by"
print " "
print " z = f(x,y) = sin(x + y)"

do

print " "

print "please input the number of x data points"

print "(a minimum of two x data points must be input)"

input nx%

loop until (nx% >= 2)

do

print " "

print "please input the number of y data points"

print "(a minimum of two y data points must be input)"

input ny%

loop until (ny% >= 2)

do

print " "

print "please input a value for delta-x"

input dx

loop until (dx > 0.0)

do
print

print "please input a value for delta-y"

input dy

loop until (dy > 0.0)

' dimension arrays

dim x(nx%) as float, y(ny%) as float, z(nx%, ny%) as float

print " "

print "please wait, computing data points ..."

' compute data points for the function defined by

' z = f(x, y) = sin(x + y)

' where x = dx * i and y = dy * j

' for i = 1, 2, ... nx and j = 1, 2, ... ny

for i% = 1 to nx%

x(i%) = i% * dx

for j% = 1 to ny%

y(j%) = j% * dy

z(i%, j%) = sin(x(i%) + y(j%))

next j%

next i%

do

print " "

print "please input the x value to interpolate"

print "(this must be a value between ", x(1), " and ", x(nx%)

input xval

loop until (xval >= x(1) and xval <= x(nx%))

do

print " "

print "please input the y value to interpolate"

print "(this must be a value between ", y(1), " and ", y(ny%)

input yval

loop until (yval >= y(1) and yval <= y(ny%))

' perform bilinear interpolation

interp2(nx%, ny%, xval, yval, fval)

print " "
print "program demo_lint2"
print " "
print "<linear interpolation of z = f(x,y)>"
print " "

print "number of x data points = ", nx%
print " "
print "number of y data points = ", ny%

print " "
print "x interpolated value = ", xval
print " "
print "y interpolated value = ", yval

print " "
print "function value = ", fval
print " "
print "interpolation error = ", abs(fval - sin(xval + yval))
print " "

end

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

sub interp2 (nx%, ny%, xval, yval, zval)

' bilinear interpolation subroutine

' z = f(x, y)

' input

' nx% = number of x data points (nx% >= 2)
' ny% = number of y data points (ny% >= 2)
' x() = vector of x data (nx% rows)
' y() = vector of y data (ny% rows)
' z() = array of z data (nx% rows by ny% columns)
' xval = x interpolation argument
' yval = y interpolation argument

' output

' zval = z function value at xval, yval

' note: the data must be ordered as follows:

' x(1) < x(2) < x(3) < .... < x(n)
' y(1) < y(2) < y(3) < .... < y(n)

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

local g1, g2, xinfac, yinfac

' compute x index and interpolation factor

for i% = 1 to nx%

if (xval <= x(i%)) then

if (i% = 1) then

xinfac = 0.0

indexx% = 1

else

j% = i% - 1

xinfac = (xval - x(j%)) / (x(i%) - x(j%))

indexx% = j%

end if

exit for

end if

next i%

' compute y index and interpolation factor

for i% = 1 to ny%

if (yval <= y(i%)) then

if (i% = 1) then

yinfac = 0.0

indexy% = 1

else

j% = i% - 1

yinfac = (yval - y(j%)) / (y(i%) - y(j%))

indexy% = j%

end if

exit for

end if

next i%

g1 = z(indexx%, indexy%) + xinfac * (z(indexx% + 1, indexy%) - z(indexx%, indexy%))

g2 = z(indexx%, indexy% + 1) + xinfac * (z(indexx% + 1, indexy% + 1) - z(indexx%, indexy% + 1))

zval = g1 + yinfac * (g2 - g1)

end sub