MODULE Sphere  ! Basic tools for operations on a sphere
               ! and/or the surface of a spherical planet.
               ! Note that these tools are algebraic, and
               ! do not provide any graphical support.
    !
    ! by Peter Bird, UCLA, May 1997 - April 2003.
    ! copyright (c) 1997, 1998, 1999, 2003 by 
    ! Peter Bird and the Regents of the University of California.
    !-----------------------------------------------------------------
    !
    !                   CONTENTS OF THIS MODULE:
    !                  --------------------------
    ! INTENDED FOR THE USER TO CALL:
    !           REAL FUNCTION  Arc
    !              SUBROUTINE  Circles_Intersect
    !        LOGICAL FUNCTION  Clockways
    !           REAL FUNCTION  Compass
    !              SUBROUTINE  Cross
    !           REAL FUNCTION  Dot
    !           REAL FUNCTION  Easting
    !           REAL FUNCTION  Length
    !              SUBROUTINE  Local_Phi
    !              SUBROUTINE  Local_Theta
    !              SUBROUTINE  LonLat_2_ThetaPhi
    !              SUBROUTINE  LonLat_2_Uvec
    !              SUBROUTINE  Magnitude
    !              SUBROUTINE  Make_Uvec
    !              SUBROUTINE  NorthEast_Convention
    !           REAL FUNCTION  Relative_Compass
    !              SUBROUTINE  ThetaPhi_2_LonLat
    !              SUBROUTINE  ThetaPhi_2_Uvec
    !              SUBROUTINE  Turn_To
    !              SUBROUTINE  Uvec_2_LonLat
    !              SUBROUTINE  Uvec_2_ThetaPhi
    !              SUBROUTINE  Uvec_2_PlungeAzimuth
    !
    ! UTILITY ROUTINE FOR DEBUGGING:
    !              SUBROUTINE  Traceback
    !-----------------------------------------------------------------
!declares

    IMPLICIT NONE

    REAL, PARAMETER :: Pi        = 3.141592654 ! use of functions not 
    REAL, PARAMETER :: Pi_over_2 = 1.570796327 ! allowed here (too bad!)
    REAL, PARAMETER :: Two_Pi    = 6.283185307
    REAL, PARAMETER :: degrees_per_radian = 57.2957795
    REAL, PARAMETER :: radians_per_degree = 0.01745329252

    ! ---------------------------------------------------------
    ! |       General Note on Unit Vectors in a Sphere        |
    ! | The Cartesian coordinate system in which these unit   |
    ! | vectors are expressed has its origin at the center    |
    ! | of the sphere. 1st axis outcrops at ( 0 E,  0 N).     |
    ! |                2nd axis outcrops at (90 E,  0 N).     |
    ! |                3rd axis outcrops at (?? E, 90 N).     | 
    ! ---------------------------------------------------------

    CONTAINS

    REAL FUNCTION Arc (from_uvec, to_uvec)
      ! Returns length of great-circle arc in radians.
        IMPLICIT NONE
        REAL, DIMENSION(3), INTENT(IN) :: from_uvec, to_uvec
        REAL :: crossed, dotted
        REAL, DIMENSION(3) :: t_vec
        CALL Cross (from_uvec, to_uvec, t_vec)
        crossed = Length(t_vec) ! >= 0.
        dotted = Dot(from_uvec, to_uvec)
        Arc = ATAN2(crossed, dotted) ! 0. to Pi
    END FUNCTION Arc

    SUBROUTINE Circles_Intersect (pole_a_uvec, dot_a, first_a_uvec, last_a_uvec, &
                                & pole_b_uvec, dot_b, first_b_uvec, last_b_uvec, & ! input
                                & overlap, number, point1_uvec, point2_uvec)   ! output
      ! Finds all the points of intersection (0, 1, or 2) between
      !   arcs of two small circles on a unit sphere.
      ! Note: The set of small circles includes great circles
      !   as a special case.  The set of arcs includes complete-
      !   circle arcs as a special case.
      ! The two arcs are generically named "a" and "b".
      ! The pole of each is given by a "pole_x_uvec" (unit vector,
      !   from center of sphere; see module data for definitions).
      ! The plane containing the small circle is specified by "dot_x", 
      !   its distance from the origin.  (If either distance is 
      !   greater than 1.00, there can be no points of intersection.)
      !   Parameters "dot_a" and "dot_b" are DOUBLE PRECISION in order
      !   to provide sufficient precision for the definition of very
      !   small circles (e.g., earthquake epicenters).
      !"First" and "last" points on each arc are specified by uvecs
      !   (unit vectors) "first_x_uvec" and "last_x_uvec", 
      !   where one goes counterclockwise around the pole
      !   first to last (assuming a viewpoint outside the
      !   unit sphere).  If "first" = "last", the arc is a full circle.
      ! NOTE: Although "first_x_uvec" and "last_x_uvec" will usually
      !   be points on the corresponding arcs, they need not be.
      !   Only their azimuths from "pole_x_uvec" are used.  In fact,
      !   if "first_x_uvec == last_x_uvec", then even the azimuths are
      !   not used, so ANY vector (even a zero vector) can be sent in
      !   these two positions to signal a complete small circle.
      ! CAUTION: Results will be strange and unpredictable if the
      !   the "unit vectors" supplied are not actually 1.00 long!
      ! Values returned are "number" (the count of intersection points)
      !   and unit vectors for as many points as were found.
      ! If only one intersection is found, it is always placed in point1_uvec.
      ! In the special case where the two small circles are the same
      !   circle, and share some common arc, the logical flag "overlap"
      !   is set, and the first and last points of the common arc are
      !   reported, and number = 2 on output.  If both of these overlapped
      !   circles are complete circles, these first and last points are
      !   the same point.
        IMPLICIT NONE
        REAL, DIMENSION(3), INTENT(IN)  :: pole_a_uvec, first_a_uvec, last_a_uvec, &
                                        &  pole_b_uvec, first_b_uvec, last_b_uvec
        DOUBLE PRECISION, INTENT(IN)    :: dot_a, dot_b
        LOGICAL, INTENT(OUT)            :: overlap
        INTEGER, INTENT(OUT)            :: number
        REAL, DIMENSION(3), INTENT(OUT) :: point1_uvec, point2_uvec

        DOUBLE PRECISION :: a_dot_b, alpha, beta, determinant, t
        DOUBLE PRECISION, DIMENSION(2,2) :: inverse
        LOGICAL :: ring_a, ring_b ! are the arcs actually complete circles?
        REAL :: along, antipole_gap, &
              & first_a_radians, first_b_radians, last_a_radians, &
              & last_b_radians, min_radius, point_wrt_a_radians, point_wrt_b_radians, &
              & pole_gap
        REAL, PARAMETER :: tolerance = 0.0014 ! NO SMALLER; otherwise,
                                              ! Determinant = 0.0 error may appear.
        REAL, DIMENSION(3) :: line_uvec, offset, t_vec

      ! Check for an easy answer (no intersections):
        number = 0
        overlap = .FALSE.
        IF (dot_a >=  1.D0) RETURN
        IF (dot_a <= -1.D0) RETURN
        IF (dot_b >=  1.D0) RETURN
        IF (dot_b <= -1.D0) RETURN
 
      ! Decide whether arcs are complete circles; otherwise, determine
      !   (relative) azimuths of endpoints:
        ring_a = (first_a_uvec(1) == last_a_uvec(1)).AND. &
               & (first_a_uvec(2) == last_a_uvec(2)).AND. &
               & (first_a_uvec(3) == last_a_uvec(3))
        ring_b = (first_b_uvec(1) == last_b_uvec(1)).AND. &
               & (first_b_uvec(2) == last_b_uvec(2)).AND. &
               & (first_b_uvec(3) == last_b_uvec(3))
        IF (.NOT.ring_a) THEN
            first_a_radians = -Relative_Compass(pole_a_uvec, first_a_uvec)
            last_a_radians  = -Relative_Compass(pole_a_uvec, last_a_uvec)
            IF ((first_a_radians - last_a_radians) > 1.E-6) last_a_radians = last_a_radians + Two_Pi
           !Note: The reason for the seemingly gratuitous tolerance is that 
           !      under Digital Visual Fortran 5.0D, a value of 
           !      0.3510835 was treated as "less than" a value of 0.3510835,
           !      turning a VERY short arc into a nearly-complete small circle,
           !      generating erroneous intersections, and putting unwanted
           !      line-segments across the plot to connect to the map boundary!
        END IF
        IF (.NOT.ring_b) THEN
            first_b_radians = -Relative_Compass(pole_b_uvec, first_b_uvec)
            last_b_radians  = -Relative_Compass(pole_b_uvec, last_b_uvec)
            IF ((first_b_radians - last_b_radians) > 1.E-6) last_b_radians = last_b_radians + Two_Pi
           !See note above.
        END IF

      ! Test for special cases of parallel and antipodal poles:
        t_vec = pole_a_uvec - pole_b_uvec
        pole_gap = Length (t_vec)
        t_vec = pole_a_uvec + pole_b_uvec
        antipole_gap = Length (t_vec)
        IF (pole_gap <= tolerance) THEN ! poles virtually identical
            IF (dot_a /= dot_b) RETURN ! no intersection
            ! From here on, assume dot_a == dot_b:
            IF (ring_a) THEN
                overlap = .TRUE.
                number = 2
                point1_uvec = first_b_uvec
                point2_uvec = last_b_uvec
            ELSE IF (ring_b) THEN
                overlap = .TRUE.
                number = 2
                point1_uvec = first_a_uvec
                point2_uvec = last_a_uvec
            ELSE ! neither circle is complete
                IF ((first_b_radians >= first_a_radians).AND. &
                   &(first_b_radians <= last_a_radians)) THEN ! 1st of b is in a.
                    point1_uvec = first_b_uvec
                    IF (last_b_radians <= last_a_radians) THEN
                        point2_uvec = last_b_uvec
                    ELSE
                        point2_uvec = last_a_uvec
                    END IF
                    IF ((point1_uvec(1) == point2_uvec(1)).AND. &
                       &(point1_uvec(2) == point2_uvec(2)).AND. &
                       &(point1_uvec(3) == point2_uvec(3))) THEN ! only one common point
                        overlap = .FALSE.
                        number = 1
                    ELSE ! a common arc
                        overlap = .TRUE.
                        number = 2 
                    END IF
                ELSE IF ((first_a_radians >= first_b_radians).AND. &
                        &(first_a_radians <= last_b_radians)) THEN ! 1st of a is in b.
                    point1_uvec = first_a_uvec
                    IF (last_a_radians <= last_b_radians) THEN
                        point2_uvec = last_a_uvec
                    ELSE
                        point2_uvec = last_b_uvec
                    END IF
                    IF ((point1_uvec(1) == point2_uvec(1)).AND. &
                       &(point1_uvec(2) == point2_uvec(2)).AND. &
                       &(point1_uvec(3) == point2_uvec(3))) THEN ! only one common point
                        overlap = .FALSE.
                        number = 1
                    ELSE ! a common arc
                        overlap = .TRUE.
                        number = 2 
                    END IF
                END IF ! separate arcs of same circle [no ELSE; we RETURN, with number = 0]
            END IF ! either circle is complete, or neither is
        ELSE IF (antipole_gap <= tolerance) THEN ! antipodal
            IF (dot_a /= -dot_b) RETURN ! no intersection
            ! From here on, assume dot_a == -dot_b:
            IF (ring_a) THEN
                overlap = .TRUE.
                number = 2
                point1_uvec = first_b_uvec
                point2_uvec = last_b_uvec
            ELSE IF (ring_b) THEN
                overlap = .TRUE.
                number = 2
                point1_uvec = first_a_uvec
                point2_uvec = last_a_uvec
            ELSE ! neither circle is complete
                ! Because of antipodal relationship, azimuths are confusing.
                ! Redefine azimuths of arc b in terms of pole a,
                !  while reversing the direction along arc b to make it
                !  counterclockwise about a:
                first_b_radians = -Relative_Compass(pole_a_uvec, last_b_uvec)
                last_b_radians  = -Relative_Compass(pole_a_uvec, first_b_uvec)
                IF (last_b_radians < first_b_radians) last_b_radians = last_b_radians + Two_Pi
                IF ((first_b_radians >= first_a_radians).AND. &
                   &(first_b_radians <= last_a_radians)) THEN ! 1st of b is in a.
                    point1_uvec = first_b_uvec
                    IF (last_b_radians <= last_a_radians) THEN
                        point2_uvec = last_b_uvec
                    ELSE
                        point2_uvec = last_a_uvec
                    END IF
                    IF ((point1_uvec(1) == point2_uvec(1)).AND. &
                       &(point1_uvec(2) == point2_uvec(2)).AND. &
                       &(point1_uvec(3) == point2_uvec(3))) THEN ! only one common point
                        overlap = .FALSE.
                        number = 1
                    ELSE ! a common arc
                        overlap = .TRUE.
                        number = 2 
                    END IF
                ELSE IF ((first_a_radians >= first_b_radians).AND. &
                        &(first_a_radians <= last_b_radians)) THEN ! 1st of a is in b.
                    point1_uvec = first_a_uvec
                    IF (last_a_radians <= last_b_radians) THEN
                        point2_uvec = last_a_uvec
                    ELSE
                        point2_uvec = last_b_uvec
                    END IF
                    IF ((point1_uvec(1) == point2_uvec(1)).AND. &
                       &(point1_uvec(2) == point2_uvec(2)).AND. &
                       &(point1_uvec(3) == point2_uvec(3))) THEN ! only one common point
                        overlap = .FALSE.
                        number = 1
                    ELSE ! a common arc
                        overlap = .TRUE.
                        number = 2 
                    END IF
                END IF ! separate arcs on same circle [no ELSE; we RETURN, with number = 0]
            END IF ! either circle is complete, or neither is

        ELSE ! **** normal case; unrelated poles *****
            ! Each small circle lies in its own plane.
            ! These two planes intersect in a line in 3-D.
            ! Find "offset" = (non-unit) vector to point on this line
            !   which is closest to origin:
            a_dot_b = (1.D0 * pole_a_uvec(1)) * (1.D0 * pole_b_uvec(1)) + &
                    & (1.D0 * pole_a_uvec(2)) * (1.D0 * pole_b_uvec(2)) + &
                    & (1.D0 * pole_a_uvec(3)) * (1.D0 * pole_b_uvec(3))
            !Note that this is the double-precision equivalent of:
            !   a_dot_b = Dot (pole_a_uvec, pole_b_uvec)
            determinant = 1.D0 - a_dot_b**2
            IF (determinant == 0.0D0) THEN
                WRITE (*,"(' ERROR: Determinant = 0.0D0')")
                CALL Traceback
            END IF
            t = 1.0D0 / determinant
            inverse(1,1) =  t           !  t * matrix(2,2) 
            inverse(1,2) = -t * a_dot_b ! -t * matrix(1,2)
            inverse(2,1) = -t * a_dot_b ! -t * matrix(2,1) 
            inverse(2,2) =  t           !  t * matrix(1,1) 
            alpha = inverse(1,1)*dot_a + inverse(1,2)*dot_b
            beta  = inverse(2,1)*dot_a + inverse(2,2)*dot_b
            offset(1) = alpha*pole_a_uvec(1) + beta*pole_b_uvec(1)
            offset(2) = alpha*pole_a_uvec(2) + beta*pole_b_uvec(2)
            offset(3) = alpha*pole_a_uvec(3) + beta*pole_b_uvec(3)
           !Note: "offset" is REAL(3), and this the end of DOUBLE PRECISON work
            min_radius = Length (offset)
            IF (min_radius == 1.0) THEN ! circles osculate at one point
                ! Check longitudes about poles to see if point is in arcs:
                IF (.NOT.ring_a) THEN
                    point_wrt_a_radians = -Relative_Compass(pole_a_uvec, offset)
                    IF (point_wrt_a_radians < first_a_radians) point_wrt_a_radians = point_wrt_a_radians + Two_Pi
                END IF
                IF (ring_a .OR. &
                   &((point_wrt_a_radians >= first_a_radians).AND. &
                   & (point_wrt_a_radians <= last_a_radians))) THEN
                    IF (.NOT.ring_b) THEN
                        point_wrt_b_radians = -Relative_Compass(pole_b_uvec, offset)
                        IF (point_wrt_b_radians < first_b_radians) point_wrt_b_radians = point_wrt_b_radians + Two_Pi
                    END IF
                    IF (ring_b .OR. &
                       &((point_wrt_b_radians >= first_b_radians).AND. &
                       & (point_wrt_b_radians <= last_b_radians))) THEN
                        number = 1
                        point1_uvec = offset
                    END IF ! in arc b [no ELSE; we RETURN with number = 0]
                END IF ! in arc a [no ELSE; we RETURN with number = 0]
            ELSE IF (min_radius < 1.) THEN ! two intersection points between circles
                ! Find vector parallel to line of intersection of the two
                !    planes which contain the two small circles:
                CALL Cross (pole_a_uvec, pole_b_uvec, t_vec)
                CALL Make_Uvec (t_vec, line_uvec)
                ! Now, points on this line are expressed by
                !   vector "offset" plus any multiple of vector "line_uvec".
                ! Call the multiple "along".
                ! Solve the hyperbolic equation that says that the
                !   radius of a point on this line is 1.00
                !  (i.e., it is also a point on the unit sphere):
                !         min_radius**2 + along**2 = 1.00**2
                along = SQRT ( 1. - min_radius**2 )
                point1_uvec = offset + along * line_uvec
                point2_uvec = offset - along * line_uvec
                ! Remember, these are only tentative solutions.

                ! Check longitudes about poles to see if point1 is in arcs:
                IF (.NOT.ring_a) THEN
                    point_wrt_a_radians = -Relative_Compass(pole_a_uvec, point1_uvec)
                    IF (point_wrt_a_radians < first_a_radians) point_wrt_a_radians = point_wrt_a_radians + Two_Pi
                END IF
                IF (ring_a .OR. &
                   &((point_wrt_a_radians >= first_a_radians).AND. &
                   & (point_wrt_a_radians <= last_a_radians))) THEN
                    IF (.NOT.ring_b) THEN
                        point_wrt_b_radians = -Relative_Compass(pole_b_uvec, point1_uvec)
                        IF (point_wrt_b_radians < first_b_radians) point_wrt_b_radians = point_wrt_b_radians + Two_Pi
                    END IF
                    IF (ring_b .OR. &
                       &((point_wrt_b_radians >= first_b_radians).AND. &
                       & (point_wrt_b_radians <= last_b_radians))) THEN
                        number = 1
                    END IF ! point1 also in arc b: a solution.
                END IF ! point1 is in arc a

                ! Check longitudes about poles to see if point2 is in arcs:
                IF (.NOT.ring_a) THEN
                    point_wrt_a_radians = -Relative_Compass(pole_a_uvec, point2_uvec)
                    IF (point_wrt_a_radians < first_a_radians) point_wrt_a_radians = point_wrt_a_radians + Two_Pi
                END IF
                IF (ring_a .OR. &
                   &((point_wrt_a_radians >= first_a_radians).AND. &
                   & (point_wrt_a_radians <= last_a_radians))) THEN
                    IF (.NOT.ring_b) THEN
                        point_wrt_b_radians = -Relative_Compass(pole_b_uvec, point2_uvec)
                        IF (point_wrt_b_radians < first_b_radians) point_wrt_b_radians = point_wrt_b_radians + Two_Pi
                    END IF
                    IF (ring_b .OR. &
                       &((point_wrt_b_radians >= first_b_radians).AND. &
                       & (point_wrt_b_radians <= last_b_radians))) THEN
                        IF (number == 1) THEN ! add a second solution
                            number = 2
                        ELSE ! point2 is the first and only solution
                            number = 1
                            point1_uvec = point2_uvec
                        END IF ! is this the first or the second solution?
                    END IF ! point2 in arc b: a solution.
                END IF ! point2 is in arc a

            END IF ! min_radius = 1., or < 1. [no ELSE; we RETURN with number = 0]
        END IF ! poles parallel, or antipodal, or unrelated
    END SUBROUTINE Circles_Intersect


    REAL FUNCTION Compass (from_uvec, to_uvec)
      ! Returns azimuth (in radians, clockwise from North)
      !   of the great-circle arc from "from_uvec" to "to_uvec",
      !   measured at location "from_uvec".
      ! Does NOT work at North or South pole!  (See Relative_Compass.)
        IMPLICIT NONE
        REAL, DIMENSION(3), INTENT(IN) :: from_uvec, to_uvec
        REAL :: ve, vn
        REAL, DIMENSION(3) :: omega_vec, omega_uvec, &
                            & Phi, Theta, &
                            & v_vec
        IF ((from_uvec(1) == 0.).AND.(from_uvec(2) == 0.)) THEN
            WRITE (*,"(' ERROR: Compass undefined at N or S pole.  Use Relative_Compass.')")
            CALL Traceback
        ELSE IF ((from_uvec(1) == to_uvec(1)).AND. &
           &(from_uvec(2) == to_uvec(2)).AND. &
           &(from_uvec(3) == to_uvec(3))) THEN
            WRITE (*,"(' ERROR: Compass bearing from point to itself undefined.')")
            CALL Traceback
        ELSE IF ((from_uvec(1) == -to_uvec(1)).AND. &
                &(from_uvec(2) == -to_uvec(2)).AND. &
                &(from_uvec(3) == -to_uvec(3))) THEN
            WRITE (*,"(' ERROR: Compass bearing from point to antipode undefined.')")
            CALL Traceback
        ELSE
          ! Normal case:
            CALL Cross (from_uvec, to_uvec, omega_vec)
            CALL Make_Uvec(omega_vec, omega_uvec)
            CALL Cross (omega_uvec, from_uvec, v_vec)
            CALL Local_Theta(from_uvec, Theta)
            CALL Local_Phi  (from_uvec, Phi)
            vn = -Dot(v_vec, Theta)
            ve =  Dot(v_vec, Phi)
           Compass = ATAN2(ve, vn)
        END IF
    END FUNCTION Compass

    SUBROUTINE Cross (a_vec, b_vec, c_vec)
      ! vector cross product: a x b = c
        IMPLICIT NONE
        REAL, DIMENSION(3), INTENT(IN)  :: a_vec, b_vec
        REAL, DIMENSION(3), INTENT(OUT) :: c_vec
        c_vec(1) = a_vec(2)*b_vec(3) - a_vec(3)*b_vec(2)
        c_vec(2) = a_vec(3)*b_vec(1) - a_vec(1)*b_vec(3)
        c_vec(3) = a_vec(1)*b_vec(2) - a_vec(2)*b_vec(1)
    END SUBROUTINE Cross

    REAL FUNCTION Dot (a_vec, b_vec)
      ! returns scalar (dot) product of two 3-component vectors
        IMPLICIT NONE
        REAL, DIMENSION(3), INTENT(IN) :: a_vec, b_vec
        Dot = a_vec(1)*b_vec(1) + a_vec(2)*b_vec(2) + a_vec(3)*b_vec(3)
    END FUNCTION Dot

    REAL FUNCTION Easting(delta_lon_degrees)
        !returns positive result, 0.-359.999
        IMPLICIT NONE
        REAL, INTENT(IN) :: delta_lon_degrees
        Easting = MOD((delta_lon_degrees + 720.),360.)
    END FUNCTION Easting

    INTEGER FUNCTION IAbove(x)
      ! returns first integer >= x, unlike INT() or NINT():
        IMPLICIT NONE
        REAL, INTENT(IN) :: x
        INTEGER answer
        REAL :: t
        answer = INT(x)
        IF (x > 0.) THEN
            t = 1.00 * answer
            IF (x > t) THEN
                answer = answer + 1
            END IF
        END IF
        IAbove = answer
    END FUNCTION IAbove

    REAL FUNCTION Length(a_vec)
        IMPLICIT NONE
        REAL, DIMENSION(3), INTENT(IN) :: a_vec
        DOUBLE PRECISION :: t
        t = a_vec(1)*1.D0*a_vec(1) + &
          & a_vec(2)*1.D0*a_vec(2) + &
          & a_vec(3)*1.D0*a_vec(3)
        IF (t == 0.D0) THEN
            Length = 0.
        ELSE
            Length = SQRT(t)
        END IF
    END FUNCTION Length

    SUBROUTINE Local_Phi (b_, Phi)
        ! returns local East-pointing unit vector in Cartesian coordinates
        ! for location b_; not intended to work at the poles!
        REAL, DIMENSION(3), INTENT(IN)  :: b_
        REAL, DIMENSION(3), INTENT(OUT) :: Phi
        REAL, DIMENSION(3)              :: temp
        IF (b_(1) == 0.) THEN
            IF (b_(2) == 0.) THEN
                WRITE (*,"(' ERROR: Local_Phi was requested for N or S pole.')")
                CALL Traceback
             END IF
        END IF
        temp(1) = - b_(2)
        temp(2) = b_(1)
        temp(3) = 0.
        CALL Make_Uvec(temp, Phi)
    END SUBROUTINE Local_Phi

    SUBROUTINE Local_Theta (b_, Theta)
        ! returns local South-pointing unit vector in Cartesian coordinates
        ! for location b_; not intended to work at the poles!
        REAL, DIMENSION(3), INTENT(IN)  :: b_
        REAL, DIMENSION(3), INTENT(OUT) :: Theta
        REAL, DIMENSION(3) :: temp
        REAL  :: equat, new_equat
        equat = SQRT(b_(1)**2 + b_(2)**2)   !equatorial component
        IF (equat == 0.) THEN
            WRITE (*,"(' ERROR: Local_Theta was requested for N or S pole.')")
            CALL Traceback
        END IF
        new_equat = b_(3)   ! swap components in a meridional plane
        temp(3) = - equat       ! "
        temp(1) = new_equat * b_(1) / equat ! partition new equatorial component
        temp(2) = new_equat * b_(2) / equat ! "
        CALL Make_Uvec (temp, Theta)
    END SUBROUTINE Local_Theta

    SUBROUTINE LonLat_2_ThetaPhi (lon, lat, theta, phi)
      ! "lon" is East longitude in degrees.
      ! "lat" is North latitude in degrees.
      ! "theta" is co-latitude, from N pole, in radians
      ! "phi" is East longitude in radians
        IMPLICIT NONE
        REAL, INTENT(IN)  :: lon, lat
        REAL, INTENT(OUT) :: theta, phi
        IF ((lat > 90.).OR.(lat < -90.)) THEN
            WRITE (*,"(' ERROR: Latitude outside legal range.')")
            CALL Traceback
        ELSE
            theta = radians_per_degree * (90. - lat)
            phi = radians_per_degree * lon
        END IF
    END SUBROUTINE LonLat_2_ThetaPhi

    SUBROUTINE LonLat_2_Uvec (lon, lat, uvec)
        IMPLICIT NONE
        REAL, INTENT(IN) :: lon, lat
        REAL, DIMENSION(3), INTENT(OUT) :: uvec
        REAL :: theta, phi
        CALL LonLat_2_ThetaPhi (lon, lat, theta, phi)
        CALL ThetaPhi_2_Uvec (theta, phi, uvec)
    END SUBROUTINE LonLat_2_Uvec

    REAL FUNCTION Magnitude (b_)
        REAL, DIMENSION(3), INTENT(IN) :: b_
        Magnitude = SQRT(b_(1)**2 +b_(2)**2 + b_(3)**2)
    END FUNCTION Magnitude

    SUBROUTINE MATHERRQQ( name, length, info, retcode)
        !Provided so that domain errors, underflows, etc.
        !  can be trapped and debugged; otherwise program crashes!
        USE DFLIB
        INTEGER(2) length, retcode
        CHARACTER(length) name
        RECORD /MTH$E_INFO/ info
        PRINT *, "Entered MATHERRQQ"
        PRINT *, "Failing function is: ", name
        PRINT *, "Error type is: ", info.errcode
        IF ((info.ftype == TY$REAL4 ).OR.(info.ftype == TY$REAL8)) THEN
            PRINT *, "Type: REAL"
            PRINT *, "Enter the desired function result: "
            READ(*,*) info.r8res
            retcode = 1
        ELSE IF ((info.ftype == TY$CMPLX8 ).OR.(info.ftype == TY$CMPLX16)) THEN
            PRINT *, "Type: COMPLEX"
            PRINT *, "Enter the desired function result: "
            READ(*,*) info.c16res
            retcode = 1
        END IF
    END SUBROUTINE MATHERRQQ

    SUBROUTINE Make_Uvec (vector, uvec)
      ! Shortens or lengthens a three-component vector to a unit vector;
      !    includes special kludge to prevent extremely small component
      !    values which result from rounding error and result in later
      !    numerical underflows.
        IMPLICIT NONE
        INTEGER :: i
        REAL, DIMENSION(3), INTENT(IN) :: vector
        REAL, DIMENSION(3), INTENT(OUT):: uvec
        REAL :: factor, size
        size = Length(vector)
        IF (size > 0.) THEN
           factor = 1. / size
           uvec = vector * factor
           DO i = 1, 3
               IF (ABS(uvec(i)) < 1.E-6) uvec(i) = 0.0
           END DO
        ELSE
           WRITE (*,"(' ERROR: Cannot Make_Uvec of (0., 0., 0.).')")
           CALL Traceback
        END IF
    END SUBROUTINE Make_Uvec

    SUBROUTINE NorthEast_Convention (location_uvec, north_uvec, east_uvec)
      ! At most positions ("location_uvec") returns "north_uvec"
      ! (same as -Local_Theta) and "east_uvec" (same as +Local_Phi).
      ! However, within a small distance from either the North or South
      ! pole, it adopts arbitrary conventional directions based on 
      ! the limiting directions along the 0E meridian as the latitude
      ! approaches +90N or -90N.  Since both FUNCTION  Relative_Compass
      ! and SUBROUTINE  Turn_To call this routine, they can work 
      ! together even when location_uvec is at, or near, a pole.
        IMPLICIT NONE
        REAL, DIMENSION(3), INTENT(IN)  :: location_uvec
        REAL, DIMENSION(3), INTENT(OUT) :: north_uvec, east_uvec
        REAL, PARAMETER :: pole_area = 7.615E-7 ! (0.05 degrees -> radians -> squared)
        REAL :: equat2
        REAL, DIMENSION(3) :: t_vec
        equat2 = location_uvec(1)**2 + location_uvec(2)**2
        IF (equat2 > pole_area) THEN
          ! normal case:
            t_vec(1) = -location_uvec(2)
            t_vec(2) = +location_uvec(1)
            t_vec(3) = 0.0
            CALL Make_Uvec (t_vec, east_uvec)
        ELSE
          ! very close to N or S pole; act as if on 0E meridian:
            east_uvec  = (/ 0., 1., 0. /)
        END IF
        CALL Cross (location_uvec, east_uvec, north_uvec)
    END SUBROUTINE NorthEast_Convention

    REAL FUNCTION Relative_Compass (from_uvec, to_uvec)
      ! At most points, works exactly like Compass:
      !   returns azimuth (in radians, clockwise from North)
      !   of the great-circle arc from "from_uvec" to "to_uvec",
      !   measured at location from_uvec.
      ! However, unlike Compass, it does not crash at the N and S poles,
      !   where Theta and Phi are undefined.  Instead, it uses
      !   SUBROUTINE  NorthEast_Convention to make an
      !   arbitrary choice of axes, so that RELATIVE azimuths can
      !   be measured from pivot points at the poles, by multiple calls.
        IMPLICIT NONE
        REAL, DIMENSION(3), INTENT(IN) :: from_uvec, to_uvec
        REAL :: ve, vn
        REAL, DIMENSION(3) :: omega_vec, omega_uvec, &
                            & north_uvec, east_uvec, &
                            & v_vec
        IF ((from_uvec(1) == to_uvec(1)).AND. &
           &(from_uvec(2) == to_uvec(2)).AND. &
           &(from_uvec(3) == to_uvec(3))) THEN
            WRITE (*,"(' ERROR: Compass bearing from point to itself undefined.')")
            CALL Traceback
        ELSE IF ((from_uvec(1) == -to_uvec(1)).AND. &
                &(from_uvec(2) == -to_uvec(2)).AND. &
                &(from_uvec(3) == -to_uvec(3))) THEN
            WRITE (*,"(' ERROR: Compass bearing from point to antipode undefined.')")
            CALL Traceback
        ELSE
          ! Normal case:
            CALL Cross (from_uvec, to_uvec, omega_vec)
            CALL Make_Uvec(omega_vec, omega_uvec)
            CALL Cross (omega_uvec, from_uvec, v_vec)
            CALL NorthEast_Convention (from_uvec, north_uvec, east_uvec)
            vn = Dot(v_vec, north_uvec)
            ve = Dot(v_vec, east_uvec)
            Relative_Compass = ATAN2(ve, vn)
        END IF
    END FUNCTION Relative_Compass

    SUBROUTINE ThetaPhi_2_LonLat (theta, phi, lon, lat)
      ! Converts from theta (co-latitude, from N pole, in radians)
      ! and phi (longitude, E from Greenwich, in radians)
      ! to "lon" (East longitude, in degrees; West is negative)
      ! and "lat" (North latitude, in degrees; South is negative).
        IMPLICIT NONE
        REAL, INTENT(IN)  :: theta, phi
        REAL, INTENT(OUT) :: lon, lat
        lat = 90. - degrees_per_radian * ABS(theta)
        lat = MAX (lat, -90.)
        lon = degrees_per_radian * phi
        IF (lon > 180.) lon = lon - 360.
        IF (lon <= -180.) lon = lon + 360.
    END SUBROUTINE ThetaPhi_2_LonLat

    SUBROUTINE ThetaPhi_2_Uvec (theta, phi, uvec)
      ! Converts from theta (co-latitude, from N pole) and
      ! phi (longitude, E from Greenwich) [both in radians]
      ! to a 3-component Cartesian unit vector, which points
      ! from the center of the unit sphere to a surface point.
      ! Its first axis outcrops at (0E, 0N).
      ! Its second axis outcrops at (90E, 0N).
      ! Its third axis outcrops at 90N.
        IMPLICIT NONE
        REAL, INTENT(IN) :: theta, phi
        REAL, DIMENSION(3), INTENT(OUT) :: uvec
        REAL :: equat
        uvec(3) = COS(theta)
        equat = SIN(theta)
        uvec(1) = equat * COS(phi)
        uvec(2) = equat * SIN(phi)
    END SUBROUTINE ThetaPhi_2_Uvec

    SUBROUTINE Traceback ()
        ! The sole function of this unit is to cause a traceable error,
        !   so that the route into the unit that called it is also traced.
        ! This unit is a good place to put a breakpoint while debugging!
        ! The intentional error must NOT be detected during compilation,
        !   but MUST cause a traceable error at run-time.
        ! If this routine has any error detected during compilation,
        !   then change its code to cause a different intentional error.
        IMPLICIT NONE
        CHARACTER*80 instring
        INTEGER :: i
        REAL, DIMENSION(3) :: y
        WRITE (*,"(' -----------------------------------------------------')")
        WRITE (*,"(' Traceback was called to execute an intentional error:')")
        WRITE (*,"(' An array subscript will be intentionally out-of-range.')")
        WRITE (*,"(/' After you read this notice, press [Enter]' &
                 & /' to stop the program (no other option): '\)")
        READ (*,"(A)") instring
        DO i = 1, 4
            y(i) = 1. * i
        END DO
        STOP ' '
    END SUBROUTINE Traceback

    SUBROUTINE Turn_To (azimuth_radians, base_uvec, far_radians, & ! inputs
                     &  omega_uvec, result_uvec)
      ! Computes uvec "result_uvec" (a 3-component Cartesian unit
      ! vector from the center of the planet) which results from 
      ! rotating along a great circle beginning at "base_uvec"
      ! for an angle of "far_radians", in the initial direction 
      ! given by azimuth_radians" (clockwise, from North).
      ! Also returned is "omega_uvec", the pole of rotation.
      ! NOTE: At the poles, azimuth is undefined.  Near the
      ! poles, it is defined but numerically unstable.  Therefore,
      ! Turn_To uses the same SUBROUTINE  NorthEast_Convention as
      ! FUNCTION  Relative_Compass does, so they can work together
      ! to find internal points on a small circle (as in Process_L4_
      ! Paths).  
        IMPLICIT NONE
        REAL,                 INTENT(IN)  :: azimuth_radians, far_radians
        REAL, DIMENSION(3),   INTENT(IN)  :: base_uvec
        REAL, DIMENSION(3),   INTENT(OUT) :: omega_uvec, result_uvec
        REAL, PARAMETER :: pole_width = 1.745E-4 ! 0.01 degrees; must match Turn_To!
        REAL :: complement, cos_size, e_part, n_part, sin_size
        REAL, DIMENSION(3)   :: east_uvec, north_uvec, t_uvec
        REAL, DIMENSION(3,3) :: rotation_matrix
        CALL NorthEast_Convention (base_uvec, north_uvec, east_uvec)
        e_part = -COS(azimuth_radians)
        n_part =  SIN(azimuth_radians)
        omega_uvec(1) = e_part*east_uvec(1) + n_part*north_uvec(1)
        omega_uvec(2) = e_part*east_uvec(2) + n_part*north_uvec(2)
        omega_uvec(3) = e_part*east_uvec(3) + n_part*north_uvec(3)
        cos_size = COS(far_radians)
        sin_size = SIN(far_radians)
        complement = 1.00 - cos_size
        rotation_matrix(1,1) = cos_size + complement*omega_uvec(1)*omega_uvec(1)
        rotation_matrix(1,2) = complement*omega_uvec(1)*omega_uvec(2) - sin_size*omega_uvec(3)
        rotation_matrix(1,3) = complement*omega_uvec(1)*omega_uvec(3) + sin_size*omega_uvec(2)
        rotation_matrix(2,1) = complement*omega_uvec(2)*omega_uvec(1) + sin_size*omega_uvec(3)
        rotation_matrix(2,2) = cos_size + complement*omega_uvec(2)*omega_uvec(2)
        rotation_matrix(2,3) = complement*omega_uvec(2)*omega_uvec(3) - sin_size*omega_uvec(1)
        rotation_matrix(3,1) = complement*omega_uvec(3)*omega_uvec(1) - sin_size*omega_uvec(2)
        rotation_matrix(3,2) = complement*omega_uvec(3)*omega_uvec(2) + sin_size*omega_uvec(1)
        rotation_matrix(3,3) = cos_size + complement*omega_uvec(3)*omega_uvec(3)
       !Copy base_uvec in case user of this routine plans to change it:
        t_uvec = base_uvec
        result_uvec(1) = rotation_matrix(1,1)*t_uvec(1) + &
                       & rotation_matrix(1,2)*t_uvec(2) + &
                       & rotation_matrix(1,3)*t_uvec(3)
        result_uvec(2) = rotation_matrix(2,1)*t_uvec(1) + &
                       & rotation_matrix(2,2)*t_uvec(2) + &
                       & rotation_matrix(2,3)*t_uvec(3)
        result_uvec(3) = rotation_matrix(3,1)*t_uvec(1) + &
                       & rotation_matrix(3,2)*t_uvec(2) + &
                       & rotation_matrix(3,3)*t_uvec(3)
    END SUBROUTINE Turn_To

    SUBROUTINE Uvec_2_LonLat (uvec, lon, lat)
        IMPLICIT NONE
        REAL, DIMENSION(3), INTENT(IN) :: uvec
        REAL, INTENT(OUT) :: lon, lat
        REAL :: theta, phi
        CALL Uvec_2_ThetaPhi (uvec, theta, phi)
        CALL ThetaPhi_2_LonLat (theta, phi, lon, lat)
    END SUBROUTINE Uvec_2_LonLat 

    SUBROUTINE Uvec_2_ThetaPhi (uvec, theta, phi)
      ! converts from Cartesian unit vector to theta (colatitude)
      ! and phi (longitude), both in radians
        IMPLICIT NONE
        REAL, DIMENSION(3), INTENT(IN) :: uvec
        REAL, INTENT(OUT) :: theta, phi
        REAL :: equat, equat2
        equat2 = uvec(1)*uvec(1) + uvec(2)*uvec(2)
        IF (equat2 == 0.) THEN
            phi = 0. ! actually undefined; default 0.
            IF (uvec(3) > 0.) THEN
                theta = 0. ! N pole
            ELSE 
                theta = Pi ! S pole
            END IF
        ELSE
            equat = SQRT(equat2)
            theta = ATAN2(equat, uvec(3))
            phi = ATAN2(uvec(2), uvec(1))
        END IF
    END SUBROUTINE Uvec_2_ThetaPhi

    SUBROUTINE Uvec_2_PlungeAzimuth (location_uvec, lineation_uvec, & ! inputs
                                  &  plunge_degrees, azimuth_degrees) ! outputs
      ! converts Cartesian unit vector lineation_uvec
      ! which is at geographic location location_uvec
      ! into plunge toward azimuth, in degrees.
        IMPLICIT NONE
        REAL, DIMENSION(3), INTENT(IN) :: location_uvec, lineation_uvec
        INTEGER, INTENT(OUT) :: plunge_degrees, azimuth_degrees
        REAL :: equat, horizontal_component, phi, phi_component, theta, theta_component, up_component
        REAL, DIMENSION(3) :: phi_uvec, theta_uvec
        equat = SQRT(location_uvec(1)*location_uvec(1) + location_uvec(2)*location_uvec(2))
        IF (equat == 0.) THEN
            WRITE (*,"(' ERROR: location_uvec sent to Uvec_2_PlungeAzimuth may not be N or S pole.')")
            CALL Traceback
        END IF
        CALL Local_Theta(location_uvec, theta_uvec)
        CALL Local_Phi  (location_uvec, phi_uvec)
        up_component    = Dot(lineation_uvec, location_uvec)
        theta_component = Dot(lineation_uvec, theta_uvec)
        phi_component   = Dot(lineation_uvec, phi_uvec)
        horizontal_component = SQRT(theta_component**2 + phi_component**2)
        IF (horizontal_component <= 0.0) THEN
            plunge_degrees = 90
            azimuth_degrees = 0
        ELSE IF (up_component <= 0.0) THEN
            plunge_degrees = NINT(degrees_per_radian * ATAN2(-up_component, horizontal_component))
            azimuth_degrees = NINT(degrees_per_radian * ATAN2(phi_component, -theta_component))
            IF (azimuth_degrees < 0) azimuth_degrees = azimuth_degrees + 360
        ELSE ! up_component is postive
            plunge_degrees = NINT(degrees_per_radian * ATAN2(up_component, horizontal_component))
            azimuth_degrees = NINT(degrees_per_radian * ATAN2(-phi_component, theta_component))
            IF (azimuth_degrees < 0) azimuth_degrees = azimuth_degrees + 360
        END IF
    END SUBROUTINE Uvec_2_PlungeAzimuth

END MODULE Sphere !===============================================
