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

Print[" "];
Print["TrigGaussRuleInfoTest.nb"];
Print["  Mathematica version"];
Print[" "];
Print["  Test the quadrature rules for a triangle, returned"];
Print["  by {{a,b,c},w} = TrigGaussRuleInfo[{rule,numer},i]"];
Print[" "];
Print["  INPUT:"];
Print[" "];
Print["  RULE is 1, 3, -3, 6, -6, 7, and 12, the order"];
Print["  of the rule.  There are two rules of order 3 and 6,"];
Print["  so the negative value returns the second one."];
Print[" "];
Print["  NUMER is True to return numerical values for the output, or"];
Print["  False to return symbolic values (except for 12)."];
Print[" "];
Print["  I is the index, between 1 and |RULE|, of the point desired."]
Print[" "];
Print["  OUTPUT:"];
Print[" "];
Print["  {A,B,C} are the barycentric coordinates of the abscissas,"]
Print["  For the unit triangle, X = A, Y = B."]
Print["  with 0 <= X, Y <= 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."]

rule =1;
numer = True;
wsum = 0;
Print[" "];
Print["Triangle rule number ", rule];
Print[" "];
Print["W"];
Print[" "];
For [ i=1, i ≤ Abs[rule], i++,
{ {a,b,c},w} = TrigGaussRuleInfo[{rule,numer},i];
Print [ w ];
wsum = wsum +w;
];
Print["sum(w) = ", wsum ];
Print[" "];
Print["X, Y"];
Print[" "];
For [ i=1, i ≤ Abs[rule], i++,
{ {a,b,c},w} = TrigGaussRuleInfo[{rule,numer},i];
Print [ a,"  ", b ];
];

rule =3;
numer = True;
wsum = 0;
Print[" "];
Print["Triangle rule number ", rule];
Print[" "];
Print["W"];
Print[" "];
For [ i=1, i ≤ Abs[rule], i++,
{ {a,b,c},w} = TrigGaussRuleInfo[{rule,numer},i];
Print [ w ];
wsum = wsum +w;
];
Print["sum(w) = ", wsum ];
Print[" "];
Print["X, Y"];
Print[" "];
For [ i=1, i ≤ Abs[rule], i++,
{ {a,b,c},w} = TrigGaussRuleInfo[{rule,numer},i];
Print [ a,"  ", b ];
];

rule =-3;
numer = True;
wsum = 0;
Print[" "];
Print["Triangle rule number ", rule];
Print[" "];
Print["W"];
Print[" "];
For [ i=1, i ≤ Abs[rule], i++,
{ {a,b,c},w} = TrigGaussRuleInfo[{rule,numer},i];
Print [ w ];
wsum = wsum +w;
];
Print["sum(w) = ", wsum ];
Print[" "];
Print["X, Y"];
Print[" "];
For [ i=1, i ≤ Abs[rule], i++,
{ {a,b,c},w} = TrigGaussRuleInfo[{rule,numer},i];
Print [ a,"  ", b ];
];

rule =6;
numer = True;
wsum = 0;
Print[" "];
Print["Triangle rule number ", rule];
Print[" "];
Print["W"];
Print[" "];
For [ i=1, i ≤ Abs[rule], i++,
{ {a,b,c},w} = TrigGaussRuleInfo[{rule,numer},i];
Print [ w ];
wsum = wsum +w;
];
Print["sum(w) = ", wsum ];
Print[" "];
Print["X, Y"];
Print[" "];
For [ i=1, i ≤ Abs[rule], i++,
{ {a,b,c},w} = TrigGaussRuleInfo[{rule,numer},i];
Print [ a,"  ", b ];
];

rule =-6;
numer = True;
wsum = 0;
Print[" "];
Print["Triangle rule number ", rule];
Print[" "];
Print["W"];
Print[" "];
For [ i=1, i ≤ Abs[rule], i++,
{ {a,b,c},w} = TrigGaussRuleInfo[{rule,numer},i];
Print [ w ];
wsum = wsum +w;
];
Print["sum(w) = ", wsum ];
Print[" "];
Print["X, Y"];
Print[" "];
For [ i=1, i ≤ Abs[rule], i++,
{ {a,b,c},w} = TrigGaussRuleInfo[{rule,numer},i];
Print [ a,"  ", b ];
];

rule =7;
numer = True;
wsum = 0;
Print[" "];
Print["Triangle rule number ", rule];
Print[" "];
Print["W"];
Print[" "];
For [ i=1, i ≤ Abs[rule], i++,
{ {a,b,c},w} = TrigGaussRuleInfo[{rule,numer},i];
Print [ w ];
wsum = wsum +w;
];
Print["sum(w) = ", wsum ];
Print[" "];
Print["X, Y"];
Print[" "];
For [ i=1, i ≤ Abs[rule], i++,
{ {a,b,c},w} = TrigGaussRuleInfo[{rule,numer},i];
Print [ a,"  ", b ];
];

rule =12;
numer = True;
wsum = 0;
Print[" "];
Print["Triangle rule number ", rule];
Print[" "];
Print["W"];
Print[" "];
For [ i=1, i ≤ Abs[rule], i++,
{ {a,b,c},w} = TrigGaussRuleInfo[{rule,numer},i];
Print [ w ];
wsum = wsum +w;
];
Print["sum(w) = ", wsum ];
Print[" "];
Print["X, Y"];
Print[" "];
For [ i=1, i ≤ Abs[rule], i++,
{ {a,b,c},w} = TrigGaussRuleInfo[{rule,numer},i];
Print [ a,"  ", b ];
];