% 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