function [xp,yp] = mapc2p(xc,yc)
%
% specifies the mapping to curvilinear coordinates -- should be consistent
% with mapc2p.f
%
%     # map (xc,yc) in the square [-1,1] x [-1,1] to
%     # a square domain with ncirc embedded circles.
%

global zout cout ncirc x0circ y0circ r1circ r2circ zin cin;

xp = xc;
yp = yc;


%--------------------------------
% Loop over inclusions:
%--------------------------------
for i=1:ncirc
  x0 = x0circ(i);
  y0 = y0circ(i);
  r1 = r1circ(i);
  r2 = r2circ(i);
  xc0 = abs(xc-x0);
  yc0 = abs(yc-y0);
  
  ij1 = find(max(xc0,yc0)<=r2);  % portion of grid deformed for this circle
  xc0 = (xc(ij1) - x0)/r2;       % subset of gridpoints on this portion
  yc0 = (yc(ij1) - y0)/r2;
  
  xc1 = abs(xc0);
  yc1 = abs(yc0);
  d = max(xc1,yc1);            % value on diagonal of computational grid
  d = max(d, 1e-10);           % to avoid divide by zero at center
  d = min(d, 0.99999);         % to avoid divide by zero at d=1
  d1 = d*r2/sqrt(2);
  R = sqrt(2) * d1;
  R = r1*ones(size(d1));
  %R = r1^2 ./ (r2*d);
  
  
  % modify d1 and R outside circle to morph back to square:
  ij = find(d>r1/r2);
  d1(ij) = r1/sqrt(2) + (d(ij)-r1/r2)*(r2-r1/sqrt(2))/(1-r1/r2);
  R(ij) = r1 * ((1-r1/r2) ./ (1-d(ij))).^(r2/r1 + .5);
  
  xp2 = d1./d .* xc1;
  yp2 = d1./d .* yc1;
  center = d1 - sqrt(R.^2 - d1.^2);
  
  ij = find(xc1 >= yc1);
  xp2(ij) = center(ij) + sqrt(R(ij).^2 - yp2(ij).^2);
  
  ij = find(yc1 >= xc1);
  yp2(ij) = center(ij) + sqrt(R(ij).^2 - xp2(ij).^2);
  
  xp2 = sign(xc0) .* xp2;
  yp2 = sign(yc0) .* yp2;
  
  xp(ij1) = x0 + xp2;      % grid is modified only on portion for this circle
  yp(ij1) = y0 + yp2;
  
  end

