mapc2m.f.html
Source file:   mapc2m.f
Directory:   /Users/rjl/git/rjleveque/clawpack-4.6.3/apps/advection/2d/sphere/rotation
Converted:   Mon Jan 21 2013 at 20:15:21   using clawcode2html
This documentation file will not reflect any later changes in the source file.

 
c
c     ==========================================
      subroutine mapc2m(x1,y1,xp,yp,zp)
c     ==========================================
c
c     Map computational point (x1,y1) to the point (xp,yp,zp)
c     on the sphere
c
c     Points in [-1,1]x[-1,1]  are mapped to the upper hemisphere
c     Points in [-3,-1]x[-1,1] are mapped to the lower hemisphere
c
      implicit real*8(a-h,o-z)

      r1 = 1.d0
      xc = x1
      yc = y1

c     # ghost cell values outside of [-3,1]x[-1,1] get mapped to other
c     # hemisphere:
      if (xc.ge.1.d0) xc = xc-4.d0
      if (xc.lt.-3.d0) xc = xc+4.d0

      if (yc.ge.1.d0) then
         yc = 2.d0 - yc
         xc = -2.d0 - xc
         endif
      if (yc.lt.-1.d0) then
         yc = -2.d0 - yc
         xc = -2.d0 - xc
         endif

      if (xc.lt.-1.d0) then
c         # points in [-3,-1] map to lower hemisphere - reflect about x=-1
c         # to compute x,y mapping and set sgnz appropriately:
          xc = -2.d0 - xc
          sgnz = -1.d0
         else
          sgnz = 1.d0
         endif
      sgnxc = dsign(1.d0,xc)
      sgnyc = dsign(1.d0,yc)

      xc1 = dabs(xc)
      yc1 = dabs(yc)
      d = dmax1(dmax1(xc1,yc1), 1.d-10)

      DD = r1*d*(2.d0 - d) / dsqrt(2.d0)
      R = r1

      center = DD - dsqrt(dmax1(R**2 - DD**2, 0.d0))
      xp = DD/d * xc1
      yp = DD/d * yc1

      if (yc1 .gt. xc1) then
          yp = center + dsqrt(dmax1(R**2 - xp**2, 0.d0))
        else
          xp = center + dsqrt(dmax1(R**2 - yp**2, 0.d0))
        endif

      zp = dsqrt(dmax1(r1**2 - (xp**2 + yp**2), 0.d0))
      xp = xp*sgnxc
      yp = yp*sgnyc
      zp = zp*sgnz
c
      return
      end