/***************************************************************
* Integration of a discrete function F(x) by Simpson's method  *
* ------------------------------------------------------------ *
* SAMPLE RUN                                                   *
*                                                              *
* (Find integral of digitalized function sin(x) from x=0 to    *
*  x=1).                                                       *
*                                                              *
* Integration of a discrete real function F(x)                 *
*             by Simpson's Method                              *
*                                                              *
* Number of points: 25                                         *
* Begin x         : 0                                          *
* End x           : 1                                          *
*                                                              *
* Value of integral: 0.459698                                  *
*                                                              *
*                                                              *
*                               C++ version by J-P Moreau.     *
*                                   (www.jpmoreau.fr)          *
***************************************************************/
#include <stdio.h>
#include <math.h>

#define  SIZE  100


double   Xi[SIZE], Yi[SIZE];
int      i, npoints;
double   dx,result,x,x1,x2;


/**********************************************************
* This function returns the integral value of a discrete  *
* real function F(x) defined by n points from x1 to xn.   *
* ------------------------------------------------------- *
* INPUTS:                                                 *
*          n      Number of points (xi, yi)               *
*          X      Table of n xi values in ascending order *
*          Y      Table of n yi values                    *
*                                                         *
* OUTPUT          The function returns the integral of    *
*                 F(x) from x1 to xn.                     *
* ------------------------------------------------------- *
* Ref.: "Mathematiques en Turbo Pascal By Alain Reverchon *
*        and Marc Ducamp, Armand Colin Editor, Paris,     *
*        1990" [BIBLI 14].                                *
**********************************************************/
double DiscreetIntegral(int n, double *X, double *Y)  {
  int i,i1,i2; double h,k,j,r;
  r=0;
  i=0;
  while (i<=n-3) {
    i1=i+1; i2=i+2; h=X[i1]-X[i];
    k=X[i2]-X[i1]; j=h+k;
    if (h==0) return 0;
    r += ((h+h-k)*Y[i]+(j*j*Y[i1]+h*(k+k-h)*Y[i2])/k)*j/6/h;
    i=i2;
  }
  if ((n%2)==0)
    r += (X[n-1]-X[n-2])*(Y[n-1]+Y[n-2])/2;
  return r;
}


void main()  {

  printf("    Integration of a discrete real function F(x)\n");
  printf("                by Simpson's Method\n\n");
  printf(" Number of points: "); scanf("%d", &npoints);
  printf(" Begin x         : "); scanf("%lf", &x1);
  printf(" End x           : "); scanf("%lf", &x2);
  // build up test discreet curve F(x)=sin(x)
  dx=(x2-x1)/(npoints-1);
  x=x1-dx;
  for (i=0; i<npoints; i++) {
    x += dx;
    Xi[i]=x; Yi[i]=sin(x);
  }

  // call integration function
  result = DiscreetIntegral(npoints,Xi,Yi);

  // print value of integral
  printf("\n\n Value of integral: %f\n\n", result);
    
}

// end of file disinteg.cpp