'********************************************************
'*   Program to demonstrate the Akima spline fitting    *
'*       of Function SIN(X) in double precision         *
'* ---------------------------------------------------- *
'*   Reference: BASIC Scientific Subroutines, Vol. II   *
'*   by F.R. Ruckdeschel, BYTE/McGRAWW-HILL, 1981 [1].  *
'*                                                      *
'* ---------------------------------------------------- *
'* SAMPLE RUN:                                          *
'*                                                      *
'* Akima spline fitting of SIN(X):                      *
'*                                                      *
'*   X   SIN(X) HANDBOOK   AKIMA INTERPOLATION   ERROR  *
'* ---------------------------------------------------- *
'* 0.00    0.0000000          0.0000000       0.0000000 *
'* 0.05    0.0499792          0.0500402      -0.0000610 *
'* 0.10    0.0998334          0.0998435      -0.0000101 *
'* 0.15    0.1494381          0.1494310       0.0000072 *
'* 0.20    0.1986693          0.1986459       0.0000235 *
'* 0.25    0.2474040          0.2474157      -0.0000118 *
'* 0.30    0.2955202          0.2955218      -0.0000016 *
'* 0.35    0.3428978          0.3428916       0.0000062 *
'* 0.40    0.3894183          0.3894265      -0.0000081 *
'* 0.45    0.4349655          0.4349655       0.0000000 *
'* 0.50    0.4794255          0.4794204       0.0000051 *
'* 0.55    0.5226872          0.5226894      -0.0000021 *
'* 0.60    0.5646425          0.5646493      -0.0000068 *
'* 0.65    0.6051864          0.6051821       0.0000043 *
'* 0.70    0.6442177          0.6442141       0.0000035 *
'* 0.75    0.6816388          0.6816405      -0.0000017 *
'* ---------------------------------------------------- *
'*                                                      *
'********************************************************
defint i-n
defdbl a-h,o-z
  iv=14  'Number of pooints in table
  dim X(iv+1),Y(iv+1),XM(iv+4),Z(iv+1)
  cls
  print
  print " Akima spline fitting of SIN(X):"
  'Input table
  for i=1 to iv
    read X(i), Y(i)
  next i
  '------------------------------------------------------------
  'Sine table values from  Handbook of mathematical
  'functions by M. Abramowitz and I.A. Stegun,
  'NBS, june 1964
  data 0.000#,0.00000000#,0.125#,0.12467473#
  data 0.217#,0.21530095#,0.299#,0.29456472#,0.376#,0.36720285#
  data 0.450#,0.43496553#,0.520#,0.49688014#,0.589#,0.55552980#
  data 0.656#,0.60995199#,0.721#,0.66013615#
  data 0.7853981634#,0.7071067812#
  data 0.849#,0.75062005#,0.911#,0.79011709#
  data 0.972#,0.82601466#
  '------------------------------------------------------------
  'print header
  print
  print " X   SIN(X) HANDBOOK   AKIMA INTERPOLATION   ERROR "
  print "----------------------------------------------------"
  'main loop
  for xx=0 to 0.80 step 0.05
    gosub 1000
    print using "#.##    #.#######          #.#######      ##.#######"; xx; SIN(xx); yy; SIN(xx)-yy
  next xx
  'print footer
  print "----------------------------------------------------"
  print
end


'********************************************************
'*          Akima spline fitting subroutine             *
'* ---------------------------------------------------- *
'* The input table is X(i), Y(i), where Y(i) is the     *
'* dependant variable. The interpolation point is xx,   *
'* which is assumed to be in the interval of the table  *
'* with at least one table value to the left, and three *
'* to the right. The interpolated returned value is yy. *
'* n is returned as an error check (n=0 implies error). *
'* It is also assumed that the X(i) are in ascending    *
'* order.                                               *
'********************************************************
1000 'Check to see if interpolation point is correct
  n=1
  if xx>=X(1) then goto 1100
  n=0 : return
1100 if xx<=X(iv-3) then goto 1200
  n=0 : return
1200 X(0)=2#*X(1)-X(2)
  'Calculate Akima coefficients
  for i=1 to iv-1
    'Shift i to i+2
    XM(i+2)=(Y(i+1)-Y(i))/(X(i+1)-X(i))
  next i
  XM(iv+2)=2#*XM(iv+1)-XM(iv)
  XM(iv+3)=2#*XM(iv+2)-XM(iv+1)
  XM(2)=2#*Xm(3)-XM(4)
  XM(1)=2#*XM(2)-XM(3)
  for i=1 to iv
    a=ABS(XM(i+3)-XM(i+2))
    b=ABS(XM(i+1)-XM(i))
    if a+b<>0 then goto 1300
    Z(i)=(XM(i+2)+XM(i+1))/2#
    goto 1350
1300 Z(i)=(a*XM(i+1)+b*XM(i+2))/(a+b)
1350 next i
  'Find relevant table interval
  i=0
1400 i=i+1
  if xx>X(i) then goto 1400
  i=i-1
  'Begin interpolation
  b=X(i+1)-X(i)
  a=xx-X(i)
  yy=Y(i)+Z(i)*a+(3#*XM(i+2)-2#*Z(i)-Z(i+1))*a*a/b
  yy=yy+(Z(i)+Z(i+1)-2#*XM(i+2))*a*a*a/(b*b)
return

'End of file Akima.bas