/**********
 * Numerical Integration
 **********/

#include <stdio.h>
#include <math.h>

/* here's the function we're going to integrate */
double f(double x){
  return( exp( -x*x ) );
}

/* perform a trapezoidal rule numerical integration
   on f() from xMin to xMax using Npoints.  Return the answer */
double TrapezoidIntegration (int Npoints, double xMin, double xMax) {
  double x, dx;
  double sum;
  int i, Nm1;
  sum = 0;
  Nm1 = Npoints - 1;
  dx = (xMax-xMin)/Nm1;
  for (i=1; i<Nm1; i++){ /* this sums over interior; not left or right end */
    x = xMin + i*dx;
    sum = sum + f(x);
  }
  return( dx* ( f(xMin)/2.0 + sum + f(xMax)/2.0 ) );
}

main () {
  double answer;

  int npoints;
  double xMin, xMax;

  printf(" Let's do some integration \n ");

  xMin = 0.0;
  xMax = 1000.0;
  npoints = 10000;

  
  xMax = 2000.0;
  answer = TrapezoidIntegration( npoints, xMin, xMax );
  printf(" with n = %i, xMax=%g ,  answer is %.16g \n", npoints,xMax,answer);

  xMax = 10000.0;
  answer = TrapezoidIntegration( npoints, xMin, xMax );
  printf(" with n = %i, xMax=%g ,  answer is %.16g \n", npoints,xMax,answer);

  xMax = 20000.0;
  answer = TrapezoidIntegration( npoints, xMin, xMax );
  printf(" with n = %i, xMax=%g ,  answer is %.16g \n", npoints,xMax,answer);

}