function [xp,yp] = mapc2p(xc,yc)
%
% Map the square [-1,1] x [-1,1] to a starlike domain.
% Radial projection interpolated with Cartesian near center.
%
% Mapping for Figure 8.4(a) as described in Section 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

xc1 = max(abs(xc),1e-10);
yc1 = max(abs(yc),1e-10);
alpha = min(1./xc1, 1./yc1);
beta = alpha.*sqrt(xc1.^2 + yc1.^2);
xp1 = xc./beta;
yp1 = yc./beta;

theta = atan(yc1./xc1);
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);

