/***************************************************
* 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 <stdio.h>
#include <math.h>

#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[1] || xx>=X[iv-3]) {
    n=0 ; return;
  }
  X[0]=2.0*X[1]-X[2];
  //Calculate Akima coefficients, a and b
  for (i=1; i<iv; i++)
    //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; i<iv+1; i++)  {
    a=fabs(XM[i+3]-XM[i+2]);
    b=fabs(XM[i+1]-XM[i]);
    if (a+b!=0) goto e100;
    Z[i]=(XM[i+2]+XM[i+1])/2.0;
    goto e200;
e100: Z[i]=(a*XM[i+1]+b*XM[i+2])/(a+b);
e200: ;}
  //Find relevant table interval
  i=0;
e300: i++;
  if (xx>X[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 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). *
*******************************************************/
void Integrale()  {
// Labels: e100,e200
  double xi1,xi2,x3; 
  zz=0 ; z1=0;
  // Check to see if end points are correct
  if (x1<X[1]) return;
  if (x2>X[iv-3]) return;
  // If x1>x2 then switch and set flag
  if (x1<x2) goto e100;
  x3=x1 ; x1=x2 ; x2=x3 ; z1=1;
e100: if (x2==x1) return;
  // Start trapezoidal integrations
  // First integration to get xi1 
  Trapez1(&xi1);
  // Second round to get xi2
  Trapez2(&xi2);
  // Richardson extrapolation
  zz=(4*xi2-xi1)/3;
  // Check to see if the end points have been reversed
  if (z1==0) goto e200;
  zz=-zz ; x2=x1 ; x1=x3;
  // Reset error flag
e200: z1=1; 
  // zz is the integral desired
  return;
}


/*****************************************************
* Routine for the first trapezoidal integration, xi1 *
* n1 keeps track of the number of intervals.         *
*****************************************************/
void Trapez1(double *xi1) {
//Labels: e100,e150,e200
  double d; int j1;
  *xi1=0 ; n1 = 0 ; xx=x1;
  // Call interpolation subroutine
  Interpol_Akima();
  // Is there at least one table interval?
  if (x2>X[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 (x2<X[j1+3]) goto e200;
  // Otherwise, keep summing
  n1++;
  *xi1=*xi1+(Y[j1+1]+Y[j1+3])*(X[j1+3]-X[j1+1])/2.0;
  j1=j1+2;
  goto e150;
e200: xx=x2;
  // Find last yy value
  Interpol_Akima();
  *xi1=*xi1+(yy+Y[j1+1])*(x2-X[j1+1])/2.0;
  n1++;
  return;
}

// ***** integration for xi2 *****
void Trapez2(double *xi2)  {
// Labels: e600,e650,e700
  double d; int j1;
  *xi2=0 ; xx=x1;
  Interpol_Akima();
  d=yy;
  if (x2>X[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<X[j1+2]) goto e700;
  *xi2=*xi2+(Y[j1+1]+Y[j1+2])*(X[j1+2]-X[j1+1])/2.0;
  j1=j1+1;
  goto e650;
e700: xx=x2-(x2-X[j1+1])/2.0;
  Interpol_Akima();
  d=yy;
  *xi2=*xi2+(Y[j1+1]+d)*(x2-xx)/2.0;
  xx=x2;
  Interpol_Akima();
  *xi2=*xi2+(d+yy)*(x2-X[j1+1])/4.0;
  return;
}


void main()  {
  iv=14;   //Number of pooints in table}
  printf("\n Integration of SIN(X)\n\n");
  /* Input sine table
  -----------------------------------------------------------------
   Sine table values from  Handbook of mathematical functions
   by M. Abramowitz and I.A. Stegun, NBS, june 1964
  ---------------------------------------------------------------*/
  X[1]=0.000; Y[1]=0.00000000; X[2]=0.125; Y[2]=0.12467473;
  X[3]=0.217; Y[3]=0.21530095; X[4]=0.299; Y[4]=0.29456472;
  X[5]=0.376; Y[5]=0.36720285; X[6]=0.450; Y[6]=0.43496553;
  X[7]=0.520; Y[7]=0.49688014; X[8]=0.589; Y[8]=0.55552980;
  X[9]=0.656; Y[9]=0.60995199; X[10]=0.721; Y[10]=0.66013615;
  X[11]=0.7853981634; Y[11]=0.7071067812;
  X[12]=0.849; Y[12]=0.75062005; X[13]=0.911; Y[13]=0.79011709;
  X[14]=0.972; Y[14]=0.82601466;
  //---------------------------------------------------------------
  printf(" Input start point  end point: "); scanf("%lf %lf",&x1,&x2);
  
  Integrale();  // Call integral function

  printf("\n Integral from %6.3f to %6.3f equals %10.7f\n",x1,x2,zz);
  printf(" Error check: %d\n\n",z1);
}

// End of file Integra.cpp