// test_gaulegf3D.java  compile    javac -cp . test_gaulegf3D.java
//                      execute    java  -cp . test_gaulegf3D
//  intergral xmin..xmax, ymin..ymax zmin..zmax f(x,y,z) dx dy dz
// f(x,y,z) = x*z+y*y+2*x
// fx(x,y,z) = 1/2 x*x*z + y*y*x + x*x                     int x
// fxy(x,y,z) = 1/2 x*x*z*y + 1/3 y*y*y*x + x*x*y          int xy
// fxyz(x,y,z) =  1/4 x*x*z*z*y + 1/3 y*y*y*x*z + x*x*y*z  int xyz

public class test_gaulegf3D
{
  double xmin = 0.1;
  double xmax = 1.0;
  double ymin = 0.0;
  double ymax = 1.1;
  double zmin = 0.2;
  double zmax = 0.9;
  int nx = 6;
  int ny = 6;
  int nz = 6;
  double xx[] = new double[10];
  double wx[] = new double[10];
  double yy[] = new double[10];
  double wy[] = new double[10];
  double zz[] = new double[10];
  double wz[] = new double[10];
  double val = 0.0;    
  double v;
  double err = 0.0;
  double vchk, intchk, intchk1, intchk2;
  double dx = (xmax-xmin)/(double)(nx-1);
  double dy = (ymax-ymin)/(double)(ny-1);
  double dz = (zmax-zmin)/(double)(nz-1);
  double dv = dx*dy*dz;
  double x1, x2, y1, y2, z1, z2;

  public test_gaulegf3D()
  {
    System.out.println("test_gaulegf3D.java running");
    System.out.println("f(x,y,z)=x*z + y*y + 2.0*x ");
    System.out.println("f(x,y,z) = 2*x+y*y+x*z");
    System.out.println("xmin="+xmin+"  xmax="+xmax+"  dx="+dx);
    System.out.println("ymin="+ymin+"  ymax="+ymax+"  dy="+dy);
    System.out.println("zmin="+zmin+"  zmax="+zmax+"  dz="+dz);
    System.out.println("integral f(x,y,z) dx dy dz =");
    System.out.println("         1/4 x*x*z*z*y + 1/3 y*y*y*x*z + x*x*z*y");
    v = (xmax-xmin)*(ymax-ymin)*(zmax-zmin);
    System.out.println("volume ="+v+"  dv="+dv);
    System.out.println(" ");

    for(int n=3; n<7; n++)
    {
      new gaulegf(xmin, xmax, xx, wx, n);
      new gaulegf(ymin, ymax, yy, wy, n);
      new gaulegf(zmin, zmax, zz, wz, n);
      val = 0.0;
      for(int i=1; i<=n; i++)
      {
        for(int j=1; j<=n; j++)
	{
       	  for(int k=1; k<=n; k++)
	  {
            val = val + wx[i]*wy[j]*wz[k]*f(xx[i],yy[j],zz[k]);
	  } // end k
	} // end j
      } // end i
      System.out.println("n="+n+"  val ="+val);
    } // end n
    System.out.println(" ");
    val = intsum();
    System.out.println("intsum()= "+val);
    val = intfxyz();
    System.out.println("intfzyz()= "+val);
    val = intfxyz2();
    System.out.println("intfzyz2()= "+val);
    val = 1.24;
    System.out.println("int gaulegf3D= "+val);
    System.out.println(" ");
    System.out.println("test_gaulegf3D.java finished ");
  } // end test_gaulegf3D

  double f(double x, double y, double z)
  {
    return x*z + y*y + 2.0*x;
  }

  // solution check

  double fxyz(double x, double y, double z) // integral f(x,y,z)
  {
    return x*x*z*z*y/4.0 + y*y*y*x*z/3.0 + x*x*z*y;
  }

  double intfxyz() // fxyz(xmax,ymax,zmax) - fxyz(xmin,ymin.zmin)
  {
      return fxyz(xmax,ymax,zmax)-fxyz(xmin,ymin,zmin);
  }

  double intfxyz2() // fxyz(xmax-xmin,ymax-ymin,zmax-zmin)
  {
      return fxyz(xmax-xmin,ymax-ymin,zmax-zmin);
  }

  double intsum() // just sum f(x,y,z) in each dx dy dz center
  {
    double sum = 0.0;
    for(int i=0; i<nx-1; i++) // vs nx-1
    {
      x1 = xmin+dx*(double)i;
      x2 = xmin+dx*(double)(i+1);
      for(int j=0; j<ny-1; j++)
      {
        y1 = ymin+dy*(double)j;
        y2 = ymin+dy*(double)(j+1);
        for(int k=0; k<nz-1; k++)
        {
          z1 = zmin+dz*(double)k;
          z2 = zmin+dz*(double)(k+1); // should be local dv if nonuniform
          sum = sum + f((x1+x2)/2.0,(y1+y2)/2.0,(z1+z2)/2.0)*dv;
	} // end k
      } // end j
    } // end i
    return sum;
  } // end intsum

  public static void main (String[] args)
  {
    new test_gaulegf3D();
  }
} // end test_gaulegf3D