/* The function int1d() calculates a 1 dimensional definite integral
   according to the bisection method. 
   The routine works with a fixed number of iterations : n. 
   For a given iteration, the length is bisected n - 1 times 
   The number of grid points along each direction = pow(2,n - 1) + 1. */



#include <int.h>

#define FUNC(x) ((*func)(x))  

/* func is a pointer to the function to be integrated */

/* x1 and x2 are the lower and upper limits of the x integration */


double int1d(double (*func)(double x), double x1, double x2, int n )
{
     double tnm, corners ,boundary = 0.0, delx;
     double s = 0.0;
     long i,k,it;

     corners = FUNC(x1) + FUNC(x2);
     if (n == 1)
     {
	s = .5 * ( x2 - x1) * corners;
	return s;
     }
     else
     {
	for (it = 1, k = 1; k < n; k++)
	   it = it << 1;
	tnm = it;
	delx = (x2 - x1) / tnm;
	
	for (i=1; i < it; i++)
	  boundary += FUNC(x1 + i * delx);
	
	s = delx * (.5 * corners + boundary);
	return s;
     }

}