
function [nodes,w]=gaussian_product_rule(n)

% Gaussian product rule of degree 2n-1.
% Nodes are in cartesian coordinates.

% Trapezoidal rule
m=2*n; h=(2*pi/m); phi=(1:m)'*h;

% Gauss-Legendre rule
ab=r_jacobi(n,0,0); xw=gauss(n,ab);

% Nodes
z=xw(:,1); 
[P,T]=meshgrid(phi,sqrt(1-z.^2));
[PZ,Z]=meshgrid(phi,z);
P=P(:); T=T(:); Z=Z(:);
nodes=[cos(P).*T sin(P).*T Z];

% Weights
hv=h*ones(m,1);
[W1,W2]=meshgrid(hv,xw(:,2));
w=W1(:).*W2(:);

