
function [nodes,weights]=cub_triangle_prodrule(ade)

%--------------------------------------------------------------------------
% Object
%--------------------------------------------------------------------------
% Cubature rule of product type on the unit simplex 
% * with vertices (0,0), (1,0), (0,1),
% * with ADE equal to ade.
%--------------------------------------------------------------------------




% Gaussian-Jacobi rule.
m=ceil((ade+1)/2);
ab_GJ=r_jacobi(m,1,0); xw_GJ=gauss(m,ab_GJ); % Gauss-Jacobi
s=xw_GJ(:,1); ws=xw_GJ(:,2);

% Gaussian-Legendre rule.
ab_GL=r_jacobi(m,0,0); xw_GL=gauss(m,ab_GL); % Gauss-Legendre
t=xw_GL(:,1); wt=xw_GL(:,2);

% Define tensorial rule
[S,T]=meshgrid(s,t); S=S(:); T=T(:);
nodes=[(S+1)/2 (1-S).*(1+T)/4]; 

[WS,WT]=meshgrid(ws,wt);
weights=(1/8)*WS(:).*WT(:);




% .......................... Gauss-Legendre ...............................

% GAUSS Gauss quadrature rule.
%
%    Given a weight function w encoded by the nx2 array ab of the
%    first n recurrence coefficients for the associated orthogonal
%    polynomials, the first column of ab containing the n alpha-
%    coefficients and the second column the n beta-coefficients,
%    the call xw=GAUSS(n,ab) generates the nodes and weights xw of
%    the n-point Gauss quadrature rule for the weight function w.
%    The nodes, in increasing order, are stored in the first
%    column, the n corresponding weights in the second column, of
%    the nx2 array xw.
%
function xw=gauss(N,ab)
N0=size(ab,1); if N0<N, error('input array ab too short'), end
J=zeros(N);
for n=1:N, J(n,n)=ab(n,1); end
for n=2:N
    J(n,n-1)=sqrt(ab(n,2));
    J(n-1,n)=J(n,n-1);
end
[V,D]=eig(J);
[D,I]=sort(diag(D));
V=V(:,I);
xw=[D ab(1,2)*V(1,:)'.^2];


function ab=r_jacobi(N,a,b)

nu=(b-a)/(a+b+2);
mu=2^(a+b+1)*gamma(a+1)*gamma(b+1)/gamma(a+b+2);
if N==1
    ab=[nu mu]; return
end

N=N-1;
n=1:N;
nab=2*n+a+b;
nuadd=(b^2-a^2)*ones(1,N)./(nab.*(nab+2));
A=[nu nuadd];
n=2:N;
nab=nab(n);
B1=4*(a+1)*(b+1)/((a+b+2)^2*(a+b+3));
B=4*(n+a).*(n+b).*n.*(n+a+b)./((nab.^2).*(nab+1).*(nab-1));
abadd=[mu; B1; B'];
ab=[A' abadd];