MODULE DSphere ! 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.
! Exactly parallel to MODULE Sphere.f90,
! except that all REAL variables have been
! replaced with DOUBLE PRECISION variables,
! and subprogram names have "D" added as a
! sign that they expect DOUBLE PRECISION
! arguments and return DOUBLE PRECISION answers.
!
! By Peter Bird, UCLA, May 1997 - April 2003, & December 2014;
! update to REAL*8 (DOUBLE PRECISION) in July 2015.
! Copyright (c) 1997, 1998, 1999, 2003, 2014, 2015 by
! Peter Bird and the Regents of the University of California.
!-----------------------------------------------------------------
!
! CONTENTS OF THIS MODULE:
! --------------------------
! INTENDED FOR THE USER TO CALL:
! DOUBLE PRECISION FUNCTION DArc
! SUBROUTINE DCircles_Intersect
! LOGICAL FUNCTION DClockways
! DOUBLE PRECISION FUNCTION DCompass
! SUBROUTINE DCross
! DOUBLE PRECISION FUNCTION DDot
! DOUBLE PRECISION FUNCTION DEasting
! DOUBLE PRECISION FUNCTION DLength
! SUBROUTINE DLocal_Phi
! SUBROUTINE DLocal_Theta
! SUBROUTINE DLonLat_2_ThetaPhi
! SUBROUTINE DLonLat_2_Uvec
! DOUBLE PRECISION FUNCTION DMagnitude
! SUBROUTINE DMake_Uvec
! SUBROUTINE DNorthEast_Convention
! DOUBLE PRECISION FUNCTION DRelative_Compass
! SUBROUTINE DSpherical_Area
! SUBROUTINE DThetaPhi_2_LonLat
! SUBROUTINE DThetaPhi_2_Uvec
! SUBROUTINE DTurn_To
! SUBROUTINE DUvec_2_LonLat
! SUBROUTINE DUvec_2_ThetaPhi
! SUBROUTINE DUvec_2_PlungeAzimuth
!
! UTILITY ROUTINE FOR DEBUGGING:
! SUBROUTINE DTraceback
!-----------------------------------------------------------------
!declares
IMPLICIT NONE
DOUBLE PRECISION, PARAMETER :: Pi = 3.14159265358979D0 ! These values were confirmed
DOUBLE PRECISION, PARAMETER :: Pi_over_2 = 1.57079632679490D0 ! with PROGRAM Check_Pi.
DOUBLE PRECISION, PARAMETER :: Two_Pi = 6.28318530717959D0
DOUBLE PRECISION, PARAMETER :: degrees_per_radian = 57.2957795130823D0
DOUBLE PRECISION, PARAMETER :: radians_per_degree = 0.0174532925199433D0
! ---------------------------------------------------------
! | 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
DOUBLE PRECISION FUNCTION DArc (from_uvec, to_uvec)
! Returns length of great-circle arc in radians.
IMPLICIT NONE
DOUBLE PRECISION, DIMENSION(3), INTENT(IN) :: from_uvec, to_uvec
DOUBLE PRECISION :: crossed, dotted
DOUBLE PRECISION, DIMENSION(3) :: t_vec
CALL DCross (from_uvec, to_uvec, t_vec)
crossed = DLength(t_vec) ! >= 0.
dotted = DDot(from_uvec, to_uvec)
DArc = DATAN2(crossed, dotted) ! 0.0D0 to Pi
END FUNCTION DArc
SUBROUTINE DCircles_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.)
!"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.00D0 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*8, DIMENSION(3), INTENT(IN) :: pole_a_uvec, first_a_uvec, last_a_uvec, &
& pole_b_uvec, first_b_uvec, last_b_uvec
!{Following DOUBLE PRECISION was present in the original MODULE Map_Projections.}
DOUBLE PRECISION, INTENT(IN) :: dot_a, dot_b
LOGICAL, INTENT(OUT) :: overlap
INTEGER, INTENT(OUT) :: number
REAL*8, DIMENSION(3), INTENT(OUT) :: point1_uvec, point2_uvec
!{Following DOUBLE PRECISIONs were present in the original MODULE Map_Projections.}
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*8 :: 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*8, PARAMETER :: tolerance = 0.0000000014D0
REAL*8, DIMENSION(3) :: line_uvec, offset, t_vec
! Check for an easy answer (no intersections):
number = 0
overlap = .FALSE.
IF (dot_a >= 1.0D0) RETURN
IF (dot_a <= -1.0D0) RETURN
IF (dot_b >= 1.0D0) RETURN
IF (dot_b <= -1.0D0) 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 = -DRelative_Compass(pole_a_uvec, first_a_uvec)
last_a_radians = -DRelative_Compass(pole_a_uvec, last_a_uvec)
IF (first_a_radians > last_a_radians) last_a_radians = last_a_radians + Two_Pi
END IF
IF (.NOT.ring_b) THEN
first_b_radians = -DRelative_Compass(pole_b_uvec, first_b_uvec)
last_b_radians = -DRelative_Compass(pole_b_uvec, last_b_uvec)
IF (first_b_radians > last_b_radians) last_b_radians = last_b_radians + Two_Pi
END IF
! Test for special cases of parallel and antipodal poles:
t_vec = pole_a_uvec - pole_b_uvec
pole_gap = DLength(t_vec)
t_vec = pole_a_uvec + pole_b_uvec
antipole_gap = DLength(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 = -DRelative_Compass(pole_a_uvec, last_b_uvec)
last_b_radians = -DRelative_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 = DDot(pole_a_uvec, pole_b_uvec)
deterMinant = 1.0D0 - a_dot_b**2
IF (deterMinant == 0.0D0) THEN
WRITE (*,"(' ERROR: DeterMinant = 0.0D0 in DCircles_Intersect')")
CALL DTraceback
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)
min_radius = DLength(offset)
IF (min_radius == 1.0D0) 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 = -DRelative_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 = -DRelative_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 DCross (pole_a_uvec, pole_b_uvec, t_vec)
CALL DMake_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 = DSQRT ( 1.0D0 - 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 = -DRelative_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 = -DRelative_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 = -DRelative_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 = -DRelative_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.0D0, or < 1.0D0 [no ELSE; we RETURN with number = 0]
END IF ! poles parallel, or antipodal, or unrelated
END SUBROUTINE DCircles_Intersect
DOUBLE PRECISION FUNCTION DCompass (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 DRelative_Compass.)
IMPLICIT NONE
DOUBLE PRECISION, DIMENSION(3), INTENT(IN) :: from_uvec, to_uvec
DOUBLE PRECISION :: ve, vn
DOUBLE PRECISION, DIMENSION(3) :: omega_vec, omega_uvec, &
& Phi, Theta, &
& v_vec
IF ((from_uvec(1) == 0.0D0).AND.(from_uvec(2) == 0.0D0)) THEN
WRITE (*,"(' ERROR: Compass undefined at N or S pole. Use Relative_Compass.')")
CALL DTraceback
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 DTraceback
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 DTraceback
ELSE
! Normal case:
CALL DCross (from_uvec, to_uvec, omega_vec)
CALL DMake_Uvec(omega_vec, omega_uvec)
CALL DCross (omega_uvec, from_uvec, v_vec)
CALL DLocal_Theta(from_uvec, Theta)
CALL DLocal_Phi (from_uvec, Phi)
vn = -DDot(v_vec, Theta)
ve = DDot(v_vec, Phi)
DCompass = DATAN2(ve, vn)
END IF
END FUNCTION DCompass
SUBROUTINE DCross (a_vec, b_vec, c_vec)
! vector cross product: a x b = c
IMPLICIT NONE
DOUBLE PRECISION, DIMENSION(3), INTENT(IN) :: a_vec, b_vec
DOUBLE PRECISION, 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 DCross
DOUBLE PRECISION FUNCTION DDot (a_vec, b_vec)
! returns scalar (dot) product of two 3-component vectors
IMPLICIT NONE
DOUBLE PRECISION, DIMENSION(3), INTENT(IN) :: a_vec, b_vec
DDot = a_vec(1)*b_vec(1) + a_vec(2)*b_vec(2) + a_vec(3)*b_vec(3)
END FUNCTION DDot
DOUBLE PRECISION FUNCTION DEasting(delta_lon_degrees)
!returns positive result, 0.0D0 ~ 359.999...D0
IMPLICIT NONE
DOUBLE PRECISION, INTENT(IN) :: delta_lon_degrees
DEasting = DMOD((delta_lon_degrees + 720.0D0),360.0D0)
END FUNCTION DEasting
INTEGER FUNCTION IAbove(x)
! returns first integer >= x, unlike INT() or NINT():
IMPLICIT NONE
DOUBLE PRECISION, INTENT(IN) :: x
INTEGER answer
DOUBLE PRECISION :: t
answer = INT(x)
IF (x > 0.0D0) THEN
t = 1.00D0 * answer
IF (x > t) THEN
answer = answer + 1
END IF
END IF
IAbove = answer
END FUNCTION IAbove
DOUBLE PRECISION FUNCTION DLength(a_vec)
IMPLICIT NONE
DOUBLE PRECISION, DIMENSION(3), INTENT(IN) :: a_vec
DOUBLE PRECISION :: t
t = a_vec(1)**2 + &
& a_vec(2)**2 + &
& a_vec(3)**2
IF (t == 0.0D0) THEN
DLength = 0.0D0
ELSE
DLength = DSQRT(t)
END IF
END FUNCTION DLength
SUBROUTINE DLocal_Phi (b_, Phi)
! returns local East-pointing unit vector in Cartesian coordinates
! for location b_; not intended to work at the poles!
DOUBLE PRECISION, DIMENSION(3), INTENT(IN) :: b_
DOUBLE PRECISION, DIMENSION(3), INTENT(OUT) :: Phi
DOUBLE PRECISION, DIMENSION(3) :: temp
IF (b_(1) == 0.0D0) THEN
IF (b_(2) == 0.0D0) THEN
WRITE (*,"(' ERROR: DLocal_Phi was requested for N or S pole.')")
CALL DTraceback
END IF
END IF
temp(1) = -b_(2)
temp(2) = b_(1)
temp(3) = 0.0D0
CALL DMake_Uvec(temp, Phi)
END SUBROUTINE DLocal_Phi
SUBROUTINE DLocal_Theta (b_, Theta)
! returns local South-pointing unit vector in Cartesian coordinates
! for location b_; not intended to work at the poles!
DOUBLE PRECISION, DIMENSION(3), INTENT(IN) :: b_
DOUBLE PRECISION, DIMENSION(3), INTENT(OUT) :: Theta
DOUBLE PRECISION, DIMENSION(3) :: temp
DOUBLE PRECISION :: equat, new_equat
equat = DSQRT(b_(1)**2 + b_(2)**2) ! equatorial component
IF (equat == 0.0D0) THEN
WRITE (*,"(' ERROR: DLocal_Theta was requested for N or S pole.')")
CALL DTraceback
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 DMake_Uvec (temp, Theta)
END SUBROUTINE DLocal_Theta
SUBROUTINE DLonLat_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
DOUBLE PRECISION, INTENT(IN) :: lon, lat
DOUBLE PRECISION, INTENT(OUT) :: theta, phi
IF ((lat > 90.0D0).OR.(lat < -90.0D0)) THEN
WRITE (*,"(' ERROR: Latitude outside legal range.')")
CALL DTraceback
ELSE
theta = radians_per_degree * (90.0D0 - lat)
phi = radians_per_degree * lon
END IF
END SUBROUTINE DLonLat_2_ThetaPhi
SUBROUTINE DLonLat_2_Uvec (lon, lat, uvec)
IMPLICIT NONE
DOUBLE PRECISION, INTENT(IN) :: lon, lat
DOUBLE PRECISION, DIMENSION(3), INTENT(OUT) :: uvec
DOUBLE PRECISION :: theta, phi
CALL DLonLat_2_ThetaPhi (lon, lat, theta, phi)
CALL DThetaPhi_2_Uvec (theta, phi, uvec)
END SUBROUTINE DLonLat_2_Uvec
DOUBLE PRECISION FUNCTION DMagnitude (b_)
DOUBLE PRECISION, DIMENSION(3), INTENT(IN) :: b_
DMagnitude = DSQRT(b_(1)**2 +b_(2)**2 + b_(3)**2)
END FUNCTION DMagnitude
!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 DMake_Uvec (vector, uvec)
! Shortens or lengthens a three-component vector to a unit vector.
IMPLICIT NONE
DOUBLE PRECISION, DIMENSION(3), INTENT(IN) :: vector
DOUBLE PRECISION, DIMENSION(3), INTENT(OUT):: uvec
DOUBLE PRECISION :: factor, size
size = DLength(vector)
IF (size > 0.0D0) THEN
factor = 1.0D0 / size
uvec = vector * factor
ELSE
WRITE (*,"(' ERROR: Cannot DMake_Uvec of (0.0D0, 0.0D0, 0.0D0).')")
CALL DTraceback
END IF
END SUBROUTINE DMake_Uvec
SUBROUTINE DNorthEast_Convention (location_uvec, north_uvec, east_uvec)
! At most positions ("location_uvec") returns "north_uvec"
! (same as -DLocal_Theta) and "east_uvec" (same as +DLocal_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 DRelative_Compass
! and SUBROUTINE DTurn_To call this routine, they can work
! together even when location_uvec is at, or near, a pole.
IMPLICIT NONE
DOUBLE PRECISION, DIMENSION(3), INTENT(IN) :: location_uvec
DOUBLE PRECISION, DIMENSION(3), INTENT(OUT) :: north_uvec, east_uvec
DOUBLE PRECISION, PARAMETER :: pole_area = 7.615D-7 ! (0.05 degrees -> radians -> squared)
DOUBLE PRECISION :: equat2
DOUBLE PRECISION, 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.0D0
CALL DMake_Uvec (t_vec, east_uvec)
ELSE
! very close to N or S pole; act as if on 0E meridian:
east_uvec = (/ 0.0D0, 1.0D0, 0.0D0 /)
END IF
CALL DCross (location_uvec, east_uvec, north_uvec)
END SUBROUTINE DNorthEast_Convention
DOUBLE PRECISION FUNCTION DRelative_Compass (from_uvec, to_uvec)
! At most points, works exactly like DCompass:
! 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 DCompass, it does not crash at the N and S poles,
! where Theta and Phi are undefined. Instead, it uses
! SUBROUTINE DNorthEast_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
DOUBLE PRECISION, DIMENSION(3), INTENT(IN) :: from_uvec, to_uvec
DOUBLE PRECISION :: ve, vn
DOUBLE PRECISION, 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 DTraceback
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 DTraceback
ELSE
! Normal case:
CALL DCross (from_uvec, to_uvec, omega_vec)
CALL DMake_Uvec(omega_vec, omega_uvec)
CALL DCross (omega_uvec, from_uvec, v_vec)
CALL DNorthEast_Convention (from_uvec, north_uvec, east_uvec)
vn = DDot(v_vec, north_uvec)
ve = DDot(v_vec, east_uvec)
DRelative_Compass = DATAN2(ve, vn)
END IF
END FUNCTION DRelative_Compass
SUBROUTINE DSpherical_Area(vertex1, vertex2, vertex3, steradians)
! Given 3 Cartesian (x, y, z) unit-vectors,
! each from the center of a unit-sphere to its surface,
! defining 3 surface points (vertex1, vertex2, vertex3),
! returns the area of the spherical triangle in steradians.
!{For dimensional area, multiply this by the square of the radius.}
! Reference is:
! I. Todhunter [1886]
! Spherical Trigonometry
! McMillan & Co., London {5th edition};
! http://www.gutenberg.org/ebooks/19770 ,
! especially sections 97 (large triangles) and
! section 109 (small triangles).
! Additional feature added in this code (only):
! If vertices are specified in counterclockwise order
!(as seen from outside the sphere), then the area is positive;
! if they are specified in clockwise order, then the area is
! returned as negative. If the 3 unit-vectors are coplanar,
! or if at least 2 of them are identical, then the area is zero.
! For testing of the ideal cross-over point between answer1 and answer2,
! see my files Test_DSpherical_Area.f90, ...txt, ...xlsx .
IMPLICIT NONE
REAL*8, DIMENSION(3), INTENT(IN) :: vertex1, vertex2, vertex3
REAL*8, INTENT(OUT) :: steradians
LOGICAL :: positive
REAL*8 :: A, alpha, answer1, answer2, azim12, azim13, azim21, azim23, azim31, azim32, &
& B, beta, C, dotted, E, fraction, gamma, plane_area
REAL*8, DIMENSION(3) :: crossed, v2_minus_v1, v3_minus_v1
!- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
!Test for sign of answer, including possibility that answer is zero:
v2_minus_v1(1:3) = vertex2(1:3) - vertex1(1:3) ! not a uvec
v3_minus_v1(1:3) = vertex3(1:3) - vertex1(1:3) ! not a uvec
CALL DCross (v2_minus_v1, v3_minus_v1, crossed) ! not a uvec
dotted = DDot(crossed, vertex1)
IF (dotted == 0.0D0) THEN
steradians = 0.0D0
RETURN
ELSE
positive = (dotted > 0.0D0)
END IF
!- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
!Method for large triangles:
! (All azimuths are in radians, clockwise from North.)
azim12 = DRelative_Compass (from_uvec = vertex1, to_uvec = vertex2)
azim13 = DRelative_Compass (from_uvec = vertex1, to_uvec = vertex3)
azim21 = DRelative_Compass (from_uvec = vertex2, to_uvec = vertex1)
azim23 = DRelative_Compass (from_uvec = vertex2, to_uvec = vertex3)
azim31 = DRelative_Compass (from_uvec = vertex3, to_uvec = vertex1)
azim32 = DRelative_Compass (from_uvec = vertex3, to_uvec = vertex2)
A = ABS(azim12 - azim13)
IF (A > Pi) A = Two_Pi - A ! where Two_Pi is a predefined global
B = ABS(azim21 - azim23)
IF (B > Pi) B = Two_Pi - B
C = ABS(azim31 - azim32)
IF (C > Pi) C = Two_Pi - C
E = A + B + C - Pi ! where Pi is a predefined global
IF (positive) THEN
answer1 = E
ELSE
answer1 = -E
END IF
!- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
!Method for small triangles:
alpha = DArc(vertex1, vertex2) ! N.B. The naming of these 3 side-arcs may not
beta = DArc(vertex2, vertex3) ! necessarily agree with naming in Todhunter[1886](?),
gamma = DArc(vertex3, vertex1) ! but this will make no difference to our answer.
fraction = (alpha**2 + beta**2 + gamma**2) / 24.0D0
plane_area = DMagnitude(crossed) / 2.0D0
IF (positive) THEN
answer2 = plane_area * (1.0D0 + fraction)
ELSE
answer2 = -plane_area * (1.0D0 + fraction)
END IF
!- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
!Selection of final answer:
IF (E > 1.0D-8) THEN ! "large" triangle formula is more accurate:
steradians = answer1
ELSE ! "small" triangle formula is more accurate:
steradians = answer2
END IF
END SUBROUTINE DSpherical_Area
SUBROUTINE DThetaPhi_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
DOUBLE PRECISION, INTENT(IN) :: theta, phi
DOUBLE PRECISION, INTENT(OUT) :: lon, lat
lat = 90.0D0 - degrees_per_radian * DABS(theta)
lat = MAX (lat, -90.0D0)
lon = degrees_per_radian * phi
IF (lon > 180.0D0) lon = lon - 360.0D0
IF (lon <= -180.0D0) lon = lon + 360.0D0
END SUBROUTINE DThetaPhi_2_LonLat
SUBROUTINE DThetaPhi_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
DOUBLE PRECISION, INTENT(IN) :: theta, phi
DOUBLE PRECISION, DIMENSION(3), INTENT(OUT) :: uvec
DOUBLE PRECISION :: equat
uvec(3) = DCOS(theta)
equat = DSIN(theta)
uvec(1) = equat * DCOS(phi)
uvec(2) = equat * DSIN(phi)
END SUBROUTINE DThetaPhi_2_Uvec
SUBROUTINE DTraceback ()
! 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
DOUBLE PRECISION, 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.0D0 * i
END DO
STOP ' '
END SUBROUTINE DTraceback
SUBROUTINE DTurn_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 DNorthEast_Convention as
! FUNCTION DRelative_Compass does, so they can work together
! to find internal points on a small circle (as in DProcess_L4_
! Paths).
IMPLICIT NONE
DOUBLE PRECISION, INTENT(IN) :: azimuth_radians, far_radians
DOUBLE PRECISION, DIMENSION(3), INTENT(IN) :: base_uvec
DOUBLE PRECISION, DIMENSION(3), INTENT(OUT) :: omega_uvec, result_uvec
DOUBLE PRECISION, PARAMETER :: pole_width = 1.745D-4 ! 0.01 degrees; must match DNortheast_Convention!
DOUBLE PRECISION :: complement, cos_size, e_part, n_part, sin_size
DOUBLE PRECISION, DIMENSION(3) :: east_uvec, north_uvec, t_uvec
DOUBLE PRECISION, DIMENSION(3,3) :: rotation_matrix
CALL DNorthEast_Convention (base_uvec, north_uvec, east_uvec)
e_part = -DCOS(azimuth_radians)
n_part = DSIN(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 = DCOS(far_radians)
sin_size = DSIN(far_radians)
complement = 1.00D0 - 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 DTurn_To
SUBROUTINE DUvec_2_LonLat (uvec, lon, lat)
IMPLICIT NONE
DOUBLE PRECISION, DIMENSION(3), INTENT(IN) :: uvec
DOUBLE PRECISION, INTENT(OUT) :: lon, lat
DOUBLE PRECISION :: theta, phi
CALL DUvec_2_ThetaPhi (uvec, theta, phi)
CALL DThetaPhi_2_LonLat (theta, phi, lon, lat)
END SUBROUTINE DUvec_2_LonLat
SUBROUTINE DUvec_2_ThetaPhi (uvec, theta, phi)
! converts from Cartesian unit vector to theta (colatitude)
! and phi (longitude), both in radians
IMPLICIT NONE
DOUBLE PRECISION, DIMENSION(3), INTENT(IN) :: uvec
DOUBLE PRECISION, INTENT(OUT) :: theta, phi
DOUBLE PRECISION :: equat, equat2
equat2 = uvec(1)*uvec(1) + uvec(2)*uvec(2)
IF (equat2 == 0.0D0) THEN
phi = 0.0D0 ! actually undefined; provide default 0.0D0
IF (uvec(3) > 0.0D0) THEN
theta = 0.0D0 ! N pole
ELSE
theta = Pi ! S pole
END IF
ELSE
equat = DSQRT(equat2)
theta = DATAN2(equat, uvec(3))
phi = DATAN2(uvec(2), uvec(1))
END IF
END SUBROUTINE DUvec_2_ThetaPhi
SUBROUTINE DUvec_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
DOUBLE PRECISION, DIMENSION(3), INTENT(IN) :: location_uvec, lineation_uvec
INTEGER, INTENT(OUT) :: plunge_degrees, azimuth_degrees
DOUBLE PRECISION :: equat, horizontal_component, phi_component, theta_component, up_component
DOUBLE PRECISION, DIMENSION(3) :: phi_uvec, theta_uvec
equat = DSQRT(location_uvec(1)*location_uvec(1) + location_uvec(2)*location_uvec(2))
IF (equat == 0.0D0) THEN
WRITE (*,"(' ERROR: location_uvec sent to DUvec_2_PlungeAzimuth may not be N or S pole.')")
CALL DTraceback
END IF
CALL DLocal_Theta(location_uvec, theta_uvec)
CALL DLocal_Phi (location_uvec, phi_uvec)
up_component = DDot(lineation_uvec, location_uvec)
theta_component = DDot(lineation_uvec, theta_uvec)
phi_component = DDot(lineation_uvec, phi_uvec)
horizontal_component = DSQRT(theta_component**2 + phi_component**2)
IF (horizontal_component <= 0.0D0) THEN
plunge_degrees = 90
azimuth_degrees = 0
ELSE IF (up_component <= 0.0D0) THEN
plunge_degrees = NINT(degrees_per_radian * DATAN2(-up_component, horizontal_component))
azimuth_degrees = NINT(degrees_per_radian * DATAN2(phi_component, -theta_component))
IF (azimuth_degrees < 0) azimuth_degrees = azimuth_degrees + 360
ELSE ! up_component is positive
plunge_degrees = NINT(degrees_per_radian * DATAN2(up_component, horizontal_component))
azimuth_degrees = NINT(degrees_per_radian * DATAN2(-phi_component, theta_component))
IF (azimuth_degrees < 0) azimuth_degrees = azimuth_degrees + 360
END IF
END SUBROUTINE DUvec_2_PlungeAzimuth
DOUBLE PRECISION FUNCTION DDot_3D (a_, b_)
! Dot product of 3-component vectors.
DOUBLE PRECISION, DIMENSION(3) :: a_, b_
DDot_3D = a_(1)*b_(1) + a_(2)*b_(2) + a_(3)*b_(3)
END FUNCTION DDot_3D
END MODULE DSphere !===============================================