function [Dx,Dy,Dxx,Dyy,ww,wx,wy]=mapping2D(X,Y,Dx,Dy,Dxx,Dyy,ww,wx,wy);
% a function the does a numerical mapping on a 2D domain
% Dx, Dy, Dxx and Dyy 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

% Derivatives of the mapping
xx=Dx*X(:); xy=Dy*X(:);   yx=Dx*Y(:);  yy=Dy*Y(:);
yyy=Dyy*Y(:);   xxx=Dxx*X(:);   yxy=Dy*yx;  yxx=Dxx*Y(:);   xyy=Dyy*X(:);   xyx=Dx*xy;
jac=xx.*yy-xy.*yx;

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

DMxx=spd((yy./jac).^2)*Dxx ...
     +spd((yx./jac).^2)*Dyy ...
     +spd((yy.^2.*yxx-2*yx.*yy.*yxy+yx.^2.*yyy)./jac.^3)*(spd(xy)*Dx-spd(xx)*Dy) ...
     +spd((yy.^2.*xxx-2*yx.*yy.*xyx+yx.^2.*xyy)./jac.^3)*(spd(yx)*Dy-spd(yy)*Dx);
DMxx2=-2*spd(yx.*yy./jac.^2)*(Dy*Dx);

DMyy=spd((xy./jac).^2)*Dxx ...
     +spd((xx./jac).^2)*Dyy ...
     +spd((xy.^2.*yxx-2*xx.*xy.*yxy+xx.^2.*yyy)./jac.^3)*(spd(xy)*Dx-spd(xx)*Dy) ...
     +spd((xy.^2.*xxx-2*xx.*xy.*xyx+xx.^2.*xyy)./jac.^3)*(spd(yx)*Dy-spd(yy)*Dx);
DMyy2=-2*spd(xx.*xy./jac.^2)*(Dy*Dx);

% the differentiation matrices
Dx=DMx; Dxx=DMxx+DMxx2; 
Dy=DMy; Dyy=DMyy+DMyy2;

% integration operators
ww=reshape(ww(:).*jac,size(X)); 
wx=reshape(wx(:).*xx,size(X)); 
wy=reshape(wy(:).*yy,size(X));