/* The function int3d() calculates a 3 dimensional definite integral  
   by overlaying the integration region with a 3 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,z) ((*func)(x,y,z))  

/* 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 */
/* z1 and z2 are the lower and upper limits of the z integration */


double int3d(double (*func)(double x, double y, double z),
double x1, double x2, double y1, double y2, double z1, double z2, int n )
{
     double tnm, corners = 0.0, sides = 0., faces = 0., interior = 0., 
     delx, dely, delz , s = 0.0;

     long i,j,k,l,it;
     
     corners = ( FUNC(x1,y1,z1) + FUNC(x1,y1,z2) + FUNC(x1,y2,z1) + 
     FUNC(x1,y2,z2) + FUNC(x2,y1,z1) + FUNC(x2,y1,z2) + 
     FUNC(x2,y2,z1) + FUNC(x2,y2,z2) );
     
     if (n == 1)
     {
	s = .125 * ( x2 - x1) * ( y2 - y1) * (z2 - z1) * corners;
	return s;
     }
     else
     {
	for (it = 1, l = 1; l < n; l++)
	   it = it << 1;
	tnm = it;
	delx = (x2 - x1) / tnm;
	dely = (y2 - y1) / tnm;
	delz = (z2 - z1) / tnm;

	for (i=1; i < it; i++)
	   sides += FUNC(x1 + i * delx,y1,z1) + FUNC(x1 + i * delx,y1,z2) +  
           FUNC(x1 + i * delx,y2,z1) + FUNC(x1 + i * delx,y2,z2)  + 
           FUNC(x1,y1 + i * dely,z1) + FUNC(x1,y1 + i * dely,z2) + 
           FUNC(x2,y1 + i * dely,z1) + FUNC(x2,y1 + i * dely,z2) + 
           FUNC(x1,y1,z1 + i * delz) + FUNC(x1,y2,z1 + i * delz) + 
           FUNC(x2,y1,z1 + i * delz) + FUNC(x2,y2,z1 + i * delz) ;

	for (i = 1; i < it; i++)
	{
	   for (j = 1; j < it; j++)
	      faces += FUNC(x1 + i * delx, y1 + j * dely,z1 ) + 
              FUNC(x1 + i * delx, y1, z1 + j * delz) + 
              FUNC(x1 + i * delx, y2, z1 + j * delz) + 
              FUNC(x1 + i * delx, y1 + j * dely,z2) + 
              FUNC(x1, y1  + i * dely, z1 + j * delz) + 
              FUNC(x2, y1  + i * dely, z1 + j * delz) ;
	}

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

	s =  delx * dely * delz * (.125 * corners + .25 * sides + 
        .5 * faces + interior);
	
        return s;
     }

}