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 a starlike domain using radial
c     # projection mapping followed by interpolation  (Fig 8.4(b))
c     # as described in Section 8.3 of
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     # radial projection mapping
      r1 = 1.d0
      d = dmax1(dabs(xc), dabs(yc))
      r = dmax1(dsqrt(xc**2 + yc**2), 1.d-10)
      xp = r1*d*xc/r
      yp = r1*d*yc/r

      theta = datan(dabs(yc)/dmax1(dabs(xc),1.d-10))
      r = 1.d0 + .2d0*dcos(6.d0*theta)
      xp = xp*r
      yp = yp*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
