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