function [xp,yp] = mapc2p_2d(xc,yc)

map_choice = 4;
r1 = 1;
alpha = 0;

% Set up input arguments
xi = min(abs(xc),abs(yc));
eta = max(abs(xc),abs(yc));

% Map north sector
[X,Y] = map_north_sector(xi,eta,map_choice,r1,alpha);

% Put arguments in the right place
minXY = min(abs(X),abs(Y));
maxXY = max(abs(X),abs(Y));

ns = abs(xc) < abs(yc);  % mask for north-south sectors
ew = abs(xc) > abs(yc);  % mask for east-west sectors
diag = abs(xc) == abs(yc);

xp = sign(xc).*(minXY.*(ns + diag) + maxXY.*ew);
yp = sign(yc).*(maxXY.*(ns + diag) + minXY.*ew);


% ---------------------------------------------------
function [X,Y] = map_north_sector(xi,eta,map_choice,r1,alpha)

% global alpha r1 map_choice;

eta(eta == 0) = 1e-13;
xi(xi == 0) = 1e-13;

switch map_choice
case 1
  % radial
  % circle1
  r = sqrt(xi.^2 + eta.^2)/sqrt(2);
  X = eta.*xi./r;
  Y = eta.*eta./r;

case 2
  % linear
  % circle2
  % Radius of curvature depends linearly on eta
  % This leads to concentric circles
  R = sqrt(2)*eta;
  X = xi;
  Y = sqrt(R.^2 - xi.^2) - sqrt(R.^2 - eta.^2) + eta;

case 3
  % constant
  % grids/circle3
  % Constant radius of curvature
  R = sqrt(2);
  X = xi;
  Y = sqrt(R.^2 - xi.^2) - sqrt(R.^2 - eta.^2) + eta;

case 4
  % convex
  % grids/circle4
  % Radius of curvature is a convex combination
  r = sqrt(xi.^2 + eta.^2)/sqrt(2);
  x = eta.*xi./r;
  y = eta.*eta./r;
  w = eta.^2;
  X = w.*x + (1-w).*xi;
  Y = w.*y + (1-w).*eta;

case 5
  % Donna's grid (parametrized by alpha)
  beta = zeros(size(xi)) + alpha;
  beta(eta > 1) = 1;
  R = sqrt(2)*(beta.*eta + (1-beta)./eta);
  X = xi;
  Y = sqrt(R.^2 - xi.^2) - sqrt(R.^2 - eta.^2) + eta;

case 6
  % Donna's new grid.
  beta = sqrt(xi.^2 + eta.^2);
  % beta = eta;
  beta(eta > 1) = 1;
  R = sqrt(2)*(beta.*eta + (1-beta)./eta);
  X = xi;
  Y = sqrt(R.^2 - xi.^2) - sqrt(R.^2 - eta.^2) + eta;

case 7
  % Pathological, but simple, mapping!
  % Watch what the ghost cells do...
  X =  xi.*sqrt(2-eta.^2);
  Y = eta.*sqrt(2 -xi.^2);

end;

% Scale the circle so it has radius r1.
sf = r1./sqrt(2);  % Scale factor
X = sf*X;
Y = sf*Y;
