(* GAUSSIAN DOUBLE INTEGRAL ALGORITHM 4.5
*
* To approximate I = double integral ( ( f(x,y) dy dx ) ) with
*   limits of integration from a to b for x and from 
*   c(x) to d(x) for y:
*
*    INPUT:   endpoints a, b; positive integers m, n.
*        (assume that the roots r(i,j) and coefficients
*         c(i,j) are available for i equals m and n for 1 <= j <= i
*
*   OUTPUT:   approximation J to I.
*)
TEMP = Input["This is the Gaussian Quadrature for double integrals.\n
   Input the function F(X,Y) in terms of x and y.\n
   \n
   For example: Sqrt[x^2+y^2]\n"];
F[x_,y_] := 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];
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 outer integral\n "];
   n = Input["Input an integer N > 1 and less than or\n 
   equal to 5 that is used for the inner integral\n "];
   If[n <= 1 || m <= 1,
      Input["Integers must be > 1.\n
      \n
      Press 1 [enter] to continue\n"],
      If[n > 5 || m > 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 in place 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++,
	 Y = K1*r[n-1,J-1]+K2;
	 Q = N[F[X,Y]];
	 JX = JX+co[n-1,J-1]*Q;
      ];
      (* Step 5 *)
      AJ = AJ+co[m-1,i-1]*K1*JX;
   ];
   (* Step 6 *)
   AJ = AJ*H1;
   (* Step 7 *)
   Print["\n"];
   Print["The double integral of F from ",SetPrecision[A,9]," to ",
        SetPrecision[B,9]," is ",SetPrecision[AJ,9]];
   Print["obtained with M = ",m," and N = ",n];
];
Quit[];