/* The function int2d() calculates a 2 dimensional definite integral  */
/* by overlaying the integration region with a 2 dim  grid. */
/* The routine works with a fixed number of iterations : n. */
/* For a given iteration, the length along each dimension is */
/* bisected n - 1 times */
/* number of grid points along each direction : pow(2,n - 1) + 1 */



#include <int.h>

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

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

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


double int2d(double (*func)(double x, double y), double x1, double x2, 
double y1, double y2, int n )
{
     double tnm, corners ,boundary = 0.0, interior = 0.0, delx, dely;
     double s = 0.0;
     long i,j,k,it;

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

	for (j=1; j < it; j++)
	   boundary += FUNC((x1 + j * delx),y1) + FUNC((x1 + j * delx),y2) + 
           FUNC(x1,(y1 + j * dely)) + FUNC(x2,(y1 + j * dely)) ;

	for (i = 1; i < it; i++)
	{
	   for (j = 1; j < it; j++)
	   {
	      interior += FUNC((x1 + i * delx), (y1 + j * dely));
	   }
	}

	s = delx * dely * (.25 * corners + .5 * boundary + interior);
	return s;
     }

}