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

      double precision xc,yc,zc,xp,yp,zp
      double precision d(6), mind, alpha, s3
      integer i,idx

c     NOTE : This is the same mapping as the radial projection
c     mapping.


      d(1) = 1.d0 - xc
      d(2) = xc + 1.d0

      d(3) = 1.d0 - yc
      d(4) = yc + 1.d0

      d(5) = 1.d0 - zc
      d(6) = zc + 1.d0

      mind = 2.d0
      do i = 1,6
         if (d(i) .le. mind) then
            mind = d(i)
            idx = i
         endif
      enddo


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

      if (idx .eq. 1) then
c        # East sector
         call top_sector(zc,yc,xc,zp,yp,xp)

      elseif (idx .eq. 2) then
c        # West sector
         call top_sector(zc,yc,-xc,zp,yp,xp)
         xp = -xp

      elseif (idx .eq. 3) then
c        # north sector
         call top_sector(xc,zc,yc,xp,zp,yp)

      elseif (idx .eq. 4) then
c        # south sector
         call top_sector(xc,zc,-yc,xp,zp,yp)
         yp = -yp

      elseif (idx .eq. 5) then
c        # top sector
         call top_sector(xc,yc,zc,xp,yp,zp)

      elseif (idx .eq. 6) then
c        # bottom sector
         call top_sector(xc,yc,-zc,xp,yp,zp)
         zp = -zp

      endif


      s3 = 1.d0/sqrt(3.d0)
      xp = s3*xp
      yp = s3*yp
      zp = s3*zp

      end


      subroutine top_sector(xc,yc,zc,xp,yp,zp)
      implicit none

      double precision xc,yc,zc,xp,yp,zp

      double precision xi,eta,zeta,alpha,beta, R
      double precision a,b,a2,b2,R2,d,d2,z2, rr

      common /mapparm/ alpha

      xi = xc
      eta = yc
      zeta = sign(max(1.d-16,abs(zc)),zc)

      a = xi/zeta
      b = eta/zeta

      if (zeta .gt. 1.d0) then
         beta = 1.d0
      else
         beta = alpha
      endif

      R = sqrt(3.d0)*(beta*zeta + (1.d0-beta)/zeta)

      R2 = R*R
      Z2 = zeta*zeta
      a2 = a*a
      b2 = b*b

      if (R2 - 2*z2 < 0) then
         write(6,*) 'r2 - 2*z2 < 0'
         write(6,'(4E16.8)') xc,yc,zc,R**2 - 2*z2
         stop
      endif
      d = sqrt(R2 - 2*z2) - zc;
      d2 = d*d

      zp = 1.d0/(1+a2+b2)*(sqrt((1+a2+b2)*R2 - (a2+b2)*d2) - d)
      xp = a*zp
      yp = b*zp


      end
