function [xp,yp] = mapc2p(xc,yc)
%
% map [-1,1] x [-1,1] to starlike domain as described in Sec 8.3 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

% Based on mapping from Figure 3.2(b) followed by projection to starlike
% domain, followed by interpolation with Cartesian grid.  


r1 = 1;                    

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/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);

xp1 = sign(xc) .* xp;
yp1 = sign(yc) .* yp;

xp1(xp1==0) = 1e-10;
theta = atan(yp1./xp1);
r = 1 + .2*cos(6*theta);
xp = xp1.*r;
yp = yp1.*r;

w = max(abs(xc),abs(yc)).^2;
xp = w.*xp + (1-w).*xc/sqrt(2);
yp = w.*yp + (1-w).*yc/sqrt(2);

