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 !===============================================