(*GAUSSIAN TRIPLE INTEGRAL ALGORITHM
*
*  To approximate I = triple integral ( ( f(x,y,z,) dz dy dx ) )
*   with limits of integration from a to b for x, from c(x) to 
*   d(x) for y, and from alpha(x,y) to beta(x,y) for z.
*
*   INPUT:  endpoints a, b; positive intgers m, n, p. (Assume that
*            the roots r(i,j) and coefficients c(i,j) are available
*            for equals m,n, and p and for 1 <= j <= i.
*)
TEMP = Input["This is the Gaussian Quadrature for triple integrals.\n
   Input the function F(x,y,z) in terms of x, y, and z.\n
   For example: Sqrt[x^2+y^2]*z\n"];
F[x_,y_,z_] := Evaluate[TEMP]; 
TEMP = Input["Input the function C(x) in terms of x\n"];
c[x_] := Evaluate[TEMP]; 
TEMP = Input["Input the function D(x) in terms of x\n"];
d[x_] := Evaluate[TEMP];
TEMP = Input["Input the function alpha(x,y) in terms\n
   of x and y\n"];
alpha[x_,y_] := Evaluate[TEMP];
TEMP = Input["Input the function beta(x,y) in terms\n
   of x and y\n"];
beta[x_,y_] := Evaluate[TEMP];
OK = 0;
While[OK == 0,
   A = Input["Input the lower limit of integration\n"];
   B = Input["Input the upper limit of integration\n"];
   If[A > B,
      Input["Lower limit must be less than upper limit\n
      \n
      Press 1 [enter] to continue\n"],
      OK = 1;
   ];
];
OK = 0;
While[OK == 0,
   m = Input["Input an integer M > 1 and less than or\n 
   equal to 5 that is used for the first dimension\n "];
   n = Input["Input an integer N > 1 and less than or\n 
   equal to 5 that is used for the second dimension\n "];
   p = Input["Input an integer P > 1 and less than or\n 
   equal to 5 that is used for the third dimension\n "];
   If[n <= 1 || m <= 1 || p <= 1,
      Input["Integers must be > 1.\n
      \n
      Press 1 [enter] to continue\n"],
      If[n > 5 || m > 5 || p > 5,
	 Input["Integers must be less than or equal to 5\n
	 \n
	 Press 1 [enter] to continue\n"],
	 OK = 1;
      ];
   ];
];
If[OK == 1,
   r[1,0] = .5773502692;
   r[1,1] = -r[1,0];
   co[1,0] = 1.0;
   co[1,1] = 1.0;
   r[2,0] = 0.7745966692;
   r[2,1] = 0.0;
   r[2,2] = -r[2,0];
   co[2,0] = 0.5555555556;
   co[2,1] = 0.8888888889;
   co[2,2] = co[2,0];
   r[3,0] = 0.8611363116;
   r[3,1] = 0.3399810436;
   r[3,2] = -r[3,1];
   r[3,3] = -r[3,0];
   co[3,0] = 0.3478548451;
   co[3,1] = 0.6521451549;
   co[3,2] = co[3,1];
   co[3,3] = co[3,0];
   r[4,0] = 0.9061798459;
   r[4,1] = 0.5384693101;
   r[4,2] = 0.0;
   r[4,3] = -r[4,1];
   r[4,4] = -r[4,0];
   co[4,0] = 0.2369268850;
   co[4,1] = 0.4786286705;
   co[4,2] = 0.5688888889;
   co[4,3] = co[4,1];
   co[4,4] = co[4,0];
   (* Step 1 *)
   H1 =(B-A)/2;
   H2 =(B+A)/2;
   (* Use AJ instead of J *)
   AJ = 0;
   (* Step 2 *)
   For[i = 1,i <= m,i++,
      (* Step 3 *)
      X = H1*r[m-1,i-1]+H2;
      JX = 0;
      C1 = N[c[X]];
      D1 = N[d[X]];
      K1 = (D1-C1)/2;
      K2 = (D1+C1)/2;
      (* Step 4 *)
      For[J = 1,J <= n,J++,
	 (* Step 5 *) 
	 Y = K1*r[n-1,J-1]+K2;
	 JY = 0;
	 (* Use Z1 for Beta and Z2 for Alpha *)
	 Z1 = N[beta[X,Y]];
	 Z2 = N[alpha[X,Y]];
	 L1 = (Z1-Z2)/2;
	 L2 = (Z1+Z2)/2;
	 (* Step 6 *)
	 For[K = 1,K <= p,K++,
	    Z = L1*r[p-1,K-1]+L2;
	    Q = N[F[X,Y,Z]];
	    JY = JY+co[p-1,K-1]*Q;
	 ];
	 (* Step 7 *)
	 JX = JX+co[n-1,J-1]*L1*JY;
      ];
      (* Step 8 *)
      AJ = AJ+co[m-1,i-1]*K1*JX;
   ];
   (* Step 9 *)
   AJ = AJ*H1;
   (* Step 10 *)
   Print["\n"];
   Print["The triple integral of F from ",SetPrecision[A,9]
        ," to ",SetPrecision[B,9]," is ",SetPrecision[AJ,9]];
   Print["obtained with M = ",m," , N = ",n," , and P = ",p];
];
Quit[];