c
c     =====================================================
      subroutine mapc2p(xc,yc,xp,yp)
c     =====================================================
c
c     # on input,  (xc,yc) is a computational grid point
c     # on output, (xp,yp) is corresponding point in physical space
c
c     # map the square [-1,1] x [-1,1] to the starlike domain
c     # using mapping from Fig 3.2(b) followed by projection to domain, 
c     # followed by interpolation with Cartesian grid,  from
c     #   Logically Rectangular Grids and Finite Volume Methods for PDEs
c     #       in Circular and Spherical Domains,
c     #   by Donna A. Calhoun, Christiane Helzel, and Randall J. LeVeque,
c     #   http://www.amath.washington.edu/~rjl/pubs/circles
c
      implicit double precision (a-h,o-z)
c
c
      xc1 = dmax1(dabs(xc), 1d-10)
      yc1 = dmax1(dabs(yc), 1d-10)

      d = dmax1(dmax1(xc1,yc1), 1.d-10)
      d1 = d / dsqrt(2.d0)
      R = 1.d0 
      center = d1 - dsqrt(R**2 - d1**2)

      if (xc1.eq.d) then
          yp = d1/d * yc1
          xp = center + dsqrt(R**2 - yp**2)
      else 
          xp = d1/d * xc1
          yp = center + dsqrt(R**2 - xp**2)
      endif

      xp1 = dsign(xp,xc)
      yp1 = dsign(yp,yc)

      theta = datan(yp1/xp1)
      r = 1.d0 + .2d0*dcos(6.d0*theta)
      xp = xp1*r
      yp = yp1*r

c     # interpolate between above mapping and Cartesian grid:
      w = dmax1(dabs(xc),dabs(yc))**2
      w = dmin1(w,1.d0)  !# for ghost cells
      xp = w*xp + (1.d0-w)*xc/dsqrt(2.d0)
      yp = w*yp + (1.d0-w)*yc/dsqrt(2.d0)

      return
      end
