/* cube_int.c  integrate a function over a cube                 */
/*             this is numerical quadrature                     */
/*             basically compute a 4D volume under 3D function  */
/*             compare methods for accuracy                     */

/* first trapezoid rule, varying step size, int f(x,y,z) dx dx dz */
/* then using Gauss-Legendre coordinates and weights            */

/* remember np points or nodes means np*np*np function evaluations */
/* or (np-1)^3  cube elements                                   */

#include <stdio.h>
#include <math.h>
#undef  abs
#define abs(x)  ((x)<0.0?(-(x)):(x))

static void gaulegf(double x1, double x2, double x[], double w[], int n);

static double f1(double x, double y, double z)
{
  return x+2.0*y+3.0*z+4.0; /* linear, first order */
}

static double f2(double x, double y, double z)
{
  return x*x+2.0*y*y+3.0*z*z+4.0*x*y+5.0*x*z+6.0*y*z+
         7.0*x+8.0*y+9.0*z+10.0; /* quadratic, second order */
}

static double f3(double x, double y, double z)
{
  return x*x*x + 2.0*y*y*y + 3.0*z*z*z + 4.0*x*x*y +
         5.0*x*x*z + 6.0*x*y*y + 7.0*x*y*z + 8.0*x*z*z +
         9.0*y*y*z + 10.0*y*z*z +
         11.0*x*x+12.0*y*y+13.0*z*z+14.0*x*y+15.0*x*z+16.0*y*z+
         17.0*x+18.0*y+19.0*z+20.0; /* cubic, third order */
}

static double fe(double x, double y, double z)
{
  return exp(x)*exp(y)*exp(z);
}

int main(int argc, char *argv[])
{
  /* simple trapezoidal integration from x=1 to x=3*/
  double volume;
  double x, y, z, hx, hy, hz;
  double xmin = 1.0;
  double xmax = 3.0;
  double ymin = 0.5;
  double ymax = 2.0;
  double zmin = 0.25;
  double zmax = 2.25;
  double f1exact = 73.5;
  double f2exact = 543.875;
  double f3exact = 2887.96875;
  double feexact = 817.8595522;
  double xx[200], wx[200]; /* for Gauss-Legendre */
  double yy[200], wy[200]; /* for Gauss-Legendre */
  double zz[200], wz[200]; /* for Gauss-Legendre */
  int i, j, k, np;

  printf("cube_int.c integrate to find 4D volume \n");
  printf("Numerical quadrature, x=%g to x=%g, y=%g to y=%g, z=%g to z=%g \n",
          xmin, xmax, ymin, ymax, zmin, zmax);
  printf("\n");
  printf("Using trapezoidal method \n");
  np = 4;
  hx = (xmax-xmin)/(double)(np-1);
  hy = (ymax-ymin)/(double)(np-1);
  hz = (zmax-zmin)/(double)(np-1);
  printf("f1(x,y,z)=x+2*y+3*z+4  np=%d,  hx=%g, hy=%g, hz=%g \n",
          np, hx, hy, hz);
  volume = 0.0;
  for(i=0; i<np; i++)
  {
    x = xmin+(double)i*hx;
    for(j=0; j<np; j++)
    {
      y = ymin+(double)j*hy;
      volume = volume + (f1(x,y,zmin) + f1(x,y,zmax))/2.0;
    }
  }
  for(i=0; i<np; i++)
  {
    x = xmin+(double)i*hx;
    for(k=0; k<np; k++)
    {
      z = zmin+(double)k*hz;
      volume = volume + (f1(x,ymin,z) + f1(x,ymax,z))/2.0;
    }
  }
  for(j=0; j<np; j++)
  {
    y = ymin+(double)j*hy;
    for(k=0; k<np; k++)
    {
      z = zmin+(double)k*hz;
      volume = volume + (f1(xmin,y,z) + f1(xmax,y,z))/2.0;
    }
  }
  for(i=1; i<np-1; i++)
  {
    x = xmin+(double)i*hx;
    for(j=1; j<np-1; j++)
    {
      y = ymin+(double)j*hy;
      for(k=1; k<np-1; k++)
      { 
        z = zmin+(double)k*hz;
        volume = volume + f1(x,y,z);
      }
    }
  }
  volume = (xmax-xmin)*(ymax-ymin)*(zmax-zmin)*volume/
           (double)((np-1)*(np-1)*(np-1));
  printf("%d function evaluations \n", np*np*np);
  printf("computed volume=%g, exact volume=%g, error=%g \n",
         volume, f1exact, volume-f1exact);

  np = 16;
  hx = (xmax-xmin)/(double)(np-1);
  hy = (ymax-ymin)/(double)(np-1);
  hz = (zmax-zmin)/(double)(np-1);
  printf("f1(x,y,z)=x+2*y+3*z+4  np=%d,  hx=%g, hy=%g, hz=%g \n",
          np, hx, hy, hz);
  volume = 0.0;
  for(i=0; i<np; i++)
  {
    x = xmin+(double)i*hx;
    for(j=0; j<np; j++)
    {
      y = ymin+(double)j*hy;
      volume = volume + (f1(x,y,zmin) + f1(x,y,zmax))/2.0;
    }
  }
  for(i=0; i<np; i++)
  {
    x = xmin+(double)i*hx;
    for(k=0; k<np; k++)
    {
      z = zmin+(double)k*hz;
      volume = volume + (f1(x,ymin,z) + f1(x,ymax,z))/2.0;
    }
  }
  for(j=0; j<np; j++)
  {
    y = ymin+(double)j*hy;
    for(k=0; k<np; k++)
    {
      z = zmin+(double)k*hz;
      volume = volume + (f1(xmin,y,z) + f1(xmax,y,z))/2.0;
    }
  }
  for(i=1; i<np-1; i++)
  {
    x = xmin+(double)i*hx;
    for(j=1; j<np-1; j++)
    {
      y = ymin+(double)j*hy;
      for(k=1; k<np-1; k++)
      {
        z = zmin+(double)k*hz;
        volume = volume + f1(x,y,z);
      }
    }
  }
  volume = (xmax-xmin)*(ymax-ymin)*(zmax-zmin)*volume/
           (double)((np-1)*(np-1)*(np-1));
  printf("%d function evaluations \n", np*np*np);
  printf("computed volume=%g, exact volume=%g, error=%g \n",
         volume, f1exact, volume-f1exact);

  np = 64;
  hx = (xmax-xmin)/(double)(np-1);
  hy = (ymax-ymin)/(double)(np-1);
  hz = (zmax-zmin)/(double)(np-1);
  printf("f1(x,y,z)=x+2*y+3*z+4  np=%d,  hx=%g, hy=%g, hz=%g \n",
          np, hx, hy, hz);
  volume = 0.0;
  for(i=0; i<np; i++)
  {
    x = xmin+(double)i*hx;
    for(j=0; j<np; j++)
    {
      y = ymin+(double)j*hy;
      volume = volume + (f1(x,y,zmin) + f1(x,y,zmax))/2.0;
    }
  }
  for(i=0; i<np; i++)
  {
    x = xmin+(double)i*hx;
    for(k=0; k<np; k++)
    {
      z = zmin+(double)k*hz;
      volume = volume + (f1(x,ymin,z) + f1(x,ymax,z))/2.0;
    }
  }
  for(j=0; j<np; j++)
  {
    y = ymin+(double)j*hy;
    for(k=0; k<np; k++)
    {
      z = zmin+(double)k*hz;
      volume = volume + (f1(xmin,y,z) + f1(xmax,y,z))/2.0;
    }
  }
  for(i=1; i<np-1; i++)
  {
    x = xmin+(double)i*hx;
    for(j=1; j<np-1; j++)
    {
      y = ymin+(double)j*hy;
      for(k=1; k<np-1; k++)
      {
        z = zmin+(double)k*hz;
        volume = volume + f1(x,y,z);
      }
    }
  }
  volume = (xmax-xmin)*(ymax-ymin)*(zmax-zmin)*volume/
           (double)((np-1)*(np-1)*(np-1));
  printf("%d function evaluations \n", np*np*np);
  printf("computed volume=%g, exact volume=%g, error=%g \n",
         volume, f1exact, volume-f1exact);
  printf("\n");


  np = 4;
  hx = (xmax-xmin)/(double)(np-1);
  hy = (ymax-ymin)/(double)(np-1);
  hz = (zmax-zmin)/(double)(np-1);
  printf("f2(x,y,z)=x*x+2*y*y+3*z*z+4*x*y+5*x*z+6*y*z+ \n");
  printf("          7*x+8*y+9*z+10  np=%d,  hx=%g, hy=%g, hz=%g \n",
          np, hx, hy, hz);
  volume = 0.0;
  for(i=0; i<np; i++)
  {
    x = xmin+(double)i*hx;
    for(j=0; j<np; j++)
    {
      y = ymin+(double)j*hy;
      volume = volume + (f2(x,y,zmin) + f2(x,y,zmax))/2.0;
    }
  }
  for(i=0; i<np; i++)
  {
    x = xmin+(double)i*hx;
    for(k=0; k<np; k++)
    {
      z = zmin+(double)k*hz;
      volume = volume + (f2(x,ymin,z) + f2(x,ymax,z))/2.0;
    }
  }
  for(j=0; j<np; j++)
  {
    y = ymin+(double)j*hy;
    for(k=0; k<np; k++)
    {
      z = zmin+(double)k*hz;
      volume = volume + (f2(xmin,y,z) + f2(xmax,y,z))/2.0;
    }
  }
  for(i=1; i<np-1; i++)
  {
    x = xmin+(double)i*hx;
    for(j=1; j<np-1; j++)
    {
      y = ymin+(double)j*hy;
      for(k=1; k<np-1; k++)
      {
        z = zmin+(double)k*hz;
        volume = volume + f2(x,y,z);
      }
    }
  }
  printf("%d function evaluations \n", np*np*np);
  volume = (xmax-xmin)*(ymax-ymin)*(zmax-zmin)*volume/
           (double)((np-1)*(np-1)*(np-1));
  printf("computed volume=%g, exact volume=%g, error=%g \n",
         volume, f2exact, volume-f2exact);

  np = 16;
  hx = (xmax-xmin)/(double)(np-1);
  hy = (ymax-ymin)/(double)(np-1);
  hz = (zmax-zmin)/(double)(np-1);
  printf("f2(x,y,z)=x*x+2*y*y+3*z*z+4*x*y+5*x*z+6*y*z+ \n");
  printf("          7*x+8*y+9*z+10  np=%d,  hx=%g, hy=%g, hz=%g \n",
          np, hx, hy, hz);
  volume = 0.0;
  for(i=0; i<np; i++)
  {
    x = xmin+(double)i*hx;
    for(j=0; j<np; j++)
    {
      y = ymin+(double)j*hy;
      volume = volume + (f2(x,y,zmin) + f2(x,y,zmax))/2.0;
    }
  }
  for(i=0; i<np; i++)
  {
    x = xmin+(double)i*hx;
    for(k=0; k<np; k++)
    {
      z = zmin+(double)k*hz;
      volume = volume + (f2(x,ymin,z) + f2(x,ymax,z))/2.0;
    }
  }
  for(j=0; j<np; j++)
  {
    y = ymin+(double)j*hy;
    for(k=0; k<np; k++)
    {
      z = zmin+(double)k*hz;
      volume = volume + (f2(xmin,y,z) + f2(xmax,y,z))/2.0;
    }
  }
  for(i=1; i<np-1; i++)
  {
    x = xmin+(double)i*hx;
    for(j=1; j<np-1; j++)
    {
      y = ymin+(double)j*hy;
      for(k=1; k<np-1; k++)
      {
        z = zmin+(double)k*hz;
        volume = volume + f2(x,y,z);
      }
    }
  }
  volume = (xmax-xmin)*(ymax-ymin)*(zmax-zmin)*volume/
           (double)((np-1)*(np-1)*(np-1));
  printf("%d function evaluations \n", np*np*np);
  printf("computed volume=%g, exact volume=%g, error=%g \n",
         volume, f2exact, volume-f2exact);

  np = 64;
  hx = (xmax-xmin)/(double)(np-1);
  hy = (ymax-ymin)/(double)(np-1);
  hz = (zmax-zmin)/(double)(np-1);
  printf("f2(x,y,z)=x*x+2*y*y+3*z*z+4*x*y+5*x*z+6*y*z+ \n");
  printf("          7*x+8*y+9*z+10  np=%d,  hx=%g, hy=%g, hz=%g \n",
          np, hx, hy, hz);
  volume = 0.0;
  for(i=0; i<np; i++)
  {
    x = xmin+(double)i*hx;
    for(j=0; j<np; j++)
    {
      y = ymin+(double)j*hy;
      volume = volume + (f2(x,y,zmin) + f2(x,y,zmax))/2.0;
    }
  }
  for(i=0; i<np; i++)
  {
    x = xmin+(double)i*hx;
    for(k=0; k<np; k++)
    {
      z = zmin+(double)k*hz;
      volume = volume + (f2(x,ymin,z) + f2(x,ymax,z))/2.0;
    }
  }
  for(j=0; j<np; j++)
  {
    y = ymin+(double)j*hy;
    for(k=0; k<np; k++)
    {
      z = zmin+(double)k*hz;
      volume = volume + (f2(xmin,y,z) + f2(xmax,y,z))/2.0;
    }
  }
  for(i=1; i<np-1; i++)
  {
    x = xmin+(double)i*hx;
    for(j=1; j<np-1; j++)
    {
      y = ymin+(double)j*hy;
      for(k=1; k<np-1; k++)
      {
        z = zmin+(double)k*hz;
        volume = volume + f2(x,y,z);
      }
    }
  }
  volume = (xmax-xmin)*(ymax-ymin)*(zmax-zmin)*volume/
           (double)((np-1)*(np-1)*(np-1));
  printf("%d function evaluations \n", np*np*np);
  printf("computed volume=%g, exact volume=%g, error=%g \n",
         volume, f2exact, volume-f2exact);

  np = 256;
  hx = (xmax-xmin)/(double)(np-1);
  hy = (ymax-ymin)/(double)(np-1);
  hz = (zmax-zmin)/(double)(np-1);
  printf("f2(x,y,z)=x*x+2*y*y+3*z*z+4*x*y+5*x*z+6*y*z+ \n");
  printf("          7*x+8*y+9*z+10  np=%d,  hx=%g, hy=%g, hz=%g \n",
          np, hx, hy, hz);
  volume = 0.0;
  for(i=0; i<np; i++)
  {
    x = xmin+(double)i*hx;
    for(j=0; j<np; j++)
    {
      y = ymin+(double)j*hy;
      volume = volume + (f2(x,y,zmin) + f2(x,y,zmax))/2.0;
    }
  }
  for(i=0; i<np; i++)
  {
    x = xmin+(double)i*hx;
    for(k=0; k<np; k++)
    {
      z = zmin+(double)k*hz;
      volume = volume + (f2(x,ymin,z) + f2(x,ymax,z))/2.0;
    }
  }
  for(j=0; j<np; j++)
  {
    y = ymin+(double)j*hy;
    for(k=0; k<np; k++)
    {
      z = zmin+(double)k*hz;
      volume = volume + (f2(xmin,y,z) + f2(xmax,y,z))/2.0;
    }
  }
  for(i=1; i<np-1; i++)
  {
    x = xmin+(double)i*hx;
    for(j=1; j<np-1; j++)
    {
      y = ymin+(double)j*hy;
      for(k=1; k<np-1; k++)
      {
        z = zmin+(double)k*hz;
        volume = volume + f2(x,y,z);
      }
    }
  }
  volume = (xmax-xmin)*(ymax-ymin)*(zmax-zmin)*volume/
           (double)((np-1)*(np-1)*(np-1));
  printf("%d function evaluations \n", np*np*np);
  printf("computed volume=%g, exact volume=%g, error=%g \n",
         volume, f2exact, volume-f2exact);
  printf("\n");


  np = 4;
  hx = (xmax-xmin)/(double)(np-1);
  hy = (ymax-ymin)/(double)(np-1);
  hz = (zmax-zmin)/(double)(np-1);
  printf("f3(x,y,z)=x*x*x+2*y*y*y+3*z*z*z+4*x*x*y+5*x*x*z+6*x*y*y + \n");
  printf("          7*x*y*z+8*x*z*z+9*y*y*z+10*y*z*z+11*x*x+12*y*y + \n");
  printf("          13*z*z+14*x*y+15*x*z+16*y*z+17*x+18*y+19*z+20 \n");
  printf("          np=%d,  hx=%g, hy=%g, hz=%g \n", np, hx, hy, hz);
  volume = 0.0;
  for(i=0; i<np; i++)
  {
    x = xmin+(double)i*hx;
    for(j=0; j<np; j++)
    {
      y = ymin+(double)j*hy;
      volume = volume + (f3(x,y,zmin) + f3(x,y,zmax))/2.0;
    }
  }
  for(i=0; i<np; i++)
  {
    x = xmin+(double)i*hx;
    for(k=0; k<np; k++)
    {
      z = zmin+(double)k*hz;
      volume = volume + (f3(x,ymin,z) + f3(x,ymax,z))/2.0;
    }
  }
  for(j=0; j<np; j++)
  {
    y = ymin+(double)j*hy;
    for(k=0; k<np; k++)
    {
      z = zmin+(double)k*hz;
      volume = volume + (f3(xmin,y,z) + f3(xmax,y,z))/2.0;
    }
  }
  for(i=1; i<np-1; i++)
  {
    x = xmin+(double)i*hx;
    for(j=1; j<np-1; j++)
    {
      y = ymin+(double)j*hy;
      for(k=1; k<np-1; k++)
      {
        z = zmin+(double)k*hz;
        volume = volume + f3(x,y,z);
      }
    }
  }
  volume = (xmax-xmin)*(ymax-ymin)*(zmax-zmin)*volume/
           (double)((np-1)*(np-1)*(np-1));
  printf("%d function evaluations \n", np*np*np);
  printf("computed volume=%g, exact volume=%g, error=%g \n",
         volume, f3exact, volume-f3exact);

  np = 16;
  hx = (xmax-xmin)/(double)(np-1);
  hy = (ymax-ymin)/(double)(np-1);
  hz = (zmax-zmin)/(double)(np-1);
  printf("f3(x,y,z)=x*x*x+2*y*y*y+3*z*z*z+4*x*x*y+5*x*x*z+6*x*y*y + \n");
  printf("          7*x*y*z+8*x*z*z+9*y*y*z+10*y*z*z+11*x*x+12*y*y + \n");
  printf("          13*z*z+14*x*y+15*x*z+16*y*z+17*x+18*y+19*z+20 \n");
  printf("          np=%d,  hx=%g, hy=%g, hz=%g \n", np, hx, hy, hz);
  volume = 0.0;
  for(i=0; i<np; i++)
  {
    x = xmin+(double)i*hx;
    for(j=0; j<np; j++)
    {
      y = ymin+(double)j*hy;
      volume = volume + (f3(x,y,zmin) + f3(x,y,zmax))/2.0;
    }
  }
  for(i=0; i<np; i++)
  {
    x = xmin+(double)i*hx;
    for(k=0; k<np; k++)
    {
      z = zmin+(double)k*hz;
      volume = volume + (f3(x,ymin,z) + f3(x,ymax,z))/2.0;
    }
  }
  for(j=0; j<np; j++)
  {
    y = ymin+(double)j*hy;
    for(k=0; k<np; k++)
    {
      z = zmin+(double)k*hz;
      volume = volume + (f3(xmin,y,z) + f3(xmax,y,z))/2.0;
    }
  }
  for(i=1; i<np-1; i++)
  {
    x = xmin+(double)i*hx;
    for(j=1; j<np-1; j++)
    {
      y = ymin+(double)j*hy;
      for(k=1; k<np-1; k++)
      {
        z = zmin+(double)k*hz;
        volume = volume + f3(x,y,z);
      }
    }
  }
  volume = (xmax-xmin)*(ymax-ymin)*(zmax-zmin)*volume/
           (double)((np-1)*(np-1)*(np-1));
  printf("%d function evaluations \n", np*np*np);
  printf("computed volume=%g, exact volume=%g, error=%g \n",
         volume, f3exact, volume-f3exact);

  np = 64;
  hx = (xmax-xmin)/(double)(np-1);
  hy = (ymax-ymin)/(double)(np-1);
  hz = (zmax-zmin)/(double)(np-1);
  printf("f3(x,y,z)=x*x*x+2*y*y*y+3*z*z*z+4*x*x*y+5*x*x*z+6*x*y*y + \n");
  printf("          7*x*y*z+8*x*z*z+9*y*y*z+10*y*z*z+11*x*x+12*y*y + \n");
  printf("          13*z*z+14*x*y+15*x*z+16*y*z+17*x+18*y+19*z+20 \n");
  printf("          np=%d,  hx=%g, hy=%g, hz=%g \n", np, hx, hy, hz);
  volume = 0.0;
  for(i=0; i<np; i++)
  {
    x = xmin+(double)i*hx;
    for(j=0; j<np; j++)
    {
      y = ymin+(double)j*hy;
      volume = volume + (f3(x,y,zmin) + f3(x,y,zmax))/2.0;
    }
  }
  for(i=0; i<np; i++)
  {
    x = xmin+(double)i*hx;
    for(k=0; k<np; k++)
    {
      z = zmin+(double)k*hz;
      volume = volume + (f3(x,ymin,z) + f3(x,ymax,z))/2.0;
    }
  }
  for(j=0; j<np; j++)
  {
    y = ymin+(double)j*hy;
    for(k=0; k<np; k++)
    {
      z = zmin+(double)k*hz;
      volume = volume + (f3(xmin,y,z) + f3(xmax,y,z))/2.0;
    }
  }
  for(i=1; i<np-1; i++)
  {
    x = xmin+(double)i*hx;
    for(j=1; j<np-1; j++)
    {
      y = ymin+(double)j*hy;
      for(k=1; k<np-1; k++)
      {
        z = zmin+(double)k*hz;
        volume = volume + f3(x,y,z);
      }
    }
  }
  volume = (xmax-xmin)*(ymax-ymin)*(zmax-zmin)*volume/
           (double)((np-1)*(np-1)*(np-1));
  printf("%d function evaluations \n", np*np*np);
  printf("computed volume=%g, exact volume=%g, error=%g \n",
         volume, f3exact, volume-f3exact);

  np = 256;
  hx = (xmax-xmin)/(double)(np-1);
  hy = (ymax-ymin)/(double)(np-1);
  hz = (zmax-zmin)/(double)(np-1);
  printf("f3(x,y,z)=x*x*x+2*y*y*y+3*z*z*z+4*x*x*y+5*x*x*z+6*x*y*y + \n");
  printf("          7*x*y*z+8*x*z*z+9*y*y*z+10*y*z*z+11*x*x+12*y*y + \n");
  printf("          13*z*z+14*x*y+15*x*z+16*y*z+17*x+18*y+19*z+20 \n");
  printf("          np=%d,  hx=%g, hy=%g, hz=%g \n", np, hx, hy, hz);
  volume = 0.0;
  for(i=0; i<np; i++)
  {
    x = xmin+(double)i*hx;
    for(j=0; j<np; j++)
    {
      y = ymin+(double)j*hy;
      volume = volume + (f3(x,y,zmin) + f3(x,y,zmax))/2.0;
    }
  }
  for(i=0; i<np; i++)
  {
    x = xmin+(double)i*hx;
    for(k=0; k<np; k++)
    {
      z = zmin+(double)k*hz;
      volume = volume + (f3(x,ymin,z) + f3(x,ymax,z))/2.0;
    }
  }
  for(j=0; j<np; j++)
  {
    y = ymin+(double)j*hy;
    for(k=0; k<np; k++)
    {
      z = zmin+(double)k*hz;
      volume = volume + (f3(xmin,y,z) + f3(xmax,y,z))/2.0;
    }
  }
  for(i=1; i<np-1; i++)
  {
    x = xmin+(double)i*hx;
    for(j=1; j<np-1; j++)
    {
      y = ymin+(double)j*hy;
      for(k=1; k<np-1; k++)
      {
        z = zmin+(double)k*hz;
        volume = volume + f3(x,y,z);
      }
    }
  }
  volume = (xmax-xmin)*(ymax-ymin)*(zmax-zmin)*volume/
           (double)((np-1)*(np-1)*(np-1));
  printf("%d function evaluations \n", np*np*np);
  printf("computed volume=%g, exact volume=%g, error=%g \n",
         volume, f3exact, volume-f3exact);
  printf("\n");

  np = 4;
  hx = (xmax-xmin)/(double)(np-1);
  hy = (ymax-ymin)/(double)(np-1);
  hz = (zmax-zmin)/(double)(np-1);
  printf("fe(x,y,z)=exp(x)*exp(y)*exp(z) \n");
  printf("          np=%d,  hx=%g, hy=%g, hz=%g \n", np, hx, hy, hz);
  volume = 0.0;
  for(i=0; i<np; i++)
  {
    x = xmin+(double)i*hx;
    for(j=0; j<np; j++)
    {
      y = ymin+(double)j*hy;
      volume = volume + (fe(x,y,zmin) + fe(x,y,zmax))/2.0;
    }
  }
  for(i=0; i<np; i++)
  {
    x = xmin+(double)i*hx;
    for(k=0; k<np; k++)
    {
      z = zmin+(double)k*hz;
      volume = volume + (fe(x,ymin,z) + fe(x,ymax,z))/2.0;
    }
  }
  for(j=0; j<np; j++)
  {
    y = ymin+(double)j*hy;
    for(k=0; k<np; k++)
    {
      z = zmin+(double)k*hz;
      volume = volume + (fe(xmin,y,z) + fe(xmax,y,z))/2.0;
    }
  }
  for(i=1; i<np-1; i++)
  {
    x = xmin+(double)i*hx;
    for(j=1; j<np-1; j++)
    {
      y = ymin+(double)j*hy;
      for(k=1; k<np-1; k++)
      {
        z = zmin+(double)k*hz;
        volume = volume + fe(x,y,z);
      }
    }
  }
  volume = (xmax-xmin)*(ymax-ymin)*(zmax-zmin)*volume/
           (double)((np-1)*(np-1)*(np-1));
  printf("%d function evaluations \n", np*np*np);
  printf("computed volume=%g, exact volume=%g, error=%g \n",
         volume, feexact, volume-feexact);

  np = 16;
  hx = (xmax-xmin)/(double)(np-1);
  hy = (ymax-ymin)/(double)(np-1);
  hz = (zmax-zmin)/(double)(np-1);
  printf("fe(x,y,z)=exp(x)*exp(y)*exp(z) \n");
  printf("          np=%d,  hx=%g, hy=%g, hz=%g \n", np, hx, hy, hz);
  volume = 0.0;
  for(i=0; i<np; i++)
  {
    x = xmin+(double)i*hx;
    for(j=0; j<np; j++)
    {
      y = ymin+(double)j*hy;
      volume = volume + (fe(x,y,zmin) + fe(x,y,zmax))/2.0;
    }
  }
  for(i=0; i<np; i++)
  {
    x = xmin+(double)i*hx;
    for(k=0; k<np; k++)
    {
      z = zmin+(double)k*hz;
      volume = volume + (fe(x,ymin,z) + fe(x,ymax,z))/2.0;
    }
  }
  for(j=0; j<np; j++)
  {
    y = ymin+(double)j*hy;
    for(k=0; k<np; k++)
    {
      z = zmin+(double)k*hz;
      volume = volume + (fe(xmin,y,z) + fe(xmax,y,z))/2.0;
    }
  }
  for(i=1; i<np-1; i++)
  {
    x = xmin+(double)i*hx;
    for(j=1; j<np-1; j++)
    {
      y = ymin+(double)j*hy;
      for(k=1; k<np-1; k++)
      {
        z = zmin+(double)k*hz;
        volume = volume + fe(x,y,z);
      }
    }
  }
  volume = (xmax-xmin)*(ymax-ymin)*(zmax-zmin)*volume/
           (double)((np-1)*(np-1)*(np-1));
  printf("%d function evaluations \n", np*np*np);
  printf("computed volume=%g, exact volume=%g, error=%g \n",
         volume, feexact, volume-feexact);

  np = 64;
  hx = (xmax-xmin)/(double)(np-1);
  hy = (ymax-ymin)/(double)(np-1);
  hz = (zmax-zmin)/(double)(np-1);
  printf("fe(x,y,z)=exp(x)*exp(y)*exp(z) \n");
  printf("          np=%d,  hx=%g, hy=%g, hz=%g \n", np, hx, hy, hz);
  volume = 0.0;
  for(i=0; i<np; i++)
  {
    x = xmin+(double)i*hx;
    for(j=0; j<np; j++)
    {
      y = ymin+(double)j*hy;
      volume = volume + (fe(x,y,zmin) + fe(x,y,zmax))/2.0;
    }
  }
  for(i=0; i<np; i++)
  {
    x = xmin+(double)i*hx;
    for(k=0; k<np; k++)
    {
      z = zmin+(double)k*hz;
      volume = volume + (fe(x,ymin,z) + fe(x,ymax,z))/2.0;
    }
  }
  for(j=0; j<np; j++)
  {
    y = ymin+(double)j*hy;
    for(k=0; k<np; k++)
    {
      z = zmin+(double)k*hz;
      volume = volume + (fe(xmin,y,z) + fe(xmax,y,z))/2.0;
    }
  }
  for(i=1; i<np-1; i++)
  {
    x = xmin+(double)i*hx;
    for(j=1; j<np-1; j++)
    {
      y = ymin+(double)j*hy;
      for(k=1; k<np-1; k++)
      {
        z = zmin+(double)k*hz;
        volume = volume + fe(x,y,z);
      }
    }
  }
  volume = (xmax-xmin)*(ymax-ymin)*(zmax-zmin)*volume/
           (double)((np-1)*(np-1)*(np-1));
  printf("%d function evaluations \n", np*np*np);
  printf("computed volume=%g, exact volume=%g, error=%g \n",
         volume, feexact, volume-feexact);

  np = 256;
  hx = (xmax-xmin)/(double)(np-1);
  hy = (ymax-ymin)/(double)(np-1);
  hz = (zmax-zmin)/(double)(np-1);
  printf("fe(x,y,z)=exp(x)*exp(y)*exp(z) \n");
  printf("          np=%d,  hx=%g, hy=%g, hz=%g \n", np, hx, hy, hz);
  volume = 0.0;
  for(i=0; i<np; i++)
  {
    x = xmin+(double)i*hx;
    for(j=0; j<np; j++)
    {
      y = ymin+(double)j*hy;
      volume = volume + (fe(x,y,zmin) + fe(x,y,zmax))/2.0;
    }
  }
  for(i=0; i<np; i++)
  {
    x = xmin+(double)i*hx;
    for(k=0; k<np; k++)
    {
      z = zmin+(double)k*hz;
      volume = volume + (fe(x,ymin,z) + fe(x,ymax,z))/2.0;
    }
  }
  for(j=0; j<np; j++)
  {
    y = ymin+(double)j*hy;
    for(k=0; k<np; k++)
    {
      z = zmin+(double)k*hz;
      volume = volume + (fe(xmin,y,z) + fe(xmax,y,z))/2.0;
    }
  }
  for(i=1; i<np-1; i++)
  {
    x = xmin+(double)i*hx;
    for(j=1; j<np-1; j++)
    {
      y = ymin+(double)j*hy;
      for(k=1; k<np-1; k++)
      {
        z = zmin+(double)k*hz;
        volume = volume + fe(x,y,z);
      }
    }
  }
  volume = (xmax-xmin)*(ymax-ymin)*(zmax-zmin)*volume/
           (double)((np-1)*(np-1)*(np-1));
  printf("%d function evaluations \n", np*np*np);
  printf("computed volume=%g, exact volume=%g, error=%g \n",
         volume, feexact, volume-feexact);
  printf("\n");


  printf("Using Gauss Legendre method \n");

  np = 4;
  printf("f1(x,y,z)=x+2y+3z+4 np=%d \n", np);
  gaulegf(xmin, xmax, xx, wx, np);
  gaulegf(ymin, ymax, yy, wy, np);
  gaulegf(zmin, zmax, zz, wz, np);
  volume = 0.0;
  for(i=1; i<=np; i++)
  {
    for(j=1; j<=np; j++)
    {
      for(k=1; k<=np; k++)
      {
        volume = volume + wx[i]*wy[j]*wz[k]*f1(xx[i], yy[j], zz[k]);
      }
    }
  }
  printf("%d function evaluations \n", np*np*np);
  printf("computed volume=%g, exact volume=%g, error=%g \n",
         volume, f1exact, volume-f1exact);
  printf("\n");

  np = 4;
  printf("f2(x,y,z)=x*x+2*y*y+3*z*z+4*x*y+5*x*z+6*y*z+ \n");
  printf("          7*x+8*y+9*z+10   np=%d \n", np);
  gaulegf(xmin, xmax, xx, wx, np);
  gaulegf(ymin, ymax, yy, wy, np);
  gaulegf(zmin, zmax, zz, wz, np);
  volume = 0.0;
  for(i=1; i<=np; i++)
  {
    for(j=1; j<=np; j++)
    {
      for(k=1; k<=np; k++)
      {
        volume = volume + wx[i]*wy[j]*wz[k]*f2(xx[i], yy[j], zz[k]);
      }
    }
  }
  printf("%d function evaluations \n", np*np*np);
  printf("computed volume=%g, exact volume=%g, error=%g \n",
         volume, f2exact, volume-f2exact);
  printf("\n");

  np = 4;
  printf("f3(x,y,z)=x*x*x+2*y*y*y+3*z*z*z+4*x*x*y+5*x*x*z+6*x*y*y + \n");
  printf("          7*x*y*z+8*x*z*z+9*y*y*z+10*y*z*z+11*x*x+12*y*y + \n");
  printf("          13*z*z+14*x*y+15*x*z+16*y*z+17*x+18*y+19*z+20  np=%d \n", np);
  gaulegf(xmin, xmax, xx, wx, np);
  gaulegf(ymin, ymax, yy, wy, np);
  gaulegf(zmin, zmax, zz, wz, np);
  volume = 0.0;
  for(i=1; i<=np; i++)
  {
    for(j=1; j<=np; j++)
    {
      for(k=1; k<=np; k++)
      {
        volume = volume + wx[i]*wy[j]*wz[k]*f3(xx[i], yy[j], zz[k]);
      }
    }
  }
  printf("%d function evaluations \n", np*np*np);
  printf("computed volume=%g, exact volume=%g, error=%g \n",
         volume, f3exact, volume-f3exact);
  printf("\n");

  np = 4;
  printf("fe(x,y,z)=exp(x)*exp(y)*exp(z)  np=%d \n", np);
  gaulegf(xmin, xmax, xx, wx, np);
  gaulegf(ymin, ymax, yy, wy, np);
  gaulegf(zmin, zmax, zz, wz, np);
  volume = 0.0;
  for(i=1; i<=np; i++)
  {
    for(j=1; j<=np; j++)
    {
      for(k=1; k<=np; k++)
      {
        volume = volume + wx[i]*wy[j]*wz[k]*fe(xx[i], yy[j], zz[k]);
      }
    }
  }
  printf("%d function evaluations \n", np*np*np);
  printf("computed volume=%g, exact volume=%g, error=%g \n",
         volume, feexact, volume-feexact);

  np = 8;
  printf("fe(x,y,z)=exp(x)*exp(y)*exp(z)  np=%d \n", np);
  gaulegf(xmin, xmax, xx, wx, np);
  gaulegf(ymin, ymax, yy, wy, np);
  gaulegf(zmin, zmax, zz, wz, np);
  volume = 0.0;
  for(i=1; i<=np; i++)
  {
    for(j=1; j<=np; j++)
    {
      for(k=1; k<=np; k++)
      {
        volume = volume + wx[i]*wy[j]*wz[k]*fe(xx[i], yy[j], zz[k]);
      }
    }
  }
  printf("%d function evaluations \n", np*np*np);
  printf("computed volume=%g, exact volume=%g, error=%g \n",
         volume, feexact, volume-feexact);
  printf("\n");

  return 0;
}

static void gaulegf(double x1, double x2, double x[], double w[], int n)
{
  int i, j, m;
  double eps = 3.0E-14;
  double p1, p2, p3, pp, xl, xm, z, z1;
  
  m = (n+1)/2;
  xm = 0.5*(x2+x1);
  xl = 0.5*(x2-x1);
  for(i=1; i<=m; i++)
  {
    z = cos(3.141592654*((double)i-0.25)/((double)n+0.5));
    while(1)
    {
      p1 = 1.0;
      p2 = 0.0;
      for(j=1; j<=n; j++)
      {
        p3 = p2;
        p2 = p1;
        p1 = ((2.0*(double)j-1.0)*z*p2-((double)j-1.0)*p3)/
              (double)j;
      }
      pp = (double)n*(z*p1-p2)/(z*z-1.0);
      z1 = z;
      z = z1 - p1/pp;
      if(abs(z-z1) <= eps) break;
    }
    x[i] = xm - xl*z;
    x[n+1-i] = xm + xl*z;
    w[i] = 2.0*xl/((1.0-z*z)*pp*pp);
    w[n+1-i] = w[i];
  }
} /* end gaulegf */
/* end cube_int.c */