LineGaussRuleInfo[ {rule_, numer_ }, point_ ] := Module [
{ g2 = {-1,1}/Sqrt[3],
w3={5/9,8/9,5/9},
g3={-Sqrt[3/5],0,Sqrt[3/5]},
w4={(1/2)-Sqrt[5/6]/6,
(1/2)+Sqrt[5/6]/6,
(1/2)+Sqrt[5/6]/6,
(1/2)-Sqrt[5/6]/6},
g4={-Sqrt[(3+2*Sqrt[6/5])/7],
-Sqrt[(3-2*Sqrt[6/5])/7],
Sqrt[(3-2*Sqrt[6/5])/7],
Sqrt[(3+2*Sqrt[6/5])/7]},
g5={-Sqrt[5+2*Sqrt[10/7]],
-Sqrt[5-2*Sqrt[10/7]],
0,
Sqrt[5-2*Sqrt[10/7]],
Sqrt[5+2*Sqrt[10/7]]}/3,
w5={322-13*Sqrt[70],
322+13*Sqrt[70],
512,
322+13*Sqrt[70],
322-13*Sqrt[70]}/900,
i=point,p=rule, info={Null,0}},
If [ p == 1,info = {0,2}];
If [ p == 2, info={g2[[i]],1}];
If [ p == 3, info={g3[[i]],w3[[i]]}];
If [ p == 4, info ={g4[[i]],w4[[i]]}];
If [ p == 5, info = {g5[[i]],w5[[i]]}];
If [ numer, Return [ N[info,20]], Return [ Simplify[info]]];
];

QuadGaussRuleInfo [ { rule_, numer_ }, point_ ] := Module [
  {xi, eta, p1, p2, i, j, w1, w2, m, info = { { Null, Null}, 0 }},
  If [ Length[rule] == 2, { p1,p2} = rule, p1=p2=rule];
  If [ Length[point]==2, {i,j}=point, m = point;
    j=Floor[(m-1)/p1]+1;i=m-p1*(j-1)];
  {xi,w1} = LineGaussRuleInfo[{p1,numer},i];
{eta,w2} = LineGaussRuleInfo[{p2,numer},j];
  info={{xi,eta},w1*w2};
  If [ numer, Return[N[info,20]], Return[Simplify[info]]];
];

Print[" "];
Print["QuadGaussRuleInfoTest.nb"];
Print["  Mathematica version"];
Print[" "];
Print["  Test the quadrature rules for a unit quadrilateral,"];
Print["  returned as {{x,y},w} = QuadGaussRuleInfo[{rule,numer},i]"];
Print[" "];
Print["  INPUT:"];
Print[" "];
Print["  RULE is 1, 2, 3, 4, 5, the order in both dimensions,"];
Print["  or RULE = { order1, order2 }, allowing the user to specify"];
Print["  different orders in each dimension."];
Print[" "];
Print["  NUMER is True to return numerical values for the output, or"];
Print["  False to return symbolic values."];
Print[" "];
Print["  I is the index, between 1 and RULE*RULE, of the point."]
Print[" The index can also be specified as {I,J}."]
Print[" "];
Print["  OUTPUT:"];
Print[" "];
Print["  {X,Y} are the coordinates of the abscissas,"]
Print["  For the unit quadrilateral, -1 <= X, Y <= 1."]
Print[" "];
Print["  W is the weight of the quadrature point."]
Print[" "]
Print["  The following cells print out each rule, and the sum"]
Print["  of the weights."]

For [ rule = 1, rule ≤ 5, rule++,
numer = True;
wsum = 0;
Print[" "];
Print["Quadrilateral rule number ", rule];
Print[" "];
Print["W"];
Print[" "];
For [ i=1, i ≤ rule * rule, i++,
{ { x, y},w} = QuadGaussRuleInfo[{rule,numer},i];
Print [ w ];
wsum = wsum +w;
];
Print["sum(w) = ", wsum ];
Print[" "];
Print["X, Y"];
Print[" "];
For [ i=1, i ≤ rule * rule, i++,
{ { x, y },w} = QuadGaussRuleInfo[{rule,numer},i];
Print [ x,"  ", y ];
];
];

rule1 = 2;
rule2 = 2;
numer = True;
wsum = 0;
Print[" "];
Print["Quadrilateral rule number {", rule1, ",", rule2, "}"];
Print[" "];
Print["W"];
Print[" "];
For [ i=1, i ≤ rule1 * rule2, i++,
{ {x,y},w} = QuadGaussRuleInfo[{{rule1,rule2},numer},i];
Print [ w ];
wsum = wsum +w;
];
Print["sum(w) = ", wsum ];
Print[" "];
Print["X, Y"];
Print[" "];
For [ i=1, i ≤ rule1 * rule2, i++,
{ {x,y},w} = QuadGaussRuleInfo[{rule,numer},i];
Print [x,"  ", y ];
];

rule1 = 2;
rule2 = 3;
rule = { rule1, rule2};
numer = True;
wsum = 0;
Print[" "];
Print["Quadrilateral rule number {", rule1, ",", rule2, "}"];
Print[" "];
Print["W"];
Print[" "];
For [ i = 1, i ≤ rule1, i++,
For [ j = 1, j ≤ rule2, j++, 
{ {x,y},w} = QuadGaussRuleInfo[{rule,numer},{i,j}];
Print [ w ];
wsum = wsum +w;
];
];
Print["sum(w) = ", wsum ];
Print[" "];
Print["X, Y"];
Print[" "];
For [ i = 1, i ≤ rule1, i++,
For [ j = 1, j ≤ rule2, j++, 
{ {x,y},w} = QuadGaussRuleInfo[{rule,numer},{i,j}];
Print [x,"  ", y ];
];
];