function [xp,yp,zp] = mapc2p(xc,yc,zc)
%
% Map [-3,1] x [-1,1] x [0,1] to spherical shell
% Mapping for Figure 5.1 as described in Section 5.1 of
%    Logically Rectangular Grids and Finite Volume Methods for PDEs
%         in Circular and Spherical Domains,
%    by Donna A. Calhoun, Christiane Helzel, and Randall J. LeVeque,
%    http://www.amath.washington.edu/~rjl/pubs/circles


r1 = 1;                    % inner radius
dr = 0.5;                  % depth of shell
r2 = r1 + dr;              % outer radius

ijlower = find(xc<-1);           % indices of points in lower hemisphere
xc(ijlower) = -2 - xc(ijlower);

d = max(abs(xc),abs(yc));  % value on diagonal of computational grid
d = max(d, 1e-10);         % to avoid divide by zero at origin

D = r1 * d.*(2-d)/sqrt(2); % mapping d to D(d) from (??)
R = r1 * ones(size(d));    % mapping d to R(d) from (??) 

center = D - sqrt(R.^2 - D.^2);
xp = D./d .* abs(xc);
yp = D./d .* abs(yc);

ij = find(abs(yc)>=abs(xc));
yp(ij) = center(ij) + sqrt(R(ij).^2 - xp(ij).^2);
ij = find(abs(xc)>=abs(yc));
xp(ij) = center(ij) + sqrt(R(ij).^2 - yp(ij).^2);

xp = sign(xc) .* xp;
yp = sign(yc) .* yp;
zp = sqrt(r1^2 - (xp.^2 + yp.^2));
zp(ijlower) = -zp(ijlower);          % negate z in lower hemisphere
Rz = r1 + zc*dr;                % radius based on zc
xp = Rz.*xp;
yp = Rz.*yp;
zp = Rz.*zp;


