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]]];
];

TrigGaussRuleInfo[{rule_,numer_},point_] := Module [
  {zeta,p=rule,i=point,g1,g2,g3,g4,s1,w2,w3,eps=10.^(-24),
  jkl={{1,2,3},{2,1,3},{1,3,2},{3,1,2},{2,3,1},{3,2,1}},
  info={{Null,Null,Null},0} },
  
  If [ p ==  1, info={{1/3,1/3,1/3},1}];
  
  If [ p ==  3, info={{1,1,1}/6,1/3}; info[[1,i]]=2/3];
  
  If [ p == -3, info={{1,1,1}/2,1/3}; info[[1,i]]=0];
  
  If [ p ==  6, g1 = (8-Sqrt[10]+Sqrt[38-44*Sqrt[2/5]])/18;
                g2 = (8-Sqrt[10]-Sqrt[38-44*Sqrt[2/5]])/18;
      If [ i < 4, info={{g1,g1,g1}, (620+Sqrt[213125-53320*Sqrt[10]])/3720};
          info[[1,i]]=1-2*g1];
      If [ i > 3, info={{g2,g2,g2}, (620-Sqrt[213125-53320*Sqrt[10]])/3720};
          info[[1,i-3]]=1-2*g2]];
          
  If [ p == -6,
      If[i<4, info={{1,1,1}/6,3/10}; info[[1,i]]=2/3];
      If[i>3, info={{1,1,1}/2,1/30}; info[[1,i-3]]=0]];
      
  If [ p == 7, g1=(6-Sqrt[15])/21;
               g2=(6+Sqrt[15])/21;
     If[i<4, info={{g1,g1,g1},(155-Sqrt[15])/1200};
             info[[1,i]] = 1-2*g1];
     If[i>3&&i<7, info={{g2,g2,g2},(155+Sqrt[15])/1200};
             info[[1,i-3]]=1-2*g2];
     If[i==7, info={{1/3,1/3,1/3},9/40}]];
     
  If [ p == 12,
     g1=0.063089014491502228340331602870819157;
     g2=0.249286745170910421291638553107019076;
     g3=0.053145049844816947353249671631398147;
     g4=0.310352451033784405416607733956552153;
     If [ !numer, {g1,g2,g3,g4}=Rationalize[{g1,g2,g3,g4},eps]];
      w1=(30*g2^3*(4*g3^2+(1-2*g4)^2+4*g3*(-1+g4))+
          g3^2*(1-15*g4)+(-1+g4)*g4-g3*(-1+g4)*(-1+15*g4)+
          2*g2+(1+60*g3*g4*(-1+g3+g4))-6*g2^2*(3+10*(-1+g4)*g4+
          10*g3^2*(1+3*g4)+10*g3*(-1+g4)*(1+3*g4)))/
         (180*(g1-g2)*(-(g2*(-1+2*g2)*(-1+g3)*g3)+(-1+g3)*(g2-2*g2^2-
          2*g3+3*g2*g3)*g4-(g2*(-1+2*g2-3*g3)+2*g3)*g4^2+
          2*g1^2*(g2*(-2+3*g2)+g3-g3^2+g4-g3*g4-g4^2)+
          g1*(-4*g2^2+(-1+g3)*g3+(-1+g3)*(1+3*g3)*g4+(1+3*g3)*
          g4^2-2*g2*(-1+g3^2+g3*(-1+g4)+(-1+g4)*g4))));
      w2=(-1+12*(2-3*g1)*g1*w1+4*g3^2*(-1+3*w1)+4*g3*(-1+g4)*(-1+3*w1)+
          4*(-1+g4)*g4*(-1+3*w1))/(12*(g2*(-2+3*g2)+g3-g3^2+g4-g3*g4-g4^2));
      w3=(1-3*w1-3*w2)/6;
      If [ i<4, info={{g1,g1,g1},w1}; info[[1,i]] = 1-2*g1];
      If [ i>3&&i<7, info={{g2,g2,g2},w2}; info[[1,i-3]]=1-2*g2];
      If [ i>6, {j,k,l}=jkl[[i-6]]; info={{0,0,0},w3};
        info[[1,j]]=g3; info[[1,k]]=g4;info[[1,l]]=1-g3-g4]];
        
  If [ numer, Return[N[info,20]], Return[Simplify[info]]];
];

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

Print[" "];
Print["WedgGaussRuleInfoTest.nb"];
Print["  Mathematica version"];
Print[" "];
Print["  Test the quadrature rules for a unit wedge,"];
Print["  returned as {{a,b,c,d},w} = WedgGaussRuleInfo[{rule,numer},i]"];
Print[" "];
Print["  INPUT:"];
Print[" "];
Print["  RULE is Rule1 x Rule2, where"];
Print["  RULE1 is a Triangle rule (1, 3, -3, 6, 7, 12)"];
Print["  RULE2 is a Line rule (1, 2, 3, 4, or 5)."];
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 RULE1*RULE2, of the point."]
Print[" The index can also be specified as {I,J}."]
Print[" "];
Print["  OUTPUT:"];
Print[" "];
Print["  {A,B,C,D},"]
Print["  {A,B,C} are barycentric coordinates of the triangle,"];
Print["  D is the lateral component."]
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."]

rule1 = 1;
rule2 = 1;
numer = True;
wsum = 0;
Print[" "];
Print["Wedge rule number {", rule1, ", ", rule2, "}"];
Print[" "];
Print["W"];
Print[" "];
For [ i=1, i ≤ rule1 * rule2, i++,
{ { a, b, c, d},w} = WedgGaussRuleInfo[{{rule1,rule2},numer},i];
Print [ w ];
wsum = wsum +w;
];
Print["sum(w) = ", wsum ];
Print[" "];
Print["X, Y, Z"];
Print[" "];
For [ i=1, i ≤ rule1 * rule2, i++,
{ { a,b,c,d},w} = WedgGaussRuleInfo[{{rule1,rule2},numer},i];
Print [ a,"  ", b, "  ", d];
];

rule1 = 3;
rule2 = 2;
numer = True;
wsum = 0;
Print[" "];
Print["Wedge rule number {", rule1, ", ", rule2, "}"];
Print[" "];
Print["W"];
Print[" "];
For [ i=1, i ≤ rule1 * rule2, i++,
{ { a, b, c, d},w} = WedgGaussRuleInfo[{{rule1,rule2},numer},i];
Print [ w ];
wsum = wsum +w;
];
Print["sum(w) = ", wsum ];
Print[" "];
Print["X, Y, Z"];
Print[" "];
For [ i=1, i ≤ rule1 * rule2, i++,
{ { a,b,c,d},w} = WedgGaussRuleInfo[{{rule1,rule2},numer},i];
Print [ a,"  ", b, "  ", d];
];

rule1 = -3;
rule2 = 2;
numer = True;
wsum = 0;
Print[" "];
Print["Wedge rule number {", rule1, ", ", rule2, "}"];
Print[" "];
Print["W"];
Print[" "];
For [ i=1, i ≤ abs ( rule1 ) * rule2, i++,
{ { a, b, c, d},w} = WedgGaussRuleInfo[{{rule1,rule2},numer},i];
Print [ w ];
wsum = wsum +w;
];
Print["sum(w) = ", wsum ];
Print[" "];
Print["X, Y, Z"];
Print[" "];
For [ i=1, i ≤ abs ( rule1 ) * rule2, i++,
{ { a,b,c,d},w} = WedgGaussRuleInfo[{{rule1,rule2},numer},i];
Print [ a,"  ", b, "  ", d];
];

rule1 = 6;
rule2 = 3;
numer = True;
wsum = 0;
Print[" "];
Print["Wedge rule number {", rule1, ", ", rule2, "}"];
Print[" "];
Print["W"];
Print[" "];
For [ i=1, i ≤ rule1 * rule2, i++,
{ { a, b, c, d},w} = WedgGaussRuleInfo[{{rule1,rule2},numer},i];
Print [ w ];
wsum = wsum +w;
];
Print["sum(w) = ", wsum ];
Print[" "];
Print["X, Y, Z"];
Print[" "];
For [ i=1, i ≤ rule1 * rule2, i++,
{ { a,b,c,d},w} = WedgGaussRuleInfo[{{rule1,rule2},numer},i];
Print [ a,"  ", b, "  ", d];
];

rule1 = 12;
rule2 = 5;
numer = True;
wsum = 0;
Print[" "];
Print["Wedge rule number {", rule1, ", ", rule2, "}"];
Print[" "];
Print["W"];
Print[" "];
For [ i=1, i ≤ rule1 * rule2, i++,
{ { a, b, c, d},w} = WedgGaussRuleInfo[{{rule1,rule2},numer},i];
Print [ w ];
wsum = wsum +w;
];
Print["sum(w) = ", wsum ];
Print[" "];
Print["X, Y, Z"];
Print[" "];
For [ i=1, i ≤ rule1 * rule2, i++,
{ { a,b,c,d},w} = WedgGaussRuleInfo[{{rule1,rule2},numer},i];
Print [ a,"  ", b, "  ", d];
];