'===========================================================================
' Subject: Parabolic Interpolation            Date: 10-19-03 (  :  )       
'  Author: David Williams                     Code: Qbasic,QB,PDS          
'  Origin: david.williams@ablelink.org      Packet: ALGOR.ABC
'===========================================================================
' PARINTER.BAS

' David Williams, 2003
' david.williams@ablelink.org

' This version dated 2003 Aug. 07

' Program demonstrates parabolic (quadratic) interpolation.

DECLARE SUB Description ()
DECLARE FUNCTION Interp# (A#)
 
DEFDBL A-Z
 
CALL Description

DIM SHARED SS  ' step size in array
 
PRINT
 
DO
  INPUT "Step size (1 - 20)"; SS
LOOP UNTIL SS >= 1 AND SS <= 20
 
L% = INT(90 / SS) + 1 ' L% defines the required size of the array
IF SS * (L% - .5) <= 90 THEN L% = L% + 1
 
DIM SHARED S(-1 TO L%)  ' set up look-up array
 
F = SS * ATN(1) / 45
 
FOR X = -1 TO L%
 S(X) = SIN(X * F)
NEXT
 
Start:
PRINT
 
DO
  INPUT "Angle (degrees, 0 - 90 [-1 to end])"; A
  IF A < 0 THEN END
LOOP UNTIL A <= 90
 
PRINT
 
G = A * F / SS
S! = SIN(G)
PRINT "True sine:"; TAB(35); S!
 
B = INT(A / SS)
C = A / SS - B
L! = S(B) + (S(B + 1) - S(B)) * C
PRINT "Linearly interpolated sine:"; TAB(35); L!

P! = Interp(A)
PRINT "Parabolically interpolated sine:"; TAB(35); P!
 
PRINT
LE! = ABS(S! - L!)
PE! = ABS(S! - P!)
PRINT "Linear error:"; TAB(35); LE!
PRINT "Parabolic error:"; TAB(35); PE!
IF PE! > 0 THEN
  PRINT "Error ratio (Linear:Parabolic):"; TAB(35); STR$(LE! / PE!); ":1"
END IF

GOTO Start

SUB Description
  CLS
  PRINT "The process of interpolation is usually done by assuming"
  PRINT "that adjacent known points are joined by straight lines,"
  PRINT "the equations of which can be used to estimate unknown"
  PRINT "quantities. However, this process takes no account of the"
  PRINT "curvature of the graph. Better results can usually be obtained"
  PRINT "by fitting a parabola to the THREE known points that are closest"
  PRINT "to the one that is to be interpolated, and using the equation"
  PRINT "of this parabola to estimate the unknown coordinate."
  PRINT
  PRINT "This program demonstrates a function that does this. It assumes"
  PRINT "that an array has been filled with the ordinates (Y-coordinates)"
  PRINT "of a set of known points whose abscissas (X-coordinates)"
  PRINT "are equally spaced (the 'step size'). The first two lines of"
  PRINT "the function calculate the array index and abscissa of the"
  PRINT "known point whose abscissa is closest to that of the interpolated"
  PRINT "point. The rest of the function does the interpolation."
  PRINT
  PRINT "For demonstration purposes, a sine function is used. The angle"
  PRINT "is entered by the user, in degrees."
  PRINT
END SUB

FUNCTION Interp (A)
  ' Performs parabolic interpolation among elements of an array.
  ' The array S() and the variable SS are shared. SS is the "step size",
  ' the difference between the abscissas of consecutive array elements.
  ' A is the abscissa of the point to be interpolated.
  ' ("Abscissa" means X-coordinate of any point (X, Y).)
  ' Array contains ordinates (Y-coordinates) of known points.
  ' Interp calculates ordinate of interpolated point, by fitting a
  ' parabola to the three closest known points.
 
  N = INT(A / SS + .5)
  ' N is index of array element with abscissa closest to interpolated point.
  ' S(N) becomes middle one of the three points used.
  B = SS * N ' Abscissa of S(N)
 
  ' If the array is set up differently than this demonstration one,
  ' the above two lines will have to be adjusted accordingly.

  S1 = S(N - 1): S2 = S(N): S3 = S(N + 1)  ' 3 nearest array entries
 
  C1 = (S1 - 2 * S2 + S3) / (2 * SS * SS)  ' )
  C2 = (S3 - S1) / (2 * SS) - 2 * C1 * B   ' ) 3 quadratic coefficients
  C3 = S2 - B * (C1 * B + C2)              ' )
 
  Interp = C1 * A * A + C2 * A + C3 ' Quadratic equation does interpolation.
END FUNCTION