/********************************/
/*    Integration               */
/*   Attempting to integrate    */
/*     exp(-x*x) over all space */
/*                              */
/*   Jon Franklin  09-24-2001   */
/********************************/

#include <stdio.h>
#include <math.h>
#define MAX 99999


/* function testf
     inputs:   h   = step size
               f[] = array of function values
     outputs:  returns sum(f)*h

   This function also prints lots of crap to the terminal.
   You can turn on debugging by uncommenting the DEBUG lines below.
*/
double testf(double h, double f[])
{ 
  int i;
  double x;

  x=0.0;

  /* DEBUGGING */
  /* printf(" x=%g \n",x); */

  for(i=0;i<MAX;i++)
  {
    x = x + h*f[i];
  }
  return(x);
}

double Simp(double h,double f[])
{
  double s0,s1,s2;
  int i;
  s0 = 0;
  s1 = 0;
  s2 = 0;
  for (i=1; i<MAX-1;i=i+2)
  {
    s1=s1+f[i-1];
    s0=s0+f[i];
    s2=s2+f[i+1];
  }
/*  if (MAX%2 == 0) */
/*return((h*(s1+4*s0+s2)/3)+(h*(5*f[MAX-1]+8*f[MAX-2]-f[MAX-3])/12));*/

  printf("\n\ns0 = %lf, s1 = %lf, s2 = %lf\n\n",s0,s1,s2);

  return(h*(s1+4*s0+s2)/3);
}

double fun(double u)
{
  return( exp(-u*u) );
}

main()
{
  int l,n; 
  double dx,pi,u,S;
  double x[MAX];
  double func[MAX];

  pi = 4.0*atan(1.0);
  printf("pi = %lf\n",pi);
  dx = pi/2/MAX;
  printf("dx = %lf\n",dx);
  for (l=0; l<MAX; l++)
  {
    x[l] = dx*l;
/*    printf("x[%i]=%lf    ",l,x[l]);  */
    u = tan(x[l]);    
    func[l] = fun(u);
/*    printf("The value of e^-%lf^2 is  %lf\n",u,func[l]); */
 }

  S= Simp(dx,func);
  printf("\n\n%16lf\n\n",2*S);

  S=testf(dx,func);
  printf("\n\n%16lf\n\n",2*S);

}