function D=map2D(X,Y,D);
% a function the does a numerical mapping on a 2D domain
% D.x, D.y, D.xx and D.yy are the differentiation matrices on the original
% cartesian grid
% and X, Y are the mapped domain (the domain for which you whish to have
% the differentiation matrices given as output

Nx=size(X,2);
Ny=size(X,1);

% Derivatives of the mapping
xx=D.x*X(:); xy=D.y*X(:);   yx=D.x*Y(:);  yy=D.y*Y(:);
yyy=D.yy*Y(:);   xxx=D.xx*X(:);   yxy=D.y*yx;  yxx=D.xx*Y(:);   xyy=D.yy*X(:);   xyx=D.x*xy;
jac=xx.*yy-xy.*yx;

% diff matrices
DMx=spd(yy./jac)*D.x-spd(yx./jac)*D.y;
DMy=-spd(xy./jac)*D.x+spd(xx./jac)*D.y;

DMxx=spd((yy./jac).^2)*D.xx ...
     +spd((yx./jac).^2)*D.yy ...
     +spd((yy.^2.*yxx-2*yx.*yy.*yxy+yx.^2.*yyy)./jac.^3)*(spd(xy)*D.x-spd(xx)*D.y) ...
     +spd((yy.^2.*xxx-2*yx.*yy.*xyx+yx.^2.*xyy)./jac.^3)*(spd(yx)*D.y-spd(yy)*D.x);
DMxx2=-2*spd(yx.*yy./jac.^2)*(D.y*D.x);

DMyy=spd((xy./jac).^2)*D.xx ...
     +spd((xx./jac).^2)*D.yy ...
     +spd((xy.^2.*yxx-2*xx.*xy.*yxy+xx.^2.*yyy)./jac.^3)*(spd(xy)*D.x-spd(xx)*D.y) ...
     +spd((xy.^2.*xxx-2*xx.*xy.*xyx+xx.^2.*xyy)./jac.^3)*(spd(yx)*D.y-spd(yy)*D.x);
DMyy2=-2*spd(xx.*xy./jac.^2)*(D.y*D.x);

% the differentiation matrices
D.x=DMx; D.xx=DMxx+DMxx2; 
D.y=DMy; D.yy=DMyy+DMyy2;

% integration operators
D.w=D.w(:).*jac; % double integration 
D.wx=D.wx(:).*xx; % integration in x
D.wy=D.wy(:).*yy; % integration in y

D.w=reshape(D.w,Ny,Nx);
D.wx=reshape(D.wx,Ny,Nx);
D.wy=reshape(D.wy,Ny,Nx);