{*************************************************** * 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].* * * * Pascal version by J-P Moreau, Paris * * (www.jpmoreau.fr) * * ------------------------------------------------ * * 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 * * * ***************************************************} PROGRAM Integra; Uses WinCrt; CONST SIZE = 14; VAR X : Array[0..SIZE] of double; Y : Array[0..SIZE] of double; XM : Array[0..SIZE+3] of double; Z : Array[0..SIZE] of double; a,b,x1,x2,xx,yy,zz : double; i,iv,n,n1,z1 : INTEGER; {******************************************************* * 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. * *******************************************************} PROCEDURE Interpol_Akima; Label 100,200,300,fin; Var i:integer; Begin n:=1; {special case xx:=0} if xx=0.0 then begin yy:=0.0; goto fin end; {Check to see if interpolation point is correct} if (xx=X[iv-3]) then begin n:=0 ; goto fin end; X[0]:=2.0*X[1]-X[2]; {Calculate Akima coefficients, a and b} for i:=1 to iv-1 do {Shift i to i+2} XM[i+2]:=(Y[i+1]-Y[i])/(X[i+1]-X[i]); XM[iv+2]:=2.0*XM[iv+1]-XM[iv]; XM[iv+3]:=2.0*XM[iv+2]-XM[iv+1]; XM[2]:=2.0*Xm[3]-XM[4]; XM[1]:=2.0*XM[2]-XM[3]; for i:=1 to iv do begin a:=ABS(XM[i+3]-XM[i+2]); b:=ABS(XM[i+1]-XM[i]); if a+b<>0 then goto 100; Z[i]:=(XM[i+2]+XM[i+1])/2.0; goto 200; 100: Z[i]:=(a*XM[i+1]+b*XM[i+2])/(a+b); 200: end; {Find relevant table interval} i:=0; 300: i:=i+1; if xx>X[i] then goto 300; i:=i-1; {Begin interpolation} b:=X[i+1]-X[i]; a:=xx-X[i]; yy:=Y[i]+Z[i]*a+(3.0*XM[i+2]-2.0*Z[i]-Z[i+1])*a*a/b; yy:=yy+(Z[i]+Z[i+1]-2.0*XM[i+2])*a*a*a/(b*b); fin: End; PROCEDURE Trapez1(VAR xi1:double); FORWARD; PROCEDURE Trapez2(VAR xi2:double); FORWARD; {******************************************************* * 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 x1error). * *******************************************************} PROCEDURE Integrale; Label 100,200,fin; Var xi1,xi2,x3 : double; Begin zz:=0 ; z1:=0; {Check to see if end points are correct} if x1X[iv-3] then goto fin; {If x1>x2 then switch and set flag} if x1X[i+1] then goto 100; {If not, integral is simple} n1:=n1+1 ; d:=yy ; xx:=x2; {Find end point yy value} Interpol_Akima; xi1:=(yy+d)*(x2-x1)/2.0; goto fin; {At least one table interval must be summed over} 100: j1:=i; xi1:=xi1+(yy+Y[i+1])*(X[i+1]-xx)/2.0; {Any more intervals? If not, finish integral with end point} 150: if x2X[i+1] then goto 600; xx:=x1+(x2-x1)/2.0; Interpol_Akima; xi2:=xi2+(d+yy)*(x2-x1)/4.0; goto fin; 600: xx:=x1+(X[i+1]-x1)/2.0; j1:=i; Interpol_Akima; xi2:=xi2+(yy+d)*(xx-x1)/2.0; xi2:=xi2+(yy+Y[j1+1])*(X[j1+1]-xx)/2.0; 650: if x2