/* 
************************************************************************
*  integ.cpp:  Integration using trapezoid, Simpson and Gauss rules    *
*  UNIX (DEC OSF, IBM AIX): cc integ.c -lm                             *
*    			                                                       *
*  comment: The derivation from the theoretical result for each method *
*           is saved in x y1 y2 y3 format.			                   * 
*  C.A. Bertulani - 03/20/2000
************************************************************************
*/

#include <stdio.h>
#include <math.h>
#include <iostream.h>
#include <iomanip.h>
#include <fstream.h>

typedef double Number;

/************************************************************************
/* function "gauss.cpp"  /* routine returns Legendre points, weights */

void gauss(int npts, int job, Number a, Number b, Number x[], Number w[])
{    
/*     npts     number of points                                       */
/*     job = 0  rescaling uniformly between (a,b)                      */
/*           1  for integral (0,b) with 50% points inside (0, ab/(a+b))*/
/*           2  for integral (a,inf) with 50% inside (a,b+2a)          */
/*     x, w     output grid points and weights.                        */ 

	int     m, i, j; 
      Number  t, t1, pp, p1, p2, p3;
      Number  pi = 3.1415926535897932385E0;
      Number  eps = 3.e-10;			/* limit for accuracy */
      
      m = (npts+1)/2;
      for(i=1; i<=m; i++)
      {  
         t  = cos(pi*(i-0.25)/(npts+0.5));
         t1 = 1;
         while((fabs(t-t1))>=eps)
         { 
            p1 = 1.0;
            p2 = 0.0;
            for(j=1; j<=npts; j++)
            {
               p3 = p2;
               p2 = p1;
               p1 = ((2*j-1)*t*p2-(j-1)*p3)/j;
	    }
            pp = npts*(t*p1-p2)/(t*t-1);
            t1 = t;
            t  = t1 - p1/pp;
         }   
         x[i-1] = -t;
         x[npts-i] = t;
         w[i-1]    = 2.0/((1-t*t)*pp*pp);
         w[npts-i] = w[i-1];
      } 

      if(job==0) 
      {
         for(i=0; i<npts ; i++)
         {
            x[i] = x[i]*(b-a)/2.0+(b+a)/2.0;
            w[i] = w[i]*(b-a)/2.0;
         }
      }
      if(job==1) 
      {
         for(i=0; i<npts; i++)
         {
            x[i] = a*b*(1+x[i]) / (b+a-(b-a)*x[i]);
            w[i] = w[i]*2*a*b*b /((b+a-(b-a)*x[i])*(b+a-(b-a)*x[i]));
         }
      }
      if(job==2) 
      {
         for(i=0; i<npts; i++)
         {
            x[i] = (b*x[i]+b+a+a) / (1-x[i]);
            w[i] =  w[i]*2*(a+b)  /((1-x[i])*(1-x[i]));
         }
      }
}
/* end function "gauss.cpp"
/************************************************************************

/* start main soubroutine
/************************************************************************
*/
#define max_in  101                     /* max number of intervals */
#define vmin 0.0                        /* ranges of integration */
#define vmax 1.0			
#define ME 2.7182818284590452354E0      /* Euler's number */

void main()
{
   int i;			 
   Number result;
   Number f(Number);		        /* subroutine containing */   
                             		/* function to integrate */
   Number trapez  (int no, Number min, Number max);    /* trapezoid rule */
   Number simpson (int no, Number min, Number max);    /* Simpson's rule */
   Number gaussint(int no, Number min, Number max);    /* Gauss' rule */
   Number romberg (int no, Number (*func)(Number), Number min, Number max);
	   /* romberg integration */

/* output */

   ofstream out("integ.dat");

   for (i=3; i<=max_in; i+=2)	    	/* Simpson's rule requires */
   { 					/* odd number of intervals */
      result = trapez(i, vmin, vmax);
	  out << setw(3) << i << setw(15) << setprecision(4) << fabs(result-1+1/ME);
      
      result = simpson(i, vmin, vmax);
	  out << setw(13) << setprecision(3) << fabs(result-1+1/ME);
      
      result = gaussint(i, vmin, vmax);
	  out << setw(13) << setprecision(3) << fabs(result-1+1/ME);
      
      result = romberg(i, f, vmin, vmax);
	  out << setw(13) << setprecision(3) << fabs(result-1+1/ME) << endl;
   }
   cout << "data stored in integ.dat\n" << endl;
}
/*------------------------end of main program-------------------------*/
/************************************************************************

/* the function we want to integrate */
Number f (Number x)
{ 
   return (exp(-x));               
}

/************************************************************************
/* Integration using trapezoid rule */
Number trapez (int no, Number min, Number max)
{
   int n;				
   Number interval, sum=0., x;		 

   interval = ((max-min) / (no-1));
   for (n=2; n<no; n++)           	/* sum the midpoints */
   {
      x    = interval * (n-1);      
      sum += f(x)*interval;
   }
   sum += 0.5 *(f(min) + f(max)) * interval;	/* add the endpoints */

   return (sum);
}      

/************************************************************************
/* Integration using Simpson's rule */ 
Number simpson (int no, Number min, Number max)
{  
   int n;				 
   Number interval, sum=0., x;
   interval = ((max -min) /(no-1));
   
   for (n=2; n<no; n+=2)                /* loop for odd points  */
   {
       x = interval * (n-1);
       sum += 4 * f(x);
   }
   for (n=3; n<no; n+=2)                /* loop for even points  */
   {
      x = interval * (n-1);
      sum += 2 * f(x);
   }   
   sum +=  f(min) + f(max);	 	/* add first and last value */          
   sum *= interval/3.;        		/* then multilpy by interval*/
   
   return (sum);
}  

/************************************************************************
/* Integration using Gauss' rule */ 
Number gaussint (int no, Number min, Number max)
{
   int n; 			 
   Number quadra = 0.;
   Number w[1000], x[1000];             /* for points and weights */
   
   gauss (no, 0, min, max, x, w);       /* returns Legendre */
   					/* points and weights */
   for (n=0; n< no; n++)
   {                               
      quadra += f(x[n])*w[n];            /* calculating the integral */
   }   
   return (quadra);                  
}
/***********************************************************************/
#include"butil.h"

Number romberg (int JMAX, Number (*func)(Number), Number min, Number max){  
   int n,j,k,ifat;
   Number sum, h, **M;
   M = matrix(JMAX,JMAX);
   
   h = max - min;
   n = 2;
   M[1][1] = h * (func(max) + func(min))/2 ;
   h = h / 2;
   M[2][1] = M[1][1]/2 + h*func(min + h) ;
   M[2][2] = (4*M[2][1] - M[1][1])/3;
   for(j=3;j<=JMAX;j++){
	   n = 2 * n;
	   h = h / 2;
	   sum = 0.;
	   for(k=1;k<=n;k+=2){sum += func(min + k*h);}
	   M[j][1] = M[j-1][1]/2 + h*sum;
	   ifat = 4;
	   for(k=2;k<=j;k++){
		   M[j][k] =(ifat*M[j][k-1] - M[j-1][k-1])/(ifat-1);
		   ifat = ifat * 4;
	   }
	   	   if(fabs(M[j][j] - M[j][j-1]) < 1.e-15*fabs(M[j][j])){
		   sum = M[j][j];
		   return sum;
		   }
   }
   sum = M[JMAX][JMAX];
   delete []M;   
   return sum;
}

/*    End Romberg integration  
/************************************************************************
*/