/*************************************************** * 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].* * * * C++ 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 * * * ***************************************************/ #include #include #define SIZE 15 double X[SIZE],Y[SIZE],XM[SIZE+3],Z[SIZE]; double a,b,x1,x2,xx,yy,zz; int i,iv,n,n1,z1; /******************************************************* * 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. * *******************************************************/ void Interpol_Akima() { //Labels: 100,200,300 int i; n=1; //special case xx=0 if (xx==0.0) { yy=0.0; return; } //Check to see if interpolation point is correct if (xx=X[iv-3]) { n=0 ; return; } X[0]=2.0*X[1]-X[2]; //Calculate Akima coefficients, a and b for (i=1; iX[i]) goto e300; i--; //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); return; } void Trapez1(double *); void Trapez2(double *); /******************************************************* * 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). * *******************************************************/ void Integrale() { // Labels: e100,e200 double xi1,xi2,x3; zz=0 ; z1=0; // Check to see if end points are correct if (x1X[iv-3]) return; // If x1>x2 then switch and set flag if (x1X[i+1]) goto e100; //If not, integral is simple n1++ ; d=yy ; xx=x2; // Find end point yy value Interpol_Akima(); *xi1=(yy+d)*(x2-x1)/2.0; return; // At least one table interval must be summed over e100: j1=i; *xi1=*xi1+(yy+Y[i+1])*(X[i+1]-xx)/2.0; // Any more intervals? If not, finish integral with end point e150: if (x2X[i+1]) goto e600; xx=x1+(x2-x1)/2.0; Interpol_Akima(); *xi2=*xi2+(d+yy)*(x2-x1)/4.0; return; e600: 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; e650: if (x2