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

      double precision xc,yc,xp,yp
      double precision tol, b, r

c     # Handle the center point -- it doesn't move.
      if (yc .eq. 0 .and. xc .eq. 0) then
         xp = xc
         yp = yc
         return
      endif

      if (abs(xc) .le. yc) then
c        # North sector
         call north_sector(xc,yc,xp,yp)

      elseif (abs(yc) .le. xc) then
c        # East sector
         call north_sector(yc,xc,yp,xp)

      elseif (abs(xc) .le. -yc) then
c        # South sector
         call north_sector(xc,-yc,xp,yp)
         yp = -yp

      elseif (abs(yc) .le. -xc) then
c        # West sector
         call north_sector(yc,-xc,yp,xp)
         xp = -xp

      endif

      xp = xp/sqrt(2.d0)
      yp = yp/sqrt(2.d0)

      end

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

      double precision xc,yc,xp,yp, xi, eta
      double precision alpha, pi, beta, R

      common /compi/ pi
      common /mapparm/ alpha

c     We can assume that y > 0 here, since we are in the north sector

      xi = xc
      eta = max(1.d-16,yc)

      if (eta > 1) then
         beta = 1.d0
      else
         beta = alpha
      endif

      R = sqrt(2.d0)*(beta*eta + (1.d0-beta)/eta)

      xp = xi
      yp = sqrt(R*R - xi*xi) - sqrt(R*R - eta*eta) + eta

      end


c     # -----------------------------------------------------
c     # Jacobian
      subroutine compute_jac(xc,yc,J,Jinv)
      implicit none

      double precision xc,yc,J(2,2),Jinv(2,2), G(2,2)


c     # Handle the center point -- it doesn't move.
      if (yc .eq. 0 .and. xc .eq. 0) then
         xc = 1.d-16
         yc = 1.d-16
      endif

      if (abs(xc) .le. yc) then
c        # North sector
         G(1,1) = 1.d0
         G(1,2) = 0.d0
         G(2,1) = 0.d0
         G(2,2) = 1.d0
         call compute_jac_north_sector(xc,yc,G,J,Jinv)

      elseif (abs(yc) .le. xc) then
c        # East sector
         G(1,1) = 0
         G(1,2) = 1
         G(2,1) = 1
         G(2,2) = 0
         call compute_jac_north_sector(yc,xc,G,J, Jinv)

      elseif (abs(xc) .le. -yc) then
c        # South sector
         G(1,1) = 1.d0
         G(1,2) = 0.d0
         G(2,1) = 0.d0
         G(2,2) = -1.d0
         call compute_jac_north_sector(xc,-yc,G,J,Jinv)

      elseif (abs(yc) .le. -xc) then
c        # West sector
         G(1,1) = 0.d0
         G(1,2) = 1.d0
         G(2,1) = -1.d0
         G(2,2) = 0.d0
         call compute_jac_north_sector(yc,-xc,G,J,Jinv)

      endif

      end


      subroutine compute_jac_north_sector(xc,yc,G,J,Jinv)
      implicit none

      double precision xc,yc, G(2,2), J(2,2), Jinv(2,2)
      double precision j1(2,2)
      double precision beta, xi, eta, R, R2, Reta
      double precision pi, alpha

      common /compi/ pi
      common /mapparm/ alpha

      xi = xc

c     # note that we don't need abs(yc) here because we are assuming
c     # that is only called for yc > = 0 (in the north sector)
      eta = max(1.d-16,yc)

      if (eta > 1) then
         beta = 1.d0
      else
         beta = alpha
      endif

      R = sqrt(2.d0)*(beta*eta + (1.d0-beta)/eta)
      Reta = sqrt(2.d0)*(beta - (1-beta)/eta**2)
      R2 = R*R

      J1(1,1) = 1
      J1(1,2) = 0

      J1(2,1) = -xi/sqrt(R2 - xi*xi)
      J1(2,2) = R*Reta/sqrt(R2 - xi*xi) -
     &      (R*Reta - eta)/sqrt(R2 - eta*eta) + 1

c     # J = G'*J1*G
      call matmult(G,J1,J)
      call construct_inverse(J,Jinv)

      end


c     # C = ABA
      subroutine matmult(A,B,C)
      implicit none

      double precision A(2,2),B(2,2), C(2,2)
      double precision C1(2,2), s
      integer i,j,k

      do i = 1,2
         do j = 1,2
            s = 0
            do k = 1,2
               s = s + B(i,k)*A(k,j)
            enddo
            C1(i,j) = s
         enddo
      enddo

      do i = 1,2
         do j = 1,2
            s = 0
            do k = 1,2
               s = s + A(k,i)*C1(k,j)
            enddo
            C(i,j) = s
         enddo
      enddo

      end

      subroutine construct_inverse(A,Ainv)
      implicit none

      double precision A(2,2), Ainv(2,2), det

      det = A(1,1)*A(2,2) - A(1,2)*A(2,1)
      if (det .eq. 0) then
         write(6,*) 'construct_inverse : det == 0'
         stop
      endif

      Ainv(1,1) = A(2,2)/det
      Ainv(1,2) = -A(1,2)/det
      Ainv(2,1) = -A(2,1)/det
      Ainv(2,2) = A(1,1)/det

      end
