      subroutine mapc2p(xc,yc,xp,yp)
      implicit none

      double precision xc,yc,xp,yp
      double precision xi,eta, X, Y
      double precision minXY,maxXY

      logical ns,ew, diag

c     Set up input arguments
      xi = min(abs(xc),abs(yc))
      eta = max(abs(xc),abs(yc))

c     Map north sector
      call map_north_sector(xi,eta,X,Y)

c     Put arguments in the right place
      minXY = min(abs(X),abs(Y))
      maxXY = max(abs(X),abs(Y))

      ns = abs(xc) .lt. abs(yc)
      ew = abs(xc) .gt. abs(yc)
      diag = abs(xc) .eq. abs(yc)

      if (ns .or. diag) then
         xp = sign(1.d0, xc)*minXY
         yp = sign(1.d0, yc)*maxXY
      else
         xp = sign(1.d0, xc)*maxXY
         yp = sign(1.d0, yc)*minXY
      endif

      end

c     ---------------------------------------------------
      subroutine map_north_sector(xi,eta,X,Y)
      implicit none

      double precision xi,eta,X,Y, beta
      double precision r, x1, y1, w, rs, theta, sf

      double precision alpha, r1

      integer map_choice

      common /mapparm/ alpha, r1
      common /mapchoice/ map_choice


      eta = max(eta, 1e-15)
      xi = max(xi, 1e-15)

      if (map_choice .eq. 1) then
c        # Radial projection grid
         r = sqrt(xi**2 + eta**2)/sqrt(2.d0)
         X = eta*xi/r
         Y = eta*eta/r
      elseif (map_choice .eq. 2) then
c        # circle2
c        # Radius of curvature depends linearly on eta
c        # This leads to concentric circles (same as
c        # map_choice == 5, with alpha = 1.
         R = sqrt(2.d0)*eta
         X = xi
         Y = sqrt(R**2 - xi**2) - sqrt(R**2 - eta**2) + eta
      elseif (map_choice .eq. 3) then
c        # grids/circle3
c        Constant radius of curvature
         R = sqrt(2.d0)
         X = xi
         Y = sqrt(R**2 - xi**2) - sqrt(R**2 - eta**2) + eta
      elseif (map_choice .eq. 4) then
c        # circle4
c        Radius of curvature is a convex combination
         r = sqrt(xi**2 + eta**2)/sqrt(2.d0)
         x1 = eta*xi/r
         y1 = eta*eta/r
         w = eta**2
         X = w*x1 + (1-w)*xi
         Y = w*y1 + (1-w)*eta
      elseif (map_choice .eq. 5) then
c        # Donna's grid (parametrized by alpha)
         if (eta .gt. 1)then
            beta = 1.d0
         else
            beta = alpha
         endif
         R = sqrt(2.d0)*(beta*eta + (1-beta)/eta)
         X = xi
         Y = sqrt(R**2 - xi**2) - sqrt(R**2 - eta**2) + eta
      elseif (map_choice .eq. 6) then
c        # Donna's new grid.
         if (eta .gt. 1) then
            beta = 1.d0
         else
            beta = sqrt(xi**2 + eta**2)/sqrt(2.d0)
         endif
         R = sqrt(2.d0)*(beta*eta + (1-beta)/eta)
         X = xi
         Y = sqrt(R**2 - xi**2) - sqrt(R**2 - eta**2) + eta
      elseif (map_choice .eq. 7) then
c        # Pathological, but simple, mapping !
         X =  xi*sqrt(2-eta**2)
         Y = eta*sqrt(2 -xi**2)
      else
         write(6,*) 'map_choice not valid'
         stop
      endif

      sf = r1/sqrt(2.d0)
      X = sf*X
      Y = sf*Y

      end
