
function [nodes,weights]=cub_disk_prodrule(ade)

% .........................................................................
% Cubature rule on the unit disk with algebraic degree of exactness equal 
% to "ade". The nodes will be in cartesian coordinates as N x 2 matrix,
% while the weights are stored as N x 1 matrix.
%
% The rule has exactly ADE equal "ade" (tested numerically).
% .........................................................................

% shifted gaussian rule: r direction (ade 1D: ade+1)
m=ceil((ade+1)/2);
abm=r_jacobi(m,0,0); yv=gauss(m,abm); % Gauss-Legendre
rnodes=(yv(:,1)+1)/2; rw=yv(:,2)/2; 

% trapezoidal rule, consider first and last node repetition
N=ade+2;
t=linspace(0,2*pi,N); t=t(1:end-1); t=t';
tw=(2*pi/(N-1))*ones(N-1,1);

% define tensorial rule
[r_mat,th_mat]=meshgrid(rnodes,t);
[r_matw,th_matw]=meshgrid(rw,tw);

x_mat=r_mat.*cos(th_mat); y_mat=r_mat.*sin(th_mat); 
w_mat=r_matw.*th_matw.*r_mat;

% assembly nodes and weights
nodes=[x_mat(:) y_mat(:)]; weights=w_mat(:);






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