function [xp,yp] = mapc2p(xc,yc)
% map [-1,1] x [-1,1] to circle of radius r1
% Mapping from Section 6.2

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 (??)
D = r1*d.*(2-d)/sqrt(2);
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;
