// integrate_volume.java
//  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 integrate_volume
{
  // problem statement
  double f(double x, double y, double z)
  {
    return x*z + y*y + 2.0*x;
  }

  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 = 7;
  int nz = 8;

  // for checking, e.g. exact calculus solution
  double intfxx;
  double intfxy;
  double intfxz;
  double intfyy;
  double intfyz;
  double intfzz;
  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;
  double diff = 0;
    
  public integrate_volume()
  {
    System.out.println("integrate_volume.java running");
    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");
    System.out.println("integral gaulegf3D = "+1.2514424999999632);
    System.out.println(" ");

    System.out.println("given F(x) = integral f(x) dx ");
    System.out.println("integral xmin to xmax f(x) dx = F(xmax)-F(xmin)");
    System.out.println("integral f(x,y,z) dx dy dz = ");
    System.out.println("integral (F(xmax)-F(xmin)) dy dz = ");
    System.out.println("integral (F(ymax)-F(ymin)) dz = ");
    System.out.println("F(zmax)-F(zmin)");
    System.out.println(" ");

    System.out.println("xmin="+xmin+", xmax="+xmax+", nx="+nx);
    System.out.println("ymin="+ymin+", ymax="+ymax+", ny="+ny);
    System.out.println("zmin="+zmin+", zmax="+zmax+", nx="+nz);
    v = (xmax-xmin)*(ymax-ymin)*(zmax-zmin);
    System.out.println("volume ="+v+"  dv="+dv);
    System.out.println(" ");

    vchk = 0.0;
    intchk = 0.0;
    x1 = xmin;
    x2 = xmin+dx;
    for(int i=0; i<nx-1; i++)
    {
      y1 = ymin;
      y2 = ymin+dy;
      for(int j=0; j<ny-1; j++)
      {
	z1 = zmin;
	z2 = zmin+dz;
        for(int k=0; k<nz-1; k++)
        {
	  vchk = vchk + (x2-x1)*(y2-y1)*(z2-z1);
          intchk = intchk + dx*dy*dz*(
		   (f(x1,y1,z1)+f(x2,y1,z1)+
                    f(x1,y2,z1)+f(x2,y2,z1)+
                    f(x1,y1,z2)+f(x2,y1,z2)+
                    f(x1,y2,z2)+f(x2,y2,z2))/8.0);
          intchk1 = intchk1 + dx*dy*dz*f(x1,y1,z1);
          intchk2 = intchk2 + dx*dy*dz*f(x2,y2,z2);
          z1 = z2;
          z2 = z1+dz;
	} // k
        y1 = y2;
        y2 = y1+dy;
      } // j
      x1 = x2;
      x2 = x1+dx;
    } // i
    System.out.println("volume check="+vchk);
    err = Math.abs(vchk-v);
    System.out.println("volume error = "+err);
    System.out.println(" ");

    System.out.println("integral check="+intchk);
    System.out.println("integral check1="+intchk1);
    System.out.println("integral check2="+intchk2);
    System.out.println(" ");

    System.out.println("f(xmin,ymin,zmin)="+f(xmin,ymin,zmin));
    System.out.println("f(xmax,ymin,zmin)="+f(xmax,ymin,zmin));
    System.out.println("f(xmin,ymax,zmin)="+f(xmin,ymax,zmin));
    System.out.println("f(xmin,ymin,zmax)="+f(xmin,ymin,zmax));
    System.out.println("f(xmax,ymax,zmax)="+f(xmax,ymax,zmax));
    System.out.println(" ");
    System.out.println("f(x) = 2*x + y*y + x*z");
    System.out.println("F(x) = x*x + y*y*x + x*x*z/2");
    System.out.println(" ");

    intfxx = Fxx(xmax,0.0,0.0)-Fxx(xmin,0.0,0.0);
    intfxy = Fxy(xmax,0.0,0.0)-Fxy(xmin,0.0,0.0);
    intfxz = Fxz(xmax,0.0,0.0)-Fxz(xmin,0.0,0.0); 
    System.out.println("intfxx="+intfxx);
    System.out.println("intfxy="+intfxy);
    System.out.println("intfxz="+intfxz);
    System.out.println("F(xmax)-F(xmin)="+intfxx+"+"+intfxy+"*y*y+"+intfxz+"*z");
    System.out.println(" ");
    intfyy = Fyy(ymax,0.0)-Fyy(ymin,0.0);
    intfyz = Fyz(ymax,0.0)-Fyz(ymin,0.0); 
    System.out.println("intfyy="+intfyy);
    System.out.println("intfyz="+intfyz);
    System.out.println("F(ymax)-F(ymin)="+intfyy+"+"+intfyz+"*z");
    System.out.println(" ");

    intfzz = Fzz(zmax)-Fzz(zmin); 
    System.out.println("intfzz="+intfzz);
    System.out.println("calculus dx dy dz F(zmax)-F(zmin)="+intfzz);
    System.out.println(" ");
    System.out.println("intfxyz = (((f(xmax,ymax,zmax)-f(xmax,ymax,zmin))+");
    System.out.println("            (f(xmax,ymin,zmax)-f(xmax,ymin,zmin)))-");
    System.out.println("           ((f(xmin,ymax,zmax)-f(xmin,ymax,zmin))+");
    System.out.println("            (f(xmin,ymin,zmax)-f(xmin,ymin,zmin))))");
    System.out.println("intfxyz = "+intfxyz());
    System.out.println(" ");
    System.out.println("f(xmax,ymax,zmax)="+f(xmax,ymax,zmax));
    System.out.println("f(xmax,ymax,zmin)="+f(xmax,ymax,zmin));
    diff=(f(xmax,ymax,zmax)-f(xmax,ymax,zmin));
    System.out.println("diff=(f(xmax,ymax,zmax)-f(xmax,ymax,zmin))="+diff);
    System.out.println("f(xmax,ymin,zmax)="+f(xmax,ymin,zmax));
    System.out.println("f(xmax,ymin,zmin)="+f(xmax,ymin,zmin));
    diff=(f(xmax,ymin,zmax)-f(xmax,ymin,zmin));
    System.out.println("diff=(f(xmax,ymin,zmax)-f(xmax,ymin,zmin))="+diff);
    System.out.println("f(xmin,ymax,zmax)="+f(xmin,ymax,zmax));
    System.out.println("f(xmin,ymax,zmin)="+f(xmin,ymax,zmin));
    diff=(f(xmin,ymax,zmax)-f(xmin,ymax,zmin));
    System.out.println("diff=(f(xmin,ymax,zmax)-f(xmin,ymax,zmin))="+diff);
    System.out.println("f(xmin,ymin,zmax)="+f(xmin,ymin,zmax));
    System.out.println("f(xmin,ymin,zmin)="+f(xmin,ymin,zmin));
    diff=(f(xmin,ymin,zmax)-f(xmin,ymin,zmin));
    System.out.println("diff=(f(xmin,ymin,zmax)-f(xmin,ymin,zmin))="+diff);
    System.out.println(" ");
    System.out.println("intsum2 sum centers dv/n ");
    System.out.println("intsum2 = "+intsum2());
    System.out.println(" ");
    System.out.println("intsum sum at grid dv/n ");
    System.out.println("intsum = "+intsum());
    System.out.println(" ");
    System.out.println("inttrap trapezoid method ");
    System.out.println("inttrap= "+inttrap());
    System.out.println(" ");
    System.out.println("gaulegf3D method ");
    System.out.println("gauleg3D= "+1.2514424999999632);
    System.out.println(" ");   
    System.out.println("integrate_volume.java finished");
  } // end integrate_volume

  // solution check

  double intfxyz() //  (((f(xmax,ymax,zmax)-f(xmax,ymax,zmin))+
                   //    (f(xmax,ymin,zmax)-f(xmax,ymin,zmin)))-
                   //   ((f(xmin,ymax,zmax)-f(xmin,ymax,zmin))+
                   //    (f(xmin,ymin,zmax)-f(xmin,ymin,zmin))))
  {
    double q;
    q = (((f(xmax,ymax,zmax)-f(xmax,ymax,zmin))+
          (f(xmax,ymin,zmax)-f(xmax,ymin,zmin)))-
         ((f(xmin,ymax,zmax)-f(xmin,ymax,zmin))+
          (f(xmin,ymin,zmax)-f(xmin,ymin,zmin))));
    return q;
  } // intfxyz
    
  double Fxx(double x,double y,double z)
  {
    //  f(x) = 2*x + y*y + x*z
    //  F(x) = x*x + y*y*x + x*x*z/2.0;
    return x*x; // constant
  }

  double Fxy(double x,double y,double z)
  {
    //  f(x) = 2*x + y*y + x*z
    //  F(x) = x*x + y*y*x + x*x*z/2.0;
    return x; // coefficient of y*y
  }

  double Fxz(double x,double y,double z)
  {
    //  f(x) = 2*x + y*y + x*z
    //  F(x) = x*x + y*y*x + x*x*z/2.0;
    return x*x/2.0; // coefficient of z
  }

  double Fyy(double y, double z)
  {
    // f(y) = intfxx +intfxy*y + intfxz*z
    return intfxx*y + intfxy*y*y/2.0; // constant
  }

  double Fyz(double y, double z)
  {
    // f(y) = intfxx +intfxy*y + intfxz*z
    return intfxz*y; // coefficient of z
  }

  double Fzz(double z)
  {
    // f(z) = intfyy + intfyz*z
    return intfyy*z + intfyz*z*z/2.0; // intrgral
  }

  double fxyz(double x, double y, double z) // analytic 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 intsum2() // 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
    System.out.println("last x center point = "+((x1+x2)/2.0));
    System.out.println("last x        point = "+x2+"  dx="+dx);
    System.out.println("last y center point = "+((y1+y2)/2.0));
    System.out.println("last y        point = "+y2+"  dy="+dy);
    System.out.println("last z center point = "+((z1+z2)/2.0));
    System.out.println("last z        point = "+z2+"  dz="+dz);
    System.out.println(" ");
    x1 = xmin;
    x2 = xmin+dx;
    System.out.println("first x        point = "+x1+"  dx="+dx);
    System.out.println("first x center point = "+((x1+x2)/2.0));
    y1 = ymin;
    y2 = ymin+dx;
    System.out.println("first y        point = "+y1+"  dy="+dy);
    System.out.println("first y center point = "+((y1+y2)/2.0));
    z1 = zmin;
    z2 = zmin+dx;
    System.out.println("first z        point = "+z1+"  dz="+dz);
    System.out.println("first z center point = "+((z1+z2)/2.0));
    System.out.println(" ");
    return sum;
  } // end intsum2

  double intsum() // just sum at points integrate f(x,y,z)
  {
    int n = 0;
    double sum = 0.0;
    double vol=(xmax-xmin)*(ymax-ymin)*(zmax-zmin);
    for(int i=0; i<nx; i++) // not boundary
    {
      x1 = xmin+dx*(double)i;
      for(int j=0; j<ny; j++)
      {
	y1 = ymin+dy*(double)j;
        for(int k=0; k<nz; k++)
        {
          z1 = zmin+dz*(double)k;
          sum = sum + f(x1,y1,z1);
	  n = n + 1;
	} // end k
      } // end j
    } // end i
    sum = sum*vol/(double)n;
    return sum;
  } // end intsum

  double inttrap() // trap sum f(x,y,z)
  {
    int n = 0;
    double sum = 0.0;
    double vol=(xmax-xmin)*(ymax-ymin)*(zmax-zmin);
    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);
          sum = sum + (f(x1,y1,z1)+f(x2,y1,z1))/2.0; // x  8 mid points
          sum = sum + (f(x1,y2,z1)+f(x2,y2,z1))/2.0;
          sum = sum + (f(x1,y1,z2)+f(x2,y1,z2))/2.0;
          sum = sum + (f(x1,y2,z2)+f(x2,y2,z2))/2.0;
          sum = sum + (f(x1,y1,z1)+f(x1,y2,z1))/2.0; // y
          sum = sum + (f(x2,y1,z1)+f(x2,y2,z1))/2.0;
          sum = sum + (f(x1,y1,z2)+f(x1,y2,z2))/2.0;
          sum = sum + (f(x2,y1,z2)+f(x2,y2,z2))/2.0;
          sum = sum + (f(x1,y1,z1)+f(x1,y1,z2))/2.0; // z
          sum = sum + (f(x2,y1,z1)+f(x2,y1,z2))/2.0; 
          sum = sum + (f(x1,y2,z1)+f(x1,y2,z2))/2.0; 
          sum = sum + (f(x2,y2,z1)+f(x2,y2,z2))/2.0;
	  n = n+12;
	} // end k
      } // end j
    } // end i
    sum = sum*vol/(double)n;
    return sum;
  } // end inttrap
    
  public static void main (String[] args)
  {
    new integrate_volume();
  }
} // end integrate_volume.java