% Demonstration of integration in polar coordinates
% We will integrate over the unit disc
%
% Samuli Siltanen November 2012



% Construct the quadrature points and weights,
% first in 1D
N      = 20;
[r,wg] = gaussint(N,0,1);
M      = 20;
phi    = 2*pi*[0:(M-1)]/M;
wphi   = (phi(2)-phi(1))*ones(size(phi));

% Now move to 2D
[R,Phi] = meshgrid(r,phi);
[W1,W2] = meshgrid(r(:).*wg(:),wphi);
Wmat    = W1.*W2;

% Construt the function to be integrated
%fun = ones(size(R));
fun = log(R);
fun = fun(:);

% Integrate by inner product
int = (Wmat(:).')*fun
