'****************************************************
'* Program to demonstrate the general integration   *
'* subroutine. Example is the integral of SIN(X)    *
'* from X1 to X2 in double precision.               *
'* ------------------------------------------------ *
'* Reference: BASIC Scientific Subroutines, Vol. II *
'* By F.R. Ruckdeschel, BYTE/McGRAWW-HILL, 1981 [1].*
'*                                                  *
'* ------------------------------------------------ *
'* SAMPLE RUN:                                      *
'*                                                  *
'*  Integration of SIN(X)                           *
'*                                                  *
'*  Input start point  end point: 0, 0.5            *
'*  Integral from 0.000 to 0.500 equals  0.1224231  *
'*  Error check: 1                                  *
'*                                                  * 
'****************************************************
defint i-n
defdbl a-h,o-z
  iv=14  'Number of pooints in table
  dim X(iv),Y(iv),XM(iv+3),Z(iv)
  cls
  print
  print " Integration of SIN(X)"
  print
  '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#
  '------------------------------------------------------------
  Input " Input start point, end point : ",x1, x2
  print
  gosub 2000
  print using " Integral from ##.### to ##.### equals ##.#######"; x1; x2; zz
  print using " Error check: #"; z1
  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

'********************************************************
'*       General integration subroutine (ITEG)          *
'* Interpolation by Akima (or other). Integration by    *
'* enhanced trapezoidal rule with Richardson extrapola- *
'* tion to give cubic accuracy.  The integration range  *
'* is (x1,x2). It is assumed that x1<x2, and that there *
'* is at least one table value to the left of x1, and   *
'* three to the right of x2. The result is returned in  *
'* zz. An error check is returned in z1 (z1=0-->error). *
'********************************************************
2000 zz=0 : z1=0
  'Check to see if end points are correct
  if x1<X(1) then return
  if x2>X(iv-3) then return
  'If x1>x2 then switch and set flag
  if x1<x2 then goto 2100
  x3=x1 : x1=x2 : x2=x3 : z1=1
2100 if x2=x1 then return
  'Start trapezoidal integrations
  'First integration to get xi1
  gosub 3000
  'Second round to get xi2
  gosub 3500
  'Richardson extrapolation
  zz=4#*xi2/3#-xi1/3#
  'Check to see if the end points have been reversed
  if z1=0 then goto 2200
  zz=-zz : x2=x1 : x1=x3
  'Reset error flag
2200 z1=1 : return
  'zz is the integral desired
return

'******************************************************
'* Routine for the first trapezoidal integration, xi1 *
'* n1 keeps track of the number of intervals.         *
'******************************************************
3000 xi1=0 : n1 = 0 : xx=x1
  'Call interpolation subroutine
  gosub 1000
  'Is there at least one table interval?
  if x2>X(i+1) then goto 3100
  'If not, integral is simple
  n1=n1+1 : d=yy : xx=x2
  'Find end point yy value
  gosub 1000
  xi1=(yy+d)*(x2-x1)/2#
  return
  'At least one table interval must be summed over
3100 j1=i
  xi1=xi1+(yy+Y(i+1))*(X(i+1)-xx)/2#
  'Any more intervals? If not, finish integral with end point
3150 if x2<X(j1+3) then goto 3200
  'Otherwise, keep summing
  n1=n1+1
  xi1=xi1+(Y(j1+1)+Y(j1+3))*(X(j1+3)-X(j1+1))/2#
  j1=j1+2
  goto 3150
3200 xx=x2
  'Find last yy value
  gosub 1000
  xi1=xi1+(yy+Y(j1+1))*(x2-X(j1+1))/2#
  n1=n1+1
return

'***** integration for xi2 *****
3500 xi2=0 : xx=x1
  gosub 1000
  d=yy
  if x2>X(i+1) then goto 3600
  xx=x1+(x2-x1)/2#
  gosub 1000
  xi2=xi2+(d+yy)*(x2-x1)/4#
  return
3600 xx=x1+(X(i+1)-x1)/2#
  j1=i
  gosub 1000
  xi2=xi2+(yy+d)*(xx-x1)/2#
  xi2=xi2+(yy+T(j1+1))*(X(j1+1)-xx)/2#
3650 if x2<X(j1+2) then goto 3700
  xi2=xi2+(Y(j1+1)+Y(j1+2))*(X(j1+2)-X(j1+1))/2#
  j1=j1+1
  goto 3650
3700 xx=x2-(x2-X(j1+1))/2#
  gosub 1000
  d=yy
  xi2=xi2+(Y(j1+1)+d)*(x2-xx)/2#
  xx=x2
  gosub 1000
  xi2=xi2+(d+yy)*(x2-X(j1+1))/4#
return

'End of file Integra.bas