function [xp,yp,zp] = mapc2m(xc,yc)
% map [-3,1] x [-1,1] to surface of a sphere of radius r1

ijlower = xc < -1;
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 = d.*(2-d)/sqrt(2); % mapping d to D(d) from (??)
R = 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(1 - (xp.^2 + yp.^2));
zp(ijlower) = -zp(ijlower);          % negate z in lower hemisphere
