
function [xyw,xyw_bar]=set_stroud_pre(ade)

% INPUT:
% ade: degree of precision.
%
% OUTPUT:
% xyw determines the cubature points and weights, in BARYCENTRIC coordinates.
% nodes: xyw(:,1:3)
% weights: xyw(:,4)

xyw=stroud_conical_rules_ref(ade);
xx=xyw(:,1); yy=xyw(:,2); ww=xyw(:,3);

L2=yy;
L3=xx-L2;
L1=1-L2-L3;

ww_bar=2*ww;

xyw_bar=[L1 L2 L3 ww_bar];

xyw=xyw_bar(:,[1 3 4]);
xyw(:,3)=xyw(:,3)/2;
xyw_bar=[xyw(:,1:2) 1-xyw(:,1)-xyw(:,2) xyw(:,3)*2];







function xyw=stroud_conical_rules_ref(ade)

% SEE LYNESS, COOLS, "A survey on numerical cubature over triangles", p.4.

N=ceil((ade+1)/2);

[t,T]=gauss_jacobi(N,0,0);
t=(t+1)/2; T=T/2;

[x,w]=gauss_jacobi(N,0,1);
x=(x+1)/2; X=w/4;

[wx,wy]=meshgrid(T,X);
ww=wx.*wy;

[yt,xx]=meshgrid(x,x);

[yt,yx]=meshgrid(t,x);
yy=yt.*yx;

xyw=[xx(:) yy(:) ww(:)];






function [x,w]=gauss_jacobi(N,a,b,gl)

% GAUSS-JACOBI (LOBATTO) RULE ON (-1,1).
% N IS ...
% a,b ARE THE GAUSS-JACOBI EXPONENTS.
% gl: 0: GAUSS POINTS. 1: GAUSS-LOBATTO POINTS.
% x, w ARE COLUMN VECTORS OF NODES AND WEIGHTS.
%      THE LENGTH OF x AND w IS "N" IF gl=0, "N+2" IF "gl=1".

if nargin < 2
    a=0; b=0;
end

if nargin < 4
    gl = 0;
end

if gl == 0
    ab=r_jacobi(N,a,b);
    xw=gauss(N,ab);
else
    xw=lobatto_jacobi(N,a,b);
end

x=xw(:,1);
w=xw(:,2);




%--------------------------------------------------------------------------
% ADDITIONAL FUNCTIONS BY D.LAURIE AND W.GAUTSCHI.
%--------------------------------------------------------------------------

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


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 xw=lobatto_jacobi(N,a,b)

if nargin<2, a=0; end;
if nargin<3, b=a; end
ab=r_jacobi(N+2,a,b);
ab(N+2,1)=(a-b)/(2*N+a+b+2);
ab(N+2,2)=4*(N+a+1)*(N+b+1)*(N+a+b+1)/((2*N+a+b+1)*(2*N+a+b+2)^2);
xw=gauss(N+2,ab);


function [xx,ww]=rules_convert(xw)

% CONVERTING RULES FROM A REFERENCE SIMPLEX TO ANOTHER.

x=xw(:,1); y=xw(:,2); ww=xw(:,3);
pts_bar=[x y 1-x-y];
delta=0;
vertices=[-1+delta -1+delta; 1-2*delta -1+delta; -1+delta 1-2*delta];
xx=pts_bar*vertices;
ww=4*ww;

