MODULE Map_Projections ! basic tools for the map-projection plane ! and/or the surface of a spherical planet: ! Level 3: curves and lines on map-projection plane (larger than planet); ! (x,y) coordinates in meters. ! Level 4: great- and small-circle arcs on planet; after any projection ! cuts and limits. Note that you cannot plot directly ! to this level; you must begin at level 5, 6, or 7. ! (expressed as Cartesian unit vectors from planet center) ! Level 5: great- and small-circle arcs on planet before (and to be subjected ! to) projection cuts and limits. ! (expressed as Cartesian unit vectors from planet center) ! Level 6: great- and small-circle arcs on planet. ! (coordinates theta and phi in radians) ! Level 7: great- and small-circle arcs on planet. ! (coordinates are East-longitude and North-latitude in degrees) ! ! Works with module Adobe_Illustrator to produce graphics files readable by ! Adobe Illustrator for Windows, versions 4, 7, ...? USE Adobe_Illustrator ! ! by Peter Bird, UCLA, May 1997 - September 1999. ! copyright (c) 1997, 1998, 1999 by ! Peter Bird and the Regents of the University of California. !----------------------------------------------------------------- ! ! USAGE NOTES FOR THIS MODULE: ! ------------------------------ ! 1. Always begin by initializing module Adobe_Illustrator, with ! CALL Select_Paper ! CALL Set_Background ! CALL Define_Margins ! CALL Begin_Page ! CALL Set_Window [ optional; Begin_Page computes best window ] ! ! 2. Before plotting any map data, you MUST: ! CALL Set_Zoom (scale_denominator, x_center_meters, & ! & y_center_meters, xy_wrt_page_radians) ! [ OR ] ! CALL Set_Zoom (scale_denominator) ! This call must occur after CALL Set_Window; if you re-call ! Set_Window, then this call must also be repeated, or else ! this module will not "know" that the window has moved! ! The scale_denominator is a large number (1.E6 to 1.E9 ?) ! which is the ratio of the size of features on the projection ! plane (roughly life-sized) to the size of their images ! on the paper map; it is dimensionless. ! If you will be plotting (x,y) data [either from flat-Earth ! simulation programs like LARAMY, FAULTS, or PLATES, or ! simply digitized as (x,y) from a map, and never converted to ! spherical coordinates], then use the remaining parameters to ! center and rotate this data. If you will be plotting only ! data which is in spherical coordinates, then you may omit ! the last 3 parameters and accept the default values of zero, ! in which case the projection point (least-distorted part of ! the map) will be in the center of the window, with North up. ! ! 3. Routines for setting the line color, line dashing, fill ! color and/or fill pattern are found in module Adobe_Illustrator, ! and can be called any time before a path is ended, at which ! time the registered values will take effect. ! ! 4. To draw directly on the map-projection plane (Level 3), ! use a combination of: ! CALL New_L3_Path [ begins a path; must be first ] ! CALL Line_To_L3 [ optional; may be repeated ] ! CALL Curve_To_L3 [ optional; may be repeated ] ! CALL End_L3_Path [ ends the path; must be last ] ! possibly in combination with: ! CALL Circle_on_L3 [ either before or after drawing path(s) ] ! possibly in combination with: ! CALL L3_Text [ either before or after drawing path(s) ] ! ! 5. Higher-level routines are required in two cases: ! (a) You want to plot data or finite-element grids in ! one of the spherical coordinate systems; OR ! (b) You want to display latitude and longitude lines ! overlaying a map of data entered in the (x,y) projection plane. ! In either case, you must also call ONE of the routines ! which defines a map projection. For example: ! CALL Set_Mercator (...) ! [ OR ] ! CALL Set_Conic (...) ! CALL Set_Stereographic (...) ! CALL Set_Orthographic (...) ! ! 6. Once the map projection is defined (above), you may draw ! on the spherical planet surface with a combination of: ! CALL New_L567_Path [ begins a path; must be first ] ! CALL Great_To_L567 [ optional; may be repeated ] ! CALL Small_To_L567 [ optional; may be repeated ] ! CALL End_L567_Path [ ends the path; must be last ] ! possibly in combination with: ! CALL L567_Text [ either before or after drawing path(s) ] ! Note that Great_To_L567 adds a great-circle arc to the path, ! while Small_To_L567 adds a small-circle arc to the path. ! ! 7. Always end by closing the Adobe_Illustrator output file: ! CALL End_Page ! ! CONTENTS OF THIS MODULE: ! -------------------------- ! ! USER SUBROUTINES (in a possible order-of-use): ! SUBROUTINE Set_Zoom [ Always required, first! ] ! SUBROUTINE New_L3_Path ! SUBROUTINE Line_To_L3 ! SUBROUTINE Curve_To_L3 ! SUBROUTINE End_L3_Path [ and/or ] ! SUBROUTINE Circle_on_L3 [ and/or ] ! SUBROUTINE L3_Text ! ---- transition from projection plane (level 3) to ! surface of spherical planet (levels 5,6,7) ------ ! SUBROUTINE Set_Mercator ! [ OR ] SUBROUTINE Set_Lambert_Conformal_Conic ! [ OR ] SUBROUTINE Set_Albers_Equal_Area_Conic ! [ OR ] SUBROUTINE Set_Polyconic ! [ OR ] SUBROUTINE Set_Geometric_Conic ! [ OR ] SUBROUTINE Set_Stereographic ! [ OR ] SUBROUTINE Set_Lambert_Azimuthal_EqualArea ! [ OR ] SUBROUTINE Set_Azimuthal_Equidistant ! [ OR ] SUBROUTINE Set_Orthographic ! [ OR ] SUBROUTINE Set_Gnomonic ! [ use ONE of the above ] ! SUBROUTINE New_L45_Path [(Level 4 for internal use) ] ! SUBROUTINE Great_To_L45 [ Level 5 accepts vectors; ] ! SUBROUTINE Small_Through_L45 ! SUBROUTINE Small_To_L45 ! SUBROUTINE End_L45_Path ! SUBROUTINE L5_Text [(There is no L4_Text.) ] ! SUBROUTINE New_L67_Path [ Level 6 accepts theta/phi;] ! SUBROUTINE Great_To_L67 [ Level 7 accepts lon/lat. ] ! SUBROUTINE Small_To_L67 ! SUBROUTINE End_L67_Path ! SUBROUTINE L67_Text ! ! UTILITY SUBROUTINES AND FUNCTIONS: ! REAL FUNCTION Arc ! SUBROUTINE Circles_Intersect ! LOGICAL FUNCTION Clockways ! REAL FUNCTION Compass ! REAL FUNCTION Conformal_Deflation ! SUBROUTINE Cross ! REAL FUNCTION Dot ! REAL FUNCTION Easting ! CHARACTER*1 FUNCTION In_Rind ! LOGICAL FUNCTION L5_In_Window ! REAL FUNCTION Length ! SUBROUTINE Local_Phi ! SUBROUTINE Local_Theta ! SUBROUTINE LonLat_2_ThetaPhi ! SUBROUTINE LonLat_2_Uvec ! SUBROUTINE Magnitude ! SUBROUTINE Make_Uvec ! SUBROUTINE Meters_2_Points ! SUBROUTINE NorthEast_Convention ! SUBROUTINE Points_2_LonLat ! SUBROUTINE Points_2_Meters ! SUBROUTINE Process_L3_Paths ! SUBROUTINE Process_L3_Text ! SUBROUTINE Process_L4_Paths ! SUBROUTINE Process_L4_Text ! SUBROUTINE Process_L5_Paths ! SUBROUTINE Process_L5_Text ! SUBROUTINE Process_L6_Paths ! SUBROUTINE Process_L6_Text ! SUBROUTINE Process_L7_Paths ! SUBROUTINE Process_L7_Text ! SUBROUTINE Project ! SUBROUTINE Reject ! REAL FUNCTION Relative_Compass ! SUBROUTINE Restore_mp_State ! SUBROUTINE Save_mp_State ! SUBROUTINE ThetaPhi_2_LonLat ! SUBROUTINE ThetaPhi_2_Uvec ! SUBROUTINE Trace_Rind ! SUBROUTINE Turn_To ! SUBROUTINE Uvec_2_LonLat ! SUBROUTINE Uvec_2_ThetaPhi ! REAL FUNCTION Vector_Azimuth !----------------------------------------------------------------- !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 ! Relation of map-projection plane (x,y in meters) to the ! map window on the graphics page: LOGICAL :: mp_zoom_defined = .FALSE. REAL :: mp_scale_denominator, & ! denominator of map scale (8.E8 ?). & mp_x_center_meters, & ! x of center of map, meters (0. ?). & mp_y_center_meters, & ! y of center of map, meters (0. ?). & mp_xy_wrt_page_radians ! counterclockwise rotation ! of the (x,y) system of ! the projection plane with ! respect to the (x,y) graphics ! coordinates of the page, ! in radians (0. ?). ! and convenience values, arrays and factors derived from these: REAL :: mp_meters_per_point REAL, DIMENSION(2,2) :: mp_m_2_pts_matrix, mp_pts_2_m_matrix REAL, DIMENSION(2) :: mp_m_2_pts_vector, mp_pts_2_m_vector ! Memory of pen position (to avoid null segments in paths): REAL :: mp_last_L3_x_meters, mp_last_L3_y_meters REAL, DIMENSION(3) :: mp_last_uvec REAL :: mp_last_theta, mp_last_phi, mp_last_lon, mp_last_lat ! Limits (on the sphere) of the area to be projected: ! There are up to 4 limits (top, left, right, bottom), each ! of which can be a small-circle, a great circle, or a half ! of a great circle. Which of these are in use for a ! particular map projection?: LOGICAL :: mp_using_top, mp_using_left, mp_using_right, mp_using_bottom ! These unit vectors ("uvec's"), pointing from the center ! of the (unit) sphere to a surface point, are used ! to define the limits of the projected area on the sphere ! when applying the clipping/windowing that separates Level 5 ! from Level 4. Typically, "bottom" and "top" are latitude ! limits, while "left" and "right" are used to define the cut ! which is needed in some projections. REAL, DIMENSION(3) :: mp_bottom_uvec, mp_right_uvec, & & mp_top_uvec, mp_left_uvec ! --------------------------------------------------------- ! | 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). | ! --------------------------------------------------------- ! These following values are the maximum values permitted for the ! dot product between a unit vector representing a position on ! the sphere and the corresponding mp_xxx_uvec above (if used): DOUBLE PRECISION :: mp_bottom_maxdot, mp_right_maxdot, mp_top_maxdot, mp_left_maxdot ! When a combination of mp_xxx_uvec and mp_xxx_maxdot is used ! to define one limit of the projected area, it is sometimes ! necessary to specify that this limit only applies on "one side" ! of the sphere. (For example, this is true of a "cut" like the ! one involved in the Mercator and conic projections.) The ! following indicator uvec's can be set to enforce this: the ! maximum value of the dot product (mp_xxx_maxdot) only applies ! if the dot product of the position with mp_xxx_applies_uvec ! is positive. (If such one-sided indicators are NOT wanted, ! just "turn them off" by setting mp_xxx_applies_uvec = mp_xxx_uvec.) REAL, DIMENSION(3) :: mp_bottom_applies_uvec, mp_right_applies_uvec, & & mp_top_applies_uvec, mp_left_applies_uvec ! To complete the definition of the "rind" (the orange-rind-shaped ! area of projectable points on the sphere) it is necessary to ! specify two points on each circle (great or small or half) which ! is in use as a limit. These "first" and "last" points are in ! order going counterclockwise around the relevant mp_x_uvec ! (as seen from outside the sphere), and they mark the beginning ! and the end of the parts of that circle which are the innermost ! boundary of the rind. That is, they mark the corners of the rind. ! (However, if one limit circle does not touch the others, then these ! two points should coincide, to show a complete circle.) REAL, DIMENSION(3) :: mp_first_top_uvec, mp_last_top_uvec, & & mp_first_left_uvec, mp_last_left_uvec, & & mp_first_right_uvec, mp_last_right_uvec, & & mp_first_bottom_uvec, mp_last_bottom_uvec ! The following flags are not set by the user or by Set_[projection], ! but are set by Corners() and read by Trace_Rind(). They describe ! the topological connections between the (up to) 4 walls of the rind. INTEGER, DIMENSION(0:1,4) :: mp_next_wall ! Choice of map projection: LOGICAL :: mp_projection_defined = .FALSE. CHARACTER*30 :: mp_projection = 'undefined ' INTEGER :: mp_projection_number = 0 REAL :: mp_radius_meters = 6371000. ! Supply of unit vectors (uvec's) used in different ways ! by different projections. However, all share the definition ! that the "mp_projpoint_uvec" points to the ! "projection point", which is a point of ! tangency between the projection plane and the sphere, ! at which North is upward on the (unrolled) projection sheet. ! This point is the default origin of the (x,y) system of level 3, ! although that may be overridden (see below). REAL, DIMENSION(3) :: mp_pole_1_uvec, mp_pole_2_uvec, mp_pole_3_uvec, & & mp_projpoint_uvec ! Additional projection parameters specific to certain projections: REAL :: mp_proj2_F, mp_proj2_lambda0, mp_proj2_n, mp_proj2_rho0 REAL :: mp_proj3_C, mp_proj3_lambda0, mp_proj3_n, mp_proj3_rho0 REAL :: mp_proj4_lat0_radians, mp_proj4_lon0_radians REAL :: mp_proj5_lambda0, mp_proj5_lat0_radians, mp_proj5_n, mp_proj5_rtan, mp_proj5_ypole ! Definition of chosen (x,y) system (in meters) in the unrolled ! map projection plane: LOGICAL :: mp_user_xy ! = F for default system, T for user system. REAL :: mp_x_projpoint_meters, mp_y_projpoint_meters, mp_y_azimuth_radians ! (Note: the default, "natural" system is all three zero, so that ! the (0.,0.) origin of the 2-D coordinate system (for level 3) ! is at the projection point (mp_pole_1_uvec) with North parallel ! to +y at that point. However, these degrees-of-freedom are ! provided in case the user has selected another (x,y) coordinate ! definition that is unrelated to the map projection in use. ! Note that with non-zero values of these 3 parameters (and also ! of mp_x_center_meters, mp_y_center_meters, and mp_xy_wrt_page_radians), ! it is possible to plot data created in any (x,y) system and ! to rotate and move the plot as wanted in the map window, and ! still plot accurate longitude and latitude lines over it! ! Convenience matrices and vectors derived from these: REAL, DIMENSION(2,2) :: mp_pp_2_xy_matrix, mp_xy_2_pp_matrix REAL, DIMENSION(2) :: mp_pp_2_xy_vector, mp_xy_2_pp_vector ! Upright and downright are used for coordinate rotation ! internal to some projections like Mercator and the ! 4 varieties of Conic: REAL, DIMENSION(2,2) :: mp_upright, mp_downright !============= variables for saving state of module Map_Projections (mp_): LOGICAL :: mp_state_saved = .FALSE. ! until CALL Save_mp_State() ! (The following are listed in same order as above, but without the comments) LOGICAL :: sp_zoom_defined REAL :: sp_scale_denominator, & & sp_x_center_meters, & & sp_y_center_meters, & & sp_xy_wrt_page_radians REAL :: sp_meters_per_point REAL, DIMENSION(2,2) :: sp_m_2_pts_matrix, sp_pts_2_m_matrix REAL, DIMENSION(2) :: sp_m_2_pts_vector, sp_pts_2_m_vector REAL :: sp_last_L3_x_meters, sp_last_L3_y_meters REAL, DIMENSION(3) :: sp_last_uvec REAL :: sp_last_theta, sp_last_phi, sp_last_lon, sp_last_lat LOGICAL :: sp_using_top, sp_using_left, sp_using_right, sp_using_bottom REAL, DIMENSION(3) :: sp_bottom_uvec, sp_right_uvec, & & sp_top_uvec, sp_left_uvec REAL :: sp_bottom_maxdot, sp_right_maxdot, sp_top_maxdot, sp_left_maxdot REAL, DIMENSION(3) :: sp_bottom_applies_uvec, sp_right_applies_uvec, & & sp_top_applies_uvec, sp_left_applies_uvec REAL, DIMENSION(3) :: sp_first_top_uvec, sp_last_top_uvec, & & sp_first_left_uvec, sp_last_left_uvec, & & sp_first_right_uvec, sp_last_right_uvec, & & sp_first_bottom_uvec, sp_last_bottom_uvec INTEGER, DIMENSION(0:1,4) :: sp_next_wall LOGICAL :: sp_projection_defined = .FALSE. CHARACTER*30 :: sp_projection = 'undefined ' INTEGER :: sp_projection_number = 0 REAL :: sp_radius_meters = 6371000. REAL, DIMENSION(3) :: sp_pole_1_uvec, sp_pole_2_uvec, sp_pole_3_uvec, & & sp_projpoint_uvec REAL :: sp_proj2_F, sp_proj2_lambda0, sp_proj2_n, sp_proj2_rho0 REAL :: sp_proj3_C, sp_proj3_lambda0, sp_proj3_n, sp_proj3_rho0 REAL :: sp_proj4_lat0_radians, sp_proj4_lon0_radians REAL :: sp_proj5_lambda0, sp_proj5_lat0_radians, sp_proj5_n, sp_proj5_rtan, sp_proj5_ypole LOGICAL :: sp_user_xy REAL :: sp_x_projpoint_meters, sp_y_projpoint_meters, sp_y_azimuth_radians REAL, DIMENSION(2,2) :: sp_pp_2_xy_matrix, sp_xy_2_pp_matrix REAL, DIMENSION(2) :: sp_pp_2_xy_vector, sp_xy_2_pp_vector REAL, DIMENSION(2,2) :: sp_upright, sp_downright !============= end of mp_ - saving variables ========================= CONTAINS ! =============================================== ! | USER ROUTINES | ! | (in the typical order of use) | ! =============================================== SUBROUTINE Set_Zoom (scale_denominator, x_center_meters, & & y_center_meters, xy_wrt_page_radians) ! Defines scaling, rotation, and centering of the (planet-sized) ! projection plane in the (page-sized) map window. ! Scale denominator is large and dimensionless (e.g., 1.E6 ~ 1.E9). ! The point (x_center_meters,y_center_meters) will be placed ! in the center of the map window (read from global data ! in module Adobe_Illustrator). ! The +x axis will appear rotated xy_wrt_page_radians ! counterclockwise with respect to the bottom edge of the page. ! Only the scale_denominator argument is required, the rest ! will default to zero (suggested values for plotting data ! that is in one of the spherical coordinate systems). IMPLICIT NONE REAL, INTENT(IN) :: scale_denominator REAL, INTENT(IN), OPTIONAL :: x_center_meters, & & y_center_meters, xy_wrt_page_radians REAL :: factor ! Store values in module-level variables: mp_zoom_defined = .TRUE. mp_scale_denominator = scale_denominator mp_meters_per_point = scale_denominator / (72. * 39.370079) IF (PRESENT(x_center_meters)) THEN mp_x_center_meters = x_center_meters ELSE mp_x_center_meters = 0.0 END IF IF (PRESENT(y_center_meters)) THEN mp_y_center_meters = y_center_meters ELSE mp_y_center_meters = 0. END IF IF (PRESENT(xy_wrt_page_radians)) THEN mp_xy_wrt_page_radians = xy_wrt_page_radians ELSE mp_xy_wrt_page_radians = 0. END IF ! Pre-compute meters -> points conversion matrix and vector: factor = 72. * 39.370079 / mp_scale_denominator mp_m_2_pts_matrix(1,1) = factor * COS(mp_xy_wrt_page_radians) mp_m_2_pts_matrix(2,2) = mp_m_2_pts_matrix(1,1) mp_m_2_pts_matrix(2,1) = factor * SIN(mp_xy_wrt_page_radians) mp_m_2_pts_matrix(1,2) = -mp_m_2_pts_matrix(2,1) mp_m_2_pts_vector(1) = ai_window_xc_points - & & mp_m_2_pts_matrix(1,1) * mp_x_center_meters - & & mp_m_2_pts_matrix(1,2) * mp_y_center_meters mp_m_2_pts_vector(2) = ai_window_yc_points - & & mp_m_2_pts_matrix(2,1) * mp_x_center_meters - & & mp_m_2_pts_matrix(2,2) * mp_y_center_meters ! Pre-compute points -> meters conversion matrix and vector: factor = mp_scale_denominator / (72. * 39.370079) mp_pts_2_m_matrix(1,1) = factor * COS(mp_xy_wrt_page_radians) mp_pts_2_m_matrix(2,2) = mp_pts_2_m_matrix(1,1) mp_pts_2_m_matrix(1,2) = factor * SIN(mp_xy_wrt_page_radians) mp_pts_2_m_matrix(2,1) = -mp_pts_2_m_matrix(1,2) mp_pts_2_m_vector(1) = mp_x_center_meters - & & mp_pts_2_m_matrix(1,1) * ai_window_xc_points - & & mp_pts_2_m_matrix(1,2) * ai_window_yc_points mp_pts_2_m_vector(2) = mp_y_center_meters - & & mp_pts_2_m_matrix(2,1) * ai_window_xc_points - & & mp_pts_2_m_matrix(2,2) * ai_window_yc_points END SUBROUTINE Set_Zoom SUBROUTINE New_L3_Path (x_meters, y_meters) ! Level 3 is data on the map-projection plane, which is ! an imaginary planet-sized plane with dimensions in meters. IMPLICIT NONE REAL, INTENT(IN) :: x_meters, y_meters INTEGER :: path IF (ai_page_open) THEN IF (mp_zoom_defined) THEN IF (.NOT. ai_in_path) THEN ! USUAL CASE: ai_in_path = .TRUE. path = Next_Free_Path() ! Note: this function will check for full library. ai_current_path_index = path ai_current_path_level = 3 ai_path_level(path) = 3 ai_Ln_paths(3) = ai_Ln_paths(3) + 1 ai_total_paths = ai_total_paths + 1 ai_segments(path) = 0 ai_pathlib(1:6,0,path) = 0.0 ai_pathlib(1,0,path) = x_meters ai_pathlib(2,0,path) = y_meters mp_last_L3_x_meters = x_meters ! (memory) mp_last_L3_y_meters = y_meters ELSE ! ai_in_path already WRITE (*,"(' ERROR: Cannot New_L3_Path with a path open.')") CALL Traceback END IF ELSE ! zoom not defined WRITE (*,"(' ERROR: Must Set_Zoom before New_L3_Path.')") CALL Traceback END IF ELSE ! no page open WRITE (*,"(' ERROR: Cannot New_L3_Path before Begin_Page.')") CALL Traceback END IF END SUBROUTINE New_L3_Path SUBROUTINE Line_To_L3 (x_meters, y_meters) IMPLICIT NONE REAL, INTENT(IN) :: x_meters, y_meters INTEGER :: next, path IF (ai_in_path) THEN IF (ai_current_path_level == 3) THEN path = ai_current_path_index IF (ai_segments(path) < ai_longest) THEN ! ignore degenerate lineto's (of zero length): IF (x_meters == mp_last_L3_x_meters) THEN IF (y_meters == mp_last_L3_y_meters) RETURN END IF ! USUAL CASE: next = ai_segments(path) + 1 ai_segments(path) = next ai_pathlib(1:6,next,path) = 0.0 ai_pathlib(1,next,path) = x_meters ai_pathlib(2,next,path) = y_meters ai_bent(next,path) = .FALSE. mp_last_L3_x_meters = x_meters ! (memory) mp_last_L3_y_meters = y_meters ELSE WRITE (*,"(' ERROR: Increase ai_longest from current ',I6)") ai_longest CALL Traceback END IF ELSE WRITE (*,"(' ERROR: Cannot Line_To_L3 in path on level ',I2)") & & ai_current_path_level CALL Traceback END IF ELSE WRITE (*,"(' ERROR: Cannot Line_To_L3 before New_L3_Path.')") CALL Traceback END IF END SUBROUTINE Line_To_L3 SUBROUTINE Curve_To_L3 (x1_meters,y1_meters, x2_meters,y2_meters, & & x3_meters,y3_meters) ! The initial control point is the current point (not input). ! The Bezier handle on the initial control point is (x1,y1). ! The Bezier handle on the final control point is (x2,y2). ! The final control point is (x3,y3). IMPLICIT NONE REAL, INTENT(IN) :: x1_meters,y1_meters, x2_meters,y2_meters, & & x3_meters,y3_meters INTEGER :: next, path IF (ai_in_path) THEN IF (ai_current_path_level == 3) THEN path = ai_current_path_index IF (ai_segments(path) < ai_longest) THEN next = ai_segments(path) + 1 ai_segments(path) = next ai_pathlib(1,next,path) = x1_meters ai_pathlib(2,next,path) = y1_meters ai_pathlib(3,next,path) = x2_meters ai_pathlib(4,next,path) = y2_meters ai_pathlib(5,next,path) = x3_meters ai_pathlib(6,next,path) = y3_meters ai_bent(next,path) = .TRUE. mp_last_L3_x_meters = x3_meters ! (memory) mp_last_L3_y_meters = y3_meters ELSE WRITE (*,"(' ERROR: Increase ai_longest from current ',I6)") ai_longest CALL Traceback END IF ELSE WRITE (*,"(' ERROR: Cannot Curve_To_L3 in path on level ',I2)") & & ai_current_path_level CALL Traceback END IF ELSE WRITE (*,"(' ERROR: Cannot Curve_To_L3 before New_L3_Path.')") CALL Traceback END IF END SUBROUTINE Curve_To_L3 SUBROUTINE End_L3_Path (close, stroke, fill) IMPLICIT NONE LOGICAL, INTENT(IN) :: close, stroke, fill INTEGER :: level, path ! Note: The .AI macros which include "closepath" ! don't close it as one would like, with an extra ! "lineto". Instead, they do a curveto, which somehow ! misses the last point you specified. ! Therefore, it is not enough to specify close = T, ! you must also program so that the path actually ! comes back to its initial point! IF (ai_in_path) THEN path = ai_current_path_index level = ai_current_path_level IF (level == 3) THEN IF (stroke .OR. fill) THEN ! USUAL ROUTE: IF (fill) THEN ! Logic of Process_L2_Paths requires these (and higher-level) ! paths to end at their initial points: IF ((mp_last_L3_x_meters /= ai_pathlib(1,0,path)).OR. & &(mp_last_L3_y_meters /= ai_pathlib(2,0,path))) THEN CALL Line_To_L3(ai_pathlib(1,0,path),ai_pathlib(2,0,path)) END IF ! completion of L3 path needed END IF ! filled L2 path ai_closed(path) = close ai_stroked(path) = stroke ai_filled(path) = fill ai_in_path = .FALSE. mp_last_L3_x_meters = -9.E25 ! (undefined) mp_last_L3_y_meters = -9.E25 CALL Process_L3_Paths ELSE !neither stroked nor filled WRITE (*,"(' ERROR: End_L3_path with neither fill nor stroke.')") CALL Traceback END IF ELSE WRITE (*,"(' ERROR: Cannot End_L3_path when other level open:',I3)") level CALL Traceback END IF ELSE WRITE (*,"(' ERROR: Cannot End_L3_Path when path not open.')") CALL Traceback END IF END SUBROUTINE End_L3_Path SUBROUTINE Circle_on_L3 (x_meters,y_meters, radius_meters, stroke, fill) ! Draws a complete circle centered at (x_meters,y_meters), ! with radius "radius_meters". ! Note that stroke width and color (if used) and ! fill color or pattern (if used) must be predefined. IMPLICIT NONE REAL, INTENT(IN) :: x_meters, y_meters, radius_meters LOGICAL, INTENT(IN) :: stroke, fill REAL :: radius_points, x_points, y_points IF (.NOT.ai_page_open) THEN WRITE (*,"(' ERROR: Cannot Circle_on_L3 before Begin_Page.')") CALL Traceback END IF IF (.NOT.mp_zoom_defined) THEN WRITE (*,"(' ERROR: Cannot Circle_on_L3 before Set_Zoom.')") CALL Traceback END IF IF (ai_in_path) THEN WRITE (*,"(' ERROR: Cannot Circle_on_L3 with another path already open.')") CALL Traceback END IF IF (radius_meters < 0.0) THEN WRITE (*,"(' ERROR: Negative radius to Circle_on_L3: ',1P,E9.2)") radius_meters CALL Traceback END IF CALL Meters_2_Points (x_meters,y_meters, x_points,y_points) radius_points = radius_meters / mp_scale_denominator CALL Circle_on_L12 (2, x_points,y_points, radius_points, stroke, fill) END SUBROUTINE Circle_on_L3 SUBROUTINE L3_Text (x_meters, y_meters, angle_radians, from_x, & & font_points, lr_fraction, ud_fraction, & & text) ! Accepts text strings at level 3 (on the planet-sized ! projection plane, before scaling/rotation/windowing. ! Note that level-3 text is either entirely plotted or ! entirely omitted, based on whether the fiducial point ! (for which coordinates are given) is in the window or not, ! after the conversion from meters to page points. ! Angle_radians is the angle of the baseline of the text, ! measured in radians, counterclockwise. ! IF (from_x) then angle_radians is measured from the +x axis ! in the planet-sized projection plane; ! otherwise it is measured from horizontal on the page. ! Font_points is the font size in points (1/72 inch); remember ! that this is NOT the I-height, but the spacing between ! consecutive lines of text (if it includes CR's). ! An integer value must be chosen. ! Lr_fraction (left/right fraction) is the relative position of ! the fiducial point in the text string: 0. for left, 0.5 ! for middle, or 1. for right-end. Values outside this range ! will also work, within reason. Note that alignment is ! most accurate at values of 0., 0.5, or 1., because Adobe ! Illustrator has better algorithms to estimate the length ! of the text string than I do! ! Ud_fraction (up/down fraction) is the relative position of ! the fiducial point from base to top of characters: ! 0. gives alignment with base of characters (normal); ! -0.4 puts fiducial point below baseline, roughly at the ! base of the "tails" of y, p, g, j; this is useful if ! there is a stroked line at the y_points of the fiducial ! point, and you want the text to clear. ! +1.0 gives fiducial point over the tops of the capitals ! (actually on the baseline of an imaginary preceding line ! of text). ! Text is the character string. ! The length of text is automatically detected; trailing ! blanks are ignored. ! IMPLICIT NONE INTEGER, INTENT(IN) :: font_points REAL, INTENT(IN) :: x_meters, y_meters, angle_radians, & & lr_fraction, ud_fraction LOGICAL, INTENT(IN) :: from_x CHARACTER*(*),INTENT(IN) :: text INTEGER :: bytes REAL :: page_angle_radians bytes = LEN_TRIM(text) IF (bytes < 1) RETURN IF (from_x) THEN page_angle_radians = angle_radians + mp_xy_wrt_page_radians ELSE page_angle_radians = angle_radians END IF CALL Process_L3_Text (x_meters, y_meters, page_angle_radians, & & font_points, lr_fraction, ud_fraction, & & text) END SUBROUTINE L3_Text SUBROUTINE L5_Text (uvec, angle_radians, from_east, & & font_points, lr_fraction, ud_fraction, & & text) ! Accepts text strings at level 5 (on the spherical surface, ! with position expressed as a 3-component unit vector). ! Note that level-5 text is either entirely plotted or ! entirely omitted, based on whether the fiducial point ! (for which unit-vector is given) is in the plotted "rind" ! of the planet, and also in the map window, or not. ! Angle_radians is the angle of the baseline of the text, ! measured in radians, counterclockwise. ! IF (from_east) then angle_radians is measured from local East, ! otherwise it is measured from horizontal on the page. ! Font_points is the font size in points (1/72 inch); remember ! that this is NOT the I-height, but the spacing between ! consecutive lines of text (if it includes CR's). ! An integer value must be chosen. ! Lr_fraction (left/right fraction) is the relative position of ! the fiducial point in the text string: 0. for left, 0.5 ! for middle, or 1. for right-end. Values outside this range ! will also work, within reason. Note that alignment is ! most accurate at values of 0., 0.5, or 1., because Adobe ! Illustrator has better algorithms to estimate the length ! of the text string than I do! ! Ud_fraction (up/down fraction) is the relative position of ! the fiducial point from base to top of characters: ! 0. gives alignment with base of characters (normal); ! -0.4 puts fiducial point below baseline, roughly at the ! base of the "tails" of y, p, g, j; this is useful if ! there is a stroked line at the y_points of the fiducial ! point, and you want the text to clear. ! +1.0 gives fiducial point over the tops of the capitals ! (actually on the baseline of an imaginary preceding line ! of text). ! Text is the character string. ! The length of text is automatically detected; trailing blanks ! are ignored. ! IMPLICIT NONE REAL, DIMENSION(3), INTENT(IN) :: uvec ! unit vector = position INTEGER, INTENT(IN) :: font_points REAL, INTENT(IN) :: angle_radians, & & lr_fraction, ud_fraction LOGICAL, INTENT(IN) :: from_east CHARACTER*(*), INTENT(IN) :: text INTEGER :: bytes bytes = LEN_TRIM(text) IF (bytes < 1) RETURN CALL Process_L5_Text (uvec, angle_radians, from_east, & & font_points, lr_fraction, ud_fraction, & & text) END SUBROUTINE L5_Text SUBROUTINE L67_Text (level, r1, r2, angle_radians, from_east, & & font_points, lr_fraction, ud_fraction, & & text) ! Accepts text strings at level 6 or 7 (on the spherical surface, ! with position expressed as either: ! Level 6: r1 = theta (colatitude from North, radians); ! r2 = phi (East longitude in radians); OR, ! Level 7: r1 = lon (East longitude in degrees); ! r2 = lat (North latitude in degrees). ! Note that level-6/7 text is either entirely plotted or ! entirely omitted, based on whether the fiducial point ! (for which unit-vector is given) is in the plotted "rind" ! of the planet, and also in the map window, or not. ! Angle_radians is the angle of the baseline of the text, ! measured in radians, counterclockwise. ! IF (from_east) then angle_radians is measured from local East, ! otherwise it is measured from horizontal on the page. ! Font_points is the font size in points (1/72 inch); remember ! that this is NOT the I-height, but the spacing between ! consecutive lines of text (if it includes CR's). ! An integer value must be chosen. ! Lr_fraction (left/right fraction) is the relative position of ! the fiducial point in the text string: 0. for left, 0.5 ! for middle, or 1. for right-end. Values outside this range ! will also work, within reason. Note that alignment is ! most accurate at values of 0., 0.5, or 1., because Adobe ! Illustrator has better algorithms to estimate the length ! of the text string than I do! ! Ud_fraction (up/down fraction) is the relative position of ! the fiducial point from base to top of characters: ! 0. gives alignment with base of characters (normal); ! -0.4 puts fiducial point below baseline, roughly at the ! base of the "tails" of y, p, g, j; this is useful if ! there is a stroked line at the y_points of the fiducial ! point, and you want the text to clear. ! +1.0 gives fiducial point over the tops of the capitals ! (actually on the baseline of an imaginary preceding line ! of text). ! Text is the character string. ! The length of text is automatically detected; trailing blanks ! are ignored. ! IMPLICIT NONE INTEGER, INTENT(IN) :: level, font_points REAL, INTENT(IN) :: r1, r2 REAL, INTENT(IN) :: angle_radians, & & lr_fraction, ud_fraction LOGICAL, INTENT(IN) :: from_east CHARACTER*(*), INTENT(IN) :: text INTEGER :: bytes REAL :: lat, lon, phi, theta bytes = LEN_TRIM(text) IF (bytes < 1) RETURN IF (level == 6) THEN theta = r1 phi = r2 IF ((theta < 0.).OR.(theta > Pi)) THEN WRITE (*,"(' ERROR: Illegal colatitude theta for L67_Text.')") CALL Traceback ELSE CALL Process_L6_Text (theta, phi, angle_radians, from_east, & & font_points, lr_fraction, ud_fraction, & & text) END IF ELSE IF (level == 7) THEN lon = r1 lat = r2 IF ((lat > 90.).OR.(lat < -90.)) THEN WRITE (*,"(' ERROR: Illegal latitude for L67_Text.')") CALL Traceback ELSE CALL Process_L7_Text (lon, lat, angle_radians, from_east, & & font_points, lr_fraction, ud_fraction, & & text) END IF ELSE WRITE (*,"(' ERROR: L67_Text cannot accept level ',I2,'.')") level CALL Traceback END IF END SUBROUTINE L67_Text SUBROUTINE Set_Mercator (radius_meters, & & projpoint_uvec, belt_azimuth_radians, & & x_projpoint_meters, y_projpoint_meters, & & y_azimuth_radians) ! First, define an "azimuth" as an angle measured clockwise from ! local North in the local horizontal plane at any point on the ! surface of the sphere. (These azimuths are all in radians!) ! Second, define a "uvec" as a unit vector from the center of the ! sphere to a surface point in a right-handed Cartesian system ! with the first axis at (0 E, 0N), the second at (90 E, 0 N), ! and the third at the North pole (?? E, 90 N). ! (See notes on uvec's in the module data above.) ! - - - - - - - - - - - - - - - - - - - - - - - - ! Projpoint_uvec is a point on the great circle of tangency ! between the sphere and the projection plane (while rolled ! into a cylinder). Not necessarily on the equator!!! ! Belt_azimuth_radians is the azimuth of this great circle ! of tangency at projpoint_uvec. Not necessarily Pi/2 !!! ! The last 3 arguments are provided in case an (x,y) system ! was already created, with no relation to the map projection: ! X_projpoint_meters is the x value desired at projpoint_uvec. ! Y_projpoint_meters is the y value desired at projpoint_uvec. ! Y_azimuth_radians is the azimuth of the +y axis desired at ! point projpoint_uvec. ! The default values for these last 3 arguments are zero. ! These will give the default definition of ! the 2-D Cartesian coordinates in the (unrolled) projection ! plane, which is (0., 0.) at projpoint_uvec, with +y Northward ! and +x Eastward at that point. IMPLICIT NONE REAL, INTENT(IN) :: radius_meters REAL, DIMENSION(3), INTENT(IN) :: projpoint_uvec REAL, INTENT(IN) :: belt_azimuth_radians REAL, INTENT(IN), OPTIONAL :: x_projpoint_meters, y_projpoint_meters, & & y_azimuth_radians REAL :: using_belt_azimuth_radians REAL, DIMENSION(3) :: omega_uvec, t_vec mp_projection_defined = .TRUE. mp_projection = 'Mercator' mp_projection_number = 1 mp_radius_meters = radius_meters CALL Make_Uvec(projpoint_uvec, mp_projpoint_uvec) ! just for sure! ! find the 3 orthogonal uvec's used for Mercator projection: CALL Make_Uvec (projpoint_uvec, mp_pole_1_uvec) ! just for insurance!!! ! Decide whether to invert belt_azimuth_radians !(to prevent whole map appearing upside-down!) IF (SIN(belt_azimuth_radians) == 0.0) THEN using_belt_azimuth_radians = 0.0 ELSE ! use Eastward direction on great circle IF (SIN(belt_azimuth_radians) > 0.0) THEN using_belt_azimuth_radians = belt_azimuth_radians ELSE using_belt_azimuth_radians = belt_azimuth_radians + Pi END IF END IF CALL Turn_To (using_belt_azimuth_radians, mp_pole_1_uvec, Pi_over_2, & ! inputs & omega_uvec, mp_pole_2_uvec) CALL Cross (mp_pole_1_uvec, mp_pole_2_uvec, mp_pole_3_uvec) ! Predefine a 2x2 matrix for the "twist" from ! [ alpha = right (along the tangent belt), beta = up (perp. to tangent) ] ! to "natural" or default (x,y) in the projection plane: mp_upright(1,1) = SIN(using_belt_azimuth_radians) mp_upright(2,2) = mp_upright(1,1) mp_upright(2,1) = COS(using_belt_azimuth_radians) mp_upright(1,2) = -mp_upright(2,1) ! and, the inverse of this matrix is: mp_downright(1,1) = mp_upright(1,1) mp_downright(1,2) = mp_upright(2,1) mp_downright(2,1) = mp_upright(1,2) mp_downright(2,2) = mp_upright(2,2) ! Set up the cut and the high- and low-latitude boundaries: ! High "latitude" (possibly in a rotated coordinate system) limit: mp_using_top = .TRUE. mp_top_uvec = mp_pole_3_uvec mp_top_applies_uvec = mp_top_uvec mp_top_maxdot = SIN(69.98 / degrees_per_radian) ! In a standard Mercator projection, the 70N and 70S parallels ! will have artificial scallops, due to extreme nonlinearity ! of the projection near the pole. We could reduce parameter ! "piece_size" in Process_L4_Paths below 0.174533 (10 degrees), ! but this would expand and slow all plots of everything! ! Instead, we just don't let the 70N and 70S parallels plot. ! Low "latitude" (possibly in a rotated coordinate system) limit: mp_using_bottom = .TRUE. mp_bottom_uvec = -mp_top_uvec ! symmetrical mp_bottom_applies_uvec = mp_bottom_uvec mp_bottom_maxdot = mp_top_maxdot ! symmetrical ! Left ("West-longitude") side of the cut: mp_using_left = .TRUE. t_vec = mp_pole_2_uvec - 0.1 * radians_per_degree * mp_pole_1_uvec CALL Make_Uvec(t_vec, mp_left_uvec) ! Notes: Because of wrap-around, it points "right". ! It is 0.2 degrees of "longitude" away ! from the corresponding "right" limit. t_vec = -mp_pole_2_uvec - 0.05 * radians_per_degree * mp_pole_1_uvec CALL Make_Uvec(t_vec, mp_left_applies_uvec) mp_left_maxdot = 0. ! cut is a great-circle arc ! Right ("East-longitude") side of the cut: mp_using_right = .TRUE. t_vec = -mp_pole_2_uvec - 0.1 * radians_per_degree * mp_pole_1_uvec CALL Make_Uvec(t_vec, mp_right_uvec) ! Notes: Because of wrap-around, it points "right". ! It is 0.2 degrees of "longitude" away ! from the corresponding "right" limit. t_vec = mp_pole_2_uvec - 0.05 * radians_per_degree * mp_pole_1_uvec CALL Make_Uvec(t_vec, mp_right_applies_uvec) mp_right_maxdot = 0. ! cut is a great-circle arc ! Complete description of rind by finding corners: CALL Corners (projpoint_uvec) ! Set up the (x,y) system in the projection plane: mp_user_xy = .FALSE. IF (PRESENT(x_projpoint_meters)) THEN mp_user_xy = .TRUE. mp_x_projpoint_meters = x_projpoint_meters ELSE mp_x_projpoint_meters = 0. END IF IF (PRESENT(y_projpoint_meters)) THEN mp_user_xy = .TRUE. mp_y_projpoint_meters = y_projpoint_meters ELSE mp_y_projpoint_meters = 0. END IF IF (PRESENT(y_azimuth_radians)) THEN mp_user_xy = .TRUE. mp_y_azimuth_radians = y_azimuth_radians ELSE mp_y_azimuth_radians = 0. END IF IF (mp_user_xy) THEN ! from "natural" (x,y) system of projection plane to user system: mp_pp_2_xy_matrix(1,1) = COS(mp_y_azimuth_radians) mp_pp_2_xy_matrix(2,2) = mp_pp_2_xy_matrix(1,1) mp_pp_2_xy_matrix(2,1) = SIN(mp_y_azimuth_radians) mp_pp_2_xy_matrix(1,2) = -mp_pp_2_xy_matrix(2,1) mp_pp_2_xy_vector(1) = mp_x_projpoint_meters mp_pp_2_xy_vector(2) = mp_y_projpoint_meters ! from user (x,y) system to "natural" system of projection plane: mp_xy_2_pp_matrix(1,1) = COS(mp_y_azimuth_radians) mp_xy_2_pp_matrix(2,2) = mp_xy_2_pp_matrix(1,1) mp_xy_2_pp_matrix(1,2) = SIN(mp_y_azimuth_radians) mp_xy_2_pp_matrix(2,1) = -mp_xy_2_pp_matrix(1,2) mp_xy_2_pp_vector(1) = -mp_xy_2_pp_matrix(1,1)*mp_x_projpoint_meters & & -mp_xy_2_pp_matrix(1,2)*mp_y_projpoint_meters mp_xy_2_pp_vector(2) = -mp_xy_2_pp_matrix(2,1)*mp_x_projpoint_meters & & -mp_xy_2_pp_matrix(2,2)*mp_y_projpoint_meters END IF END SUBROUTINE Set_Mercator SUBROUTINE Set_Lambert_Conformal_Conic & & (radius_meters, & & projpoint_uvec, cone_pole_uvec, & & standard_parallel_gap_radians, & & x_projpoint_meters, y_projpoint_meters, & & y_azimuth_radians) ! Radius_meters is the radius of the planet, in meters. ! - - - - - - - - - - - - - - - - - - - - - - - - ! A "uvec" is a unit vector from the center of the ! sphere to a surface point in a right-handed Cartesian system ! with the first axis at (0 E, 0N), the second at (90 E, 0 N), ! and the third at the North pole (?? E, 90 N). ! (See notes on uvec's in the module data above.) ! - - - - - - - - - - - - - - - - - - - - - - - - ! Projpoint_uvec is a point on the small circle of tangency ! between the sphere and the projection plane (while rolled ! into a cone). At this origin point, +x will be Eastward ! and +y will be Northward in the "native" (x,y) of the ! projection plane (as opposed to any user coordinates). ! Cone_pole_uvec is aligned with the axis of the cone, pointing ! toward the closed end. ! Standard_parallel_gap_radians is the difference in latitude ! (in a system centered on cone_pole_uvec as the "North pole") ! between the two standard meridians enclosing projpoint_uvec. ! The last 3 arguments are provided in case an (x,y) system ! was already created, with no relation to the map projection: ! X_projpoint_meters is the x value desired at projpoint_uvec. ! Y_projpoint_meters is the y value desired at projpoint_uvec. ! Y_azimuth_radians is the azimuth of the +y axis desired at ! point projpoint_uvec. ! The default values for these last 3 arguments are zero. ! These will give the default definition of ! the 2-D Cartesian coordinates in the (unrolled) projection ! plane, which is (0., 0.) at projpoint_uvec, with +y Northward ! and +x Eastward at that point. IMPLICIT NONE REAL, INTENT(IN) :: radius_meters, standard_parallel_gap_radians REAL, DIMENSION(3), INTENT(IN) :: projpoint_uvec, cone_pole_uvec REAL, INTENT(IN), OPTIONAL :: x_projpoint_meters, y_projpoint_meters, & & y_azimuth_radians REAL :: latitude_0_radians, latitude_1_radians, latitude_2_radians, & & Pi_over_4, toconepole_azimuth_radians REAL, DIMENSION(3) :: t_vec mp_projection_defined = .TRUE. mp_projection = 'Lambert Conformal Conic' mp_projection_number = 2 mp_radius_meters = radius_meters CALL Make_Uvec(projpoint_uvec, mp_projpoint_uvec) ! just for sure! IF (Arc(mp_projpoint_uvec, cone_pole_uvec) > 1.396) THEN WRITE (*,"(' ERROR: Distance from projection point to pole of cone exceeds 80 degrees.')") CALL Traceback END IF ! find the 3 orthogonal uvec's used for all conic projections: CALL Make_Uvec (cone_pole_uvec, mp_pole_3_uvec) ! just for insurance!!! CALL Cross (mp_pole_3_uvec, mp_projpoint_uvec, t_vec) CALL Make_Uvec (t_vec, mp_pole_2_uvec) CALL Cross (mp_pole_2_uvec, mp_pole_3_uvec, mp_pole_1_uvec) ! Note: Pole_1 is to S of projpoint, if cone_pole is N pole. ! Parameters corresponding to page 105 of John P. Snyder (1983), ! US Geological Survey Bulletin 1532: latitude_0_radians = Pi_over_2 - Arc(mp_projpoint_uvec, cone_pole_uvec) ! Note: "Latitude" is now with respect to conic_pole_uvec as "North pole"! Pi_over_4 = 0.5 * Pi_over_2 IF (standard_parallel_gap_radians > 0.0) THEN latitude_1_radians = latitude_0_radians - 0.5 * standard_parallel_gap_radians latitude_2_radians = latitude_0_radians + 0.5 * standard_parallel_gap_radians IF (latitude_2_radians >= Pi_over_2) THEN WRITE (*,"(' ERROR: Excessive gap between standard parallels.')") CALL Traceback END IF mp_proj2_n = LOG(COS(latitude_1_radians)/COS(latitude_2_radians)) / & & LOG(TAN(Pi_over_4 + latitude_2_radians/2.)/TAN(Pi_over_4 + latitude_1_radians/2.)) mp_proj2_F = COS(latitude_1_radians) * TAN(Pi_over_4 + latitude_1_radians/2.)**mp_proj2_n / mp_proj2_n ELSE ! only one standard parallel, assigned to be latitude_1 = latitude_0. mp_proj2_n = SIN(latitude_0_radians) mp_proj2_F = COS(latitude_0_radians) * TAN(Pi_over_4 + latitude_0_radians/2.)**mp_proj2_n / mp_proj2_n END IF mp_proj2_rho0 = mp_radius_meters * mp_proj2_F / & & TAN(Pi_over_4 + latitude_0_radians/2.)**mp_proj2_n mp_proj2_lambda0 = 0. ! by definition of "longitude" in pole1/pole2/pole3 system. ! Predefine a 2x2 matrix for the "twist" from ! [ up = toward the pole of the cone ] TO [ up = North ], ! the "natural" or default (x,y) in the projection plane: toconepole_azimuth_radians = Relative_Compass(mp_projpoint_uvec, mp_pole_3_uvec) mp_upright(1,1) = COS(toconepole_azimuth_radians) mp_upright(2,2) = mp_upright(1,1) mp_upright(1,2) = SIN(toconepole_azimuth_radians) mp_upright(2,1) = -mp_upright(1,2) ! and, the inverse of this matrix is: mp_downright(1,1) = mp_upright(1,1) mp_downright(1,2) = mp_upright(2,1) mp_downright(2,1) = mp_upright(1,2) mp_downright(2,2) = mp_upright(2,2) ! Set up the cut and the high- and low-latitude boundaries: ! High "latitude" (in the possibly-rotated coordinate system) limit: mp_using_top = (latitude_0_radians < Pi_over_4) ! limit at 45 deg. N from projpoint? IF (mp_using_top) THEN mp_top_uvec = mp_pole_3_uvec mp_top_applies_uvec = mp_top_uvec mp_top_maxdot = SIN(latitude_0_radians + Pi_over_4) ! < 1. END IF ! Low "latitude" (in the possibly-rotated coordinate system) limit: mp_using_bottom = .TRUE. mp_bottom_uvec = -mp_pole_3_uvec mp_bottom_applies_uvec = mp_bottom_uvec mp_bottom_maxdot = MAX(0., -SIN(latitude_0_radians - Pi_over_4)) ! Left ("West-longitude") side of the cut: mp_using_left = .TRUE. t_vec = mp_pole_2_uvec - 0.1 * radians_per_degree * mp_pole_1_uvec CALL Make_Uvec(t_vec, mp_left_uvec) ! Notes: Because of wrap-around, it points "right". ! It is 0.2 degrees of "longitude" away ! from the corresponding "right" limit. t_vec = -mp_pole_2_uvec - 0.05 * radians_per_degree * mp_pole_1_uvec CALL Make_Uvec(t_vec, mp_left_applies_uvec) mp_left_maxdot = 0. ! cut is a great-circle arc ! Right ("East-longitude") side of the cut: mp_using_right = .TRUE. t_vec = -mp_pole_2_uvec - 0.1 * radians_per_degree * mp_pole_1_uvec CALL Make_Uvec(t_vec, mp_right_uvec) ! Notes: Because of wrap-around, it points "right". ! It is 0.2 degrees of "longitude" away ! from the corresponding "right" limit. t_vec = mp_pole_2_uvec - 0.05 * radians_per_degree * mp_pole_1_uvec CALL Make_Uvec(t_vec, mp_right_applies_uvec) mp_right_maxdot = 0. ! cut is a great-circle arc ! Complete description of rind by finding corners: CALL Corners (projpoint_uvec) ! Set up the (x,y) system in the projection plane: mp_user_xy = .FALSE. IF (PRESENT(x_projpoint_meters)) THEN mp_user_xy = .TRUE. mp_x_projpoint_meters = x_projpoint_meters ELSE mp_x_projpoint_meters = 0. END IF IF (PRESENT(y_projpoint_meters)) THEN mp_user_xy = .TRUE. mp_y_projpoint_meters = y_projpoint_meters ELSE mp_y_projpoint_meters = 0. END IF IF (PRESENT(y_azimuth_radians)) THEN mp_user_xy = .TRUE. mp_y_azimuth_radians = y_azimuth_radians ELSE mp_y_azimuth_radians = 0. END IF IF (mp_user_xy) THEN ! from "natural" (x,y) system of projection plane to user system: mp_pp_2_xy_matrix(1,1) = COS(mp_y_azimuth_radians) mp_pp_2_xy_matrix(2,2) = mp_pp_2_xy_matrix(1,1) mp_pp_2_xy_matrix(2,1) = SIN(mp_y_azimuth_radians) mp_pp_2_xy_matrix(1,2) = -mp_pp_2_xy_matrix(2,1) mp_pp_2_xy_vector(1) = mp_x_projpoint_meters mp_pp_2_xy_vector(2) = mp_y_projpoint_meters ! from user (x,y) system to "natural" system of projection plane: mp_xy_2_pp_matrix(1,1) = COS(mp_y_azimuth_radians) mp_xy_2_pp_matrix(2,2) = mp_xy_2_pp_matrix(1,1) mp_xy_2_pp_matrix(1,2) = SIN(mp_y_azimuth_radians) mp_xy_2_pp_matrix(2,1) = -mp_xy_2_pp_matrix(1,2) mp_xy_2_pp_vector(1) = -mp_xy_2_pp_matrix(1,1)*mp_x_projpoint_meters & & -mp_xy_2_pp_matrix(1,2)*mp_y_projpoint_meters mp_xy_2_pp_vector(2) = -mp_xy_2_pp_matrix(2,1)*mp_x_projpoint_meters & & -mp_xy_2_pp_matrix(2,2)*mp_y_projpoint_meters END IF END SUBROUTINE Set_Lambert_Conformal_Conic SUBROUTINE Set_Albers_Equal_Area_Conic & & (radius_meters, & & projpoint_uvec, cone_pole_uvec, & & standard_parallel_gap_radians, & & x_projpoint_meters, y_projpoint_meters, & & y_azimuth_radians) ! Radius_meters is the radius of the planet, in meters. ! - - - - - - - - - - - - - - - - - - - - - - - - ! A "uvec" is a unit vector from the center of the ! sphere to a surface point in a right-handed Cartesian system ! with the first axis at (0 E, 0N), the second at (90 E, 0 N), ! and the third at the North pole (?? E, 90 N). ! (See notes on uvec's in the module data above.) ! - - - - - - - - - - - - - - - - - - - - - - - - ! Projpoint_uvec is a point on the small circle of tangency ! between the sphere and the projection plane (while rolled ! into a cone). At this origin point, +x will be Eastward ! and +y will be Northward in the "native" (x,y) of the ! projection plane (as opposed to any user coordinates). ! Cone_pole_uvec is aligned with the axis of the cone, pointing ! toward the closed end. ! Standard_parallel_gap_radians is the difference in latitude ! (in a system centered on cone_pole_uvec as the "North pole") ! between the two standard meridians enclosing projpoint_uvec. ! The last 3 arguments are provided in case an (x,y) system ! was already created, with no relation to the map projection: ! X_projpoint_meters is the x value desired at projpoint_uvec. ! Y_projpoint_meters is the y value desired at projpoint_uvec. ! Y_azimuth_radians is the azimuth of the +y axis desired at ! point projpoint_uvec. ! The default values for these last 3 arguments are zero. ! These will give the default definition of ! the 2-D Cartesian coordinates in the (unrolled) projection ! plane, which is (0., 0.) at projpoint_uvec, with +y Northward ! and +x Eastward at that point. IMPLICIT NONE REAL, INTENT(IN) :: radius_meters, standard_parallel_gap_radians REAL, DIMENSION(3), INTENT(IN) :: projpoint_uvec, cone_pole_uvec REAL, INTENT(IN), OPTIONAL :: x_projpoint_meters, y_projpoint_meters, & & y_azimuth_radians REAL :: latitude_0_radians, latitude_1_radians, latitude_2_radians, & & Pi_over_4, toconepole_azimuth_radians REAL, DIMENSION(3) :: t_vec mp_projection_defined = .TRUE. mp_projection = 'Albers Equal-Area Conic' mp_projection_number = 3 mp_radius_meters = radius_meters CALL Make_Uvec(projpoint_uvec, mp_projpoint_uvec) ! just for sure! IF (Arc(mp_projpoint_uvec, cone_pole_uvec) > 1.396) THEN WRITE (*,"(' ERROR: Distance from projection point to pole of cone exceeds 80 degrees.')") CALL Traceback END IF ! find the 3 orthogonal uvec's used for all conic projections: CALL Make_Uvec (cone_pole_uvec, mp_pole_3_uvec) ! just for insurance!!! CALL Cross (mp_pole_3_uvec, mp_projpoint_uvec, t_vec) CALL Make_Uvec (t_vec, mp_pole_2_uvec) CALL Cross (mp_pole_2_uvec, mp_pole_3_uvec, mp_pole_1_uvec) ! Note: Pole_1 is to S of projpoint, if cone_pole is N pole. ! Parameters corresponding to pages 95-96 of John P. Snyder (1983), ! US Geological Survey Bulletin 1532: latitude_0_radians = Pi_over_2 - Arc(mp_projpoint_uvec, cone_pole_uvec) ! Note: "Latitude" is now with respect to conic_pole_uvec as "North pole"! Pi_over_4 = 0.5 * Pi_over_2 IF (standard_parallel_gap_radians > 0.0) THEN latitude_1_radians = latitude_0_radians - 0.5 * standard_parallel_gap_radians latitude_2_radians = latitude_0_radians + 0.5 * standard_parallel_gap_radians IF (latitude_2_radians > Pi_over_2) THEN IF (ABS(latitude_2_radians - Pi_over_2) <= 0.01) THEN latitude_2_radians = Pi_over_2 ELSE WRITE (*,"(' ERROR: More poleward standard parallel has latitude over 90.')") CALL Traceback END IF END IF IF (ABS(latitude_1_radians) < 0.01) THEN WRITE (*,"(' ERROR: Neither standard parallel may be on equator.')") CALL Traceback END IF mp_proj3_n = (SIN(latitude_1_radians) + SIN(latitude_2_radians)) / 2. mp_proj3_C = COS(latitude_1_radians)**2 + 2. * mp_proj3_n * SIN(latitude_1_radians) ELSE ! only one standard parallel, assumed to be latitude_1 = latitude_0 (projection point). mp_proj3_n = SIN(latitude_0_radians) mp_proj3_C = COS(latitude_0_radians)**2 + 2. * mp_proj3_n * SIN(latitude_0_radians) END IF mp_proj3_rho0 = mp_radius_meters * SQRT(mp_proj3_C - 2. * mp_proj3_n * SIN(latitude_0_radians)) / mp_proj3_n mp_proj3_lambda0 = 0. ! by definition of "longitude" in pole1/pole2/pole3 system. ! Predefine a 2x2 matrix for the "twist" from ! [ up = toward the pole of the cone ] TO [ up = North ], ! the "natural" or default (x,y) in the projection plane: toconepole_azimuth_radians = Relative_Compass(mp_projpoint_uvec, mp_pole_3_uvec) mp_upright(1,1) = COS(toconepole_azimuth_radians) mp_upright(2,2) = mp_upright(1,1) mp_upright(1,2) = SIN(toconepole_azimuth_radians) mp_upright(2,1) = -mp_upright(1,2) ! and, the inverse of this matrix is: mp_downright(1,1) = mp_upright(1,1) mp_downright(1,2) = mp_upright(2,1) mp_downright(2,1) = mp_upright(1,2) mp_downright(2,2) = mp_upright(2,2) ! Set up the cut and the high- and low-latitude boundaries: ! High "latitude" (in the possibly-rotated coordinate system) limit: mp_using_top = (latitude_0_radians < Pi_over_4).OR. & ! limit at 45 deg. N from projpoint & (.NOT.((latitude_2_radians == Pi_over_2).OR.(latitude_0_radians == Pi_over_2))) ! pole is an arc, except in these two special cases IF (mp_using_top) THEN mp_top_uvec = mp_pole_3_uvec mp_top_applies_uvec = mp_top_uvec mp_top_maxdot = SIN(MIN(1.4835,latitude_0_radians + Pi_over_4)) ! < 1. END IF ! Low "latitude" (in the possibly-rotated coordinate system) limit: mp_using_bottom = .TRUE. mp_bottom_uvec = -mp_pole_3_uvec mp_bottom_applies_uvec = mp_bottom_uvec mp_bottom_maxdot = MAX(0., -SIN(latitude_0_radians - Pi_over_4)) ! Left ("West-longitude") side of the cut: mp_using_left = .TRUE. t_vec = mp_pole_2_uvec - 0.1 * radians_per_degree * mp_pole_1_uvec CALL Make_Uvec(t_vec, mp_left_uvec) ! Notes: Because of wrap-around, it points "right". ! It is 0.2 degrees of "longitude" away ! from the corresponding "right" limit. t_vec = -mp_pole_2_uvec - 0.05 * radians_per_degree * mp_pole_1_uvec CALL Make_Uvec(t_vec, mp_left_applies_uvec) mp_left_maxdot = 0. ! cut is a great-circle arc ! Right ("East-longitude") side of the cut: mp_using_right = .TRUE. t_vec = -mp_pole_2_uvec - 0.1 * radians_per_degree * mp_pole_1_uvec CALL Make_Uvec(t_vec, mp_right_uvec) ! Notes: Because of wrap-around, it points "right". ! It is 0.2 degrees of "longitude" away ! from the corresponding "right" limit. t_vec = mp_pole_2_uvec - 0.05 * radians_per_degree * mp_pole_1_uvec CALL Make_Uvec(t_vec, mp_right_applies_uvec) mp_right_maxdot = 0. ! cut is a great-circle arc ! Complete description of rind by finding corners: CALL Corners (projpoint_uvec) ! Set up the (x,y) system in the projection plane: mp_user_xy = .FALSE. IF (PRESENT(x_projpoint_meters)) THEN mp_user_xy = .TRUE. mp_x_projpoint_meters = x_projpoint_meters ELSE mp_x_projpoint_meters = 0. END IF IF (PRESENT(y_projpoint_meters)) THEN mp_user_xy = .TRUE. mp_y_projpoint_meters = y_projpoint_meters ELSE mp_y_projpoint_meters = 0. END IF IF (PRESENT(y_azimuth_radians)) THEN mp_user_xy = .TRUE. mp_y_azimuth_radians = y_azimuth_radians ELSE mp_y_azimuth_radians = 0. END IF IF (mp_user_xy) THEN ! from "natural" (x,y) system of projection plane to user system: mp_pp_2_xy_matrix(1,1) = COS(mp_y_azimuth_radians) mp_pp_2_xy_matrix(2,2) = mp_pp_2_xy_matrix(1,1) mp_pp_2_xy_matrix(2,1) = SIN(mp_y_azimuth_radians) mp_pp_2_xy_matrix(1,2) = -mp_pp_2_xy_matrix(2,1) mp_pp_2_xy_vector(1) = mp_x_projpoint_meters mp_pp_2_xy_vector(2) = mp_y_projpoint_meters ! from user (x,y) system to "natural" system of projection plane: mp_xy_2_pp_matrix(1,1) = COS(mp_y_azimuth_radians) mp_xy_2_pp_matrix(2,2) = mp_xy_2_pp_matrix(1,1) mp_xy_2_pp_matrix(1,2) = SIN(mp_y_azimuth_radians) mp_xy_2_pp_matrix(2,1) = -mp_xy_2_pp_matrix(1,2) mp_xy_2_pp_vector(1) = -mp_xy_2_pp_matrix(1,1)*mp_x_projpoint_meters & & -mp_xy_2_pp_matrix(1,2)*mp_y_projpoint_meters mp_xy_2_pp_vector(2) = -mp_xy_2_pp_matrix(2,1)*mp_x_projpoint_meters & & -mp_xy_2_pp_matrix(2,2)*mp_y_projpoint_meters END IF END SUBROUTINE Set_Albers_Equal_Area_Conic SUBROUTINE Set_Polyconic (radius_meters, & & projpoint_uvec, cone_pole_uvec, & & x_projpoint_meters, y_projpoint_meters, & & y_azimuth_radians) ! Radius_meters is the radius of the planet, in meters. ! - - - - - - - - - - - - - - - - - - - - - - - - ! A "uvec" is a unit vector from the center of the ! sphere to a surface point in a right-handed Cartesian system ! with the first axis at (0 E, 0N), the second at (90 E, 0 N), ! and the third at the North pole (?? E, 90 N). ! (See notes on uvec's in the module data above.) ! - - - - - - - - - - - - - - - - - - - - - - - - ! Projpoint_uvec is a point on the central meridian ! (the one projected as straight, if cone_pole_uvec is N or S pole). ! At this origin point, +x will be Eastward ! and +y will be Northward in the "native" (x,y) of the ! projection plane (as opposed to any user coordinates). ! Cone_pole_uvec is aligned with the axis of the cone, pointing ! toward the closed end. ! The last 3 arguments are provided in case an (x,y) system ! was already created, with no relation to the map projection: ! X_projpoint_meters is the x value desired at projpoint_uvec. ! Y_projpoint_meters is the y value desired at projpoint_uvec. ! Y_azimuth_radians is the azimuth of the +y axis desired at ! point projpoint_uvec. ! The default values for these last 3 arguments are zero. ! These will give the default definition of ! the 2-D Cartesian coordinates in the (unrolled) projection ! plane, which is (0., 0.) at projpoint_uvec, with +y Northward ! and +x Eastward at that point. IMPLICIT NONE REAL, INTENT(IN) :: radius_meters REAL, DIMENSION(3), INTENT(IN) :: projpoint_uvec, cone_pole_uvec REAL, INTENT(IN), OPTIONAL :: x_projpoint_meters, y_projpoint_meters, & & y_azimuth_radians REAL :: toconepole_azimuth_radians REAL, DIMENSION(3) :: t_vec mp_projection_defined = .TRUE. mp_projection = 'Polyconic' mp_projection_number = 4 mp_radius_meters = radius_meters CALL Make_Uvec(projpoint_uvec, mp_projpoint_uvec) ! just for sure! IF (Arc(mp_projpoint_uvec, cone_pole_uvec) > 1.396) THEN WRITE (*,"(' ERROR: Distance from projection point to pole of cone exceeds 80 degrees.')") CALL Traceback END IF ! find the 3 orthogonal uvec's used for all conic projections: CALL Make_Uvec (cone_pole_uvec, mp_pole_3_uvec) ! just for insurance!!! CALL Cross (mp_pole_3_uvec, mp_projpoint_uvec, t_vec) CALL Make_Uvec (t_vec, mp_pole_2_uvec) CALL Cross (mp_pole_2_uvec, mp_pole_3_uvec, mp_pole_1_uvec) ! Note: Pole_1 is to S of projpoint, if cone_pole is N pole. ! Parameters corresponding to pages 128-129 of John P. Snyder (1983), ! US Geological Survey Bulletin 1532: mp_proj4_lat0_radians = Pi_over_2 - Arc(mp_projpoint_uvec, cone_pole_uvec) ! Note: "Latitude" is now with respect to conic_pole_uvec as "North pole"! mp_proj4_lon0_radians = 0. ! because of rotated, custom definition of "longitude" in pole1/2/3 system. ! Predefine a 2x2 matrix for the "twist" from ! [ up = toward the pole of the cone ] TO [ up = North ], ! the "natural" or default (x,y) in the projection plane: toconepole_azimuth_radians = Relative_Compass(mp_projpoint_uvec, mp_pole_3_uvec) mp_upright(1,1) = COS(toconepole_azimuth_radians) mp_upright(2,2) = mp_upright(1,1) mp_upright(1,2) = SIN(toconepole_azimuth_radians) mp_upright(2,1) = -mp_upright(1,2) ! and, the inverse of this matrix is: mp_downright(1,1) = mp_upright(1,1) mp_downright(1,2) = mp_upright(2,1) mp_downright(2,1) = mp_upright(1,2) mp_downright(2,2) = mp_upright(2,2) ! Set up the cut and the high- and low-latitude boundaries: ! High "latitude" (in the possibly-rotated coordinate system) limit: mp_using_top = .FALSE. ! Low "latitude" (in the possibly-rotated coordinate system) limit: ! Since polyconic equation are ill-conditioned for relative latitudes ! close to 0., bottom is tilted forward to cut out the last few degrees. mp_using_bottom = .TRUE. t_vec = -mp_pole_3_uvec + 0.036 * mp_pole_1_uvec ! 2 degree tilt CALL Make_Uvec (t_vec, mp_bottom_uvec) mp_bottom_applies_uvec = mp_bottom_uvec mp_bottom_maxdot = 0. ! Left ("West-longitude") side of the cut: ! Since inverse formulas for Polyconic rejection do not converge ! at or beyond 90 degrees longitude from the central meridian ! (the one with the projection point), both side boundaries ! are placed 80 degrees away. Actually, the distortion is ! already quite bad at this distance! mp_using_left = .TRUE. t_vec = -SIN(80.01*radians_per_degree)*mp_pole_1_uvec & & -COS(80.01*radians_per_degree)*mp_pole_2_uvec CALL Make_Uvec(t_vec, mp_left_uvec) mp_left_applies_uvec = mp_left_uvec mp_left_maxdot = 0. ! cut is a great-circle arc ! Right ("East-longitude") side of the cut: mp_using_right = .TRUE. t_vec = -SIN(80.01*radians_per_degree)*mp_pole_1_uvec & & +COS(80.01*radians_per_degree)*mp_pole_2_uvec CALL Make_Uvec(t_vec, mp_right_uvec) mp_right_applies_uvec = mp_right_uvec mp_right_maxdot = 0. ! cut is a great-circle arc ! Complete description of rind by finding corners: CALL Corners (projpoint_uvec) ! Set up the (x,y) system in the projection plane: mp_user_xy = .FALSE. IF (PRESENT(x_projpoint_meters)) THEN mp_user_xy = .TRUE. mp_x_projpoint_meters = x_projpoint_meters ELSE mp_x_projpoint_meters = 0. END IF IF (PRESENT(y_projpoint_meters)) THEN mp_user_xy = .TRUE. mp_y_projpoint_meters = y_projpoint_meters ELSE mp_y_projpoint_meters = 0. END IF IF (PRESENT(y_azimuth_radians)) THEN mp_user_xy = .TRUE. mp_y_azimuth_radians = y_azimuth_radians ELSE mp_y_azimuth_radians = 0. END IF IF (mp_user_xy) THEN ! from "natural" (x,y) system of projection plane to user system: mp_pp_2_xy_matrix(1,1) = COS(mp_y_azimuth_radians) mp_pp_2_xy_matrix(2,2) = mp_pp_2_xy_matrix(1,1) mp_pp_2_xy_matrix(2,1) = SIN(mp_y_azimuth_radians) mp_pp_2_xy_matrix(1,2) = -mp_pp_2_xy_matrix(2,1) mp_pp_2_xy_vector(1) = mp_x_projpoint_meters mp_pp_2_xy_vector(2) = mp_y_projpoint_meters ! from user (x,y) system to "natural" system of projection plane: mp_xy_2_pp_matrix(1,1) = COS(mp_y_azimuth_radians) mp_xy_2_pp_matrix(2,2) = mp_xy_2_pp_matrix(1,1) mp_xy_2_pp_matrix(1,2) = SIN(mp_y_azimuth_radians) mp_xy_2_pp_matrix(2,1) = -mp_xy_2_pp_matrix(1,2) mp_xy_2_pp_vector(1) = -mp_xy_2_pp_matrix(1,1)*mp_x_projpoint_meters & & -mp_xy_2_pp_matrix(1,2)*mp_y_projpoint_meters mp_xy_2_pp_vector(2) = -mp_xy_2_pp_matrix(2,1)*mp_x_projpoint_meters & & -mp_xy_2_pp_matrix(2,2)*mp_y_projpoint_meters END IF END SUBROUTINE Set_Polyconic SUBROUTINE Set_Geometric_Conic & & (radius_meters, & & projpoint_uvec, cone_pole_uvec, & & x_projpoint_meters, y_projpoint_meters, & & y_azimuth_radians) ! This conic projection is purely geometric: ! * A plane is rolled into a cone. ! * The cone is tangent to the sphere along one small circle, ! passing through location "projpoint_uvec". The axis of ! the cone passes through "cone_pole_uvec". ! * Surface points are projected to the cone along extended radii. ! * The cone is cut opposite to projpoint_uvec, and unrolled. ! NOTE that this projection is rarely used, and has no standard ! name. It is included here for continuity with my older ! software: when it specifies a "conic projection", this ! is what is meant. !============================================================= ! Radius_meters is the radius of the planet, in meters. ! - - - - - - - - - - - - - - - - - - - - - - - - ! A "uvec" is a unit vector from the center of the ! sphere to a surface point in a right-handed Cartesian system ! with the first axis at (0 E, 0N), the second at (90 E, 0 N), ! and the third at the North pole (?? E, 90 N). ! (See notes on uvec's in the module data above.) ! - - - - - - - - - - - - - - - - - - - - - - - - ! Projpoint_uvec is a point on the small circle of tangency ! between the sphere and the projection plane (while rolled ! into a cone). At this origin point, +x will be Eastward ! and +y will be Northward in the "native" (x,y) of the ! projection plane (as opposed to any user coordinates). ! Cone_pole_uvec is aligned with the axis of the cone, pointing ! toward the closed end. ! The last 3 arguments are provided in case an (x,y) system ! was already created, with no relation to the map projection: ! X_projpoint_meters is the x value desired at projpoint_uvec. ! Y_projpoint_meters is the y value desired at projpoint_uvec. ! Y_azimuth_radians is the azimuth of the +y axis desired at ! point projpoint_uvec. ! The default values for these last 3 arguments are zero. ! These will give the default definition of ! the 2-D Cartesian coordinates in the (unrolled) projection ! plane, which is (0., 0.) at projpoint_uvec, with +y Northward ! and +x Eastward at that point. IMPLICIT NONE REAL, INTENT(IN) :: radius_meters REAL, DIMENSION(3), INTENT(IN) :: projpoint_uvec, cone_pole_uvec REAL, INTENT(IN), OPTIONAL :: x_projpoint_meters, y_projpoint_meters, & & y_azimuth_radians REAL :: Pi_over_4, toconepole_azimuth_radians REAL, DIMENSION(3) :: t_vec mp_projection_defined = .TRUE. mp_projection = 'Geometric Conic' mp_projection_number = 5 mp_radius_meters = radius_meters CALL Make_Uvec(projpoint_uvec, mp_projpoint_uvec) ! just for sure! IF (Arc(mp_projpoint_uvec, cone_pole_uvec) > 1.396) THEN WRITE (*,"(' ERROR: Distance from projection point to pole of cone exceeds 80 degrees.')") CALL Traceback END IF ! find the 3 orthogonal uvec's used for all conic projections: CALL Make_Uvec (cone_pole_uvec, mp_pole_3_uvec) ! just for insurance!!! CALL Cross (mp_pole_3_uvec, mp_projpoint_uvec, t_vec) CALL Make_Uvec (t_vec, mp_pole_2_uvec) CALL Cross (mp_pole_2_uvec, mp_pole_3_uvec, mp_pole_1_uvec) ! Note: Pole_1 is to S of projpoint, if cone_pole is N pole. Pi_over_4 = Pi_over_2 / 2. mp_proj5_lambda0 = 0. ! by definition of rotated coordinate system mp_proj5_lat0_radians = Pi_over_2 - Arc(mp_projpoint_uvec, mp_pole_3_uvec) mp_proj5_n = SIN(mp_proj5_lat0_radians) ! ratio of "longitude" angles in unrolled cone to sphere mp_proj5_rtan = mp_radius_meters * TAN(Arc(mp_projpoint_uvec, mp_pole_3_uvec)) mp_proj5_ypole = mp_proj5_rtan ! because tangent "latitude" parallel passes through ! projpoint_uvec. If it didn't, formula would be (in Fortran77): ! YPOLE=RTAN-RADIUS*TANDEG(Y0NLAT-CPNLAT) ! Predefine a 2x2 matrix for the "twist" from ! [ up = toward the pole of the cone ] TO [ up = North ], ! the "natural" or default (x,y) in the projection plane: toconepole_azimuth_radians = Relative_Compass(mp_projpoint_uvec, mp_pole_3_uvec) mp_upright(1,1) = COS(toconepole_azimuth_radians) mp_upright(2,2) = mp_upright(1,1) mp_upright(1,2) = SIN(toconepole_azimuth_radians) mp_upright(2,1) = -mp_upright(1,2) ! and, the inverse of this matrix is: mp_downright(1,1) = mp_upright(1,1) mp_downright(1,2) = mp_upright(2,1) mp_downright(2,1) = mp_upright(1,2) mp_downright(2,2) = mp_upright(2,2) ! Set up the cut and the high- and low-latitude boundaries: ! High "latitude" (in the possibly-rotated coordinate system) limit: mp_using_top = (mp_proj5_lat0_radians < Pi_over_4) ! limit at 45 deg. N from projpoint? IF (mp_using_top) THEN mp_top_uvec = mp_pole_3_uvec mp_top_applies_uvec = mp_top_uvec mp_top_maxdot = SIN(mp_proj5_lat0_radians + Pi_over_4) ! < 1. END IF ! Low "latitude" (in the possibly-rotated coordinate system) limit: mp_using_bottom = .TRUE. mp_bottom_uvec = -mp_pole_3_uvec mp_bottom_applies_uvec = mp_bottom_uvec mp_bottom_maxdot = MAX(0., -SIN(mp_proj5_lat0_radians - Pi_over_4)) ! Left ("West-longitude") side of the cut: mp_using_left = .TRUE. t_vec = mp_pole_2_uvec - 0.1 * radians_per_degree * mp_pole_1_uvec CALL Make_Uvec(t_vec, mp_left_uvec) ! Notes: Because of wrap-around, it points "right". ! It is 0.2 degrees of "longitude" away ! from the corresponding "right" limit. t_vec = -mp_pole_2_uvec - 0.05 * radians_per_degree * mp_pole_1_uvec CALL Make_Uvec(t_vec, mp_left_applies_uvec) mp_left_maxdot = 0. ! cut is a great-circle arc ! Right ("East-longitude") side of the cut: mp_using_right = .TRUE. t_vec = -mp_pole_2_uvec - 0.1 * radians_per_degree * mp_pole_1_uvec CALL Make_Uvec(t_vec, mp_right_uvec) ! Notes: Because of wrap-around, it points "right". ! It is 0.2 degrees of "longitude" away ! from the corresponding "right" limit. t_vec = mp_pole_2_uvec - 0.05 * radians_per_degree * mp_pole_1_uvec CALL Make_Uvec(t_vec, mp_right_applies_uvec) mp_right_maxdot = 0. ! cut is a great-circle arc ! Complete description of rind by finding corners: CALL Corners (projpoint_uvec) ! Set up the (x,y) system in the projection plane: mp_user_xy = .FALSE. IF (PRESENT(x_projpoint_meters)) THEN mp_user_xy = .TRUE. mp_x_projpoint_meters = x_projpoint_meters ELSE mp_x_projpoint_meters = 0. END IF IF (PRESENT(y_projpoint_meters)) THEN mp_user_xy = .TRUE. mp_y_projpoint_meters = y_projpoint_meters ELSE mp_y_projpoint_meters = 0. END IF IF (PRESENT(y_azimuth_radians)) THEN mp_user_xy = .TRUE. mp_y_azimuth_radians = y_azimuth_radians ELSE mp_y_azimuth_radians = 0. END IF IF (mp_user_xy) THEN ! from "natural" (x,y) system of projection plane to user system: mp_pp_2_xy_matrix(1,1) = COS(mp_y_azimuth_radians) mp_pp_2_xy_matrix(2,2) = mp_pp_2_xy_matrix(1,1) mp_pp_2_xy_matrix(2,1) = SIN(mp_y_azimuth_radians) mp_pp_2_xy_matrix(1,2) = -mp_pp_2_xy_matrix(2,1) mp_pp_2_xy_vector(1) = mp_x_projpoint_meters mp_pp_2_xy_vector(2) = mp_y_projpoint_meters ! from user (x,y) system to "natural" system of projection plane: mp_xy_2_pp_matrix(1,1) = COS(mp_y_azimuth_radians) mp_xy_2_pp_matrix(2,2) = mp_xy_2_pp_matrix(1,1) mp_xy_2_pp_matrix(1,2) = SIN(mp_y_azimuth_radians) mp_xy_2_pp_matrix(2,1) = -mp_xy_2_pp_matrix(1,2) mp_xy_2_pp_vector(1) = -mp_xy_2_pp_matrix(1,1)*mp_x_projpoint_meters & & -mp_xy_2_pp_matrix(1,2)*mp_y_projpoint_meters mp_xy_2_pp_vector(2) = -mp_xy_2_pp_matrix(2,1)*mp_x_projpoint_meters & & -mp_xy_2_pp_matrix(2,2)*mp_y_projpoint_meters END IF END SUBROUTINE Set_Geometric_Conic SUBROUTINE Set_Stereographic (radius_meters, & & projpoint_uvec, & & x_projpoint_meters, y_projpoint_meters, & & y_azimuth_radians) ! Radius_meters is the radius of the planet. ! - - - - - - - - - - - - - - - - - - - - - - - - ! A "uvec" is a unit vector from the center of the ! sphere to a surface point in a right-handed Cartesian system ! with the first axis at (0 E, 0N), the second at (90 E, 0 N), ! and the third at the North pole (?? E, 90 N). ! (See notes on uvec's in the module data above.) ! - - - - - - - - - - - - - - - - - - - - - - - - ! Projpoint_uvec is the point of tangency ! between the sphere and the projection plane. ! Not necessarily on the equator!!! ! The last 3 arguments are provided in case an (x,y) system ! was already created, with no relation to the map projection: ! X_projpoint_meters is the x value desired at projpoint_uvec. ! Y_projpoint_meters is the y value desired at projpoint_uvec. ! Y_azimuth_radians is the azimuth of the +y axis desired at ! point projpoint_uvec. ! The default values for these last 3 arguments are zero. ! These will give the default definition of ! the 2-D Cartesian coordinates in the (unrolled) projection ! plane, which is (0., 0.) at projpoint_uvec, with +y Northward ! and +x Eastward at that point. IMPLICIT NONE REAL, INTENT(IN) :: radius_meters REAL, DIMENSION(3), INTENT(IN) :: projpoint_uvec REAL, INTENT(IN), OPTIONAL :: x_projpoint_meters, y_projpoint_meters, & & y_azimuth_radians REAL, DIMENSION(3) :: east_uvec, north_uvec mp_projection_defined = .TRUE. mp_projection = 'Stereographic' mp_projection_number = 6 mp_radius_meters = radius_meters CALL Make_Uvec(projpoint_uvec, mp_projpoint_uvec) ! just for sure! ! find the 3 orthogonal uvec's used for Stereographic projection: mp_pole_1_uvec = mp_projpoint_uvec CALL NorthEast_Convention (mp_pole_1_uvec, north_uvec, east_uvec) mp_pole_2_uvec = east_uvec mp_pole_3_uvec = north_uvec ! Set up the single great-circle boundary of the projected area: mp_using_top = .FALSE. mp_using_left = .FALSE. mp_using_right = .FALSE. ! NOTE: With this projection, whole planet CAN be shown, ! but far hemisphere is terribly distorted! Here, I allow ! only the near hemisphere. mp_using_bottom = .TRUE. mp_bottom_uvec = -mp_pole_1_uvec mp_bottom_applies_uvec = mp_bottom_uvec mp_bottom_maxdot = 0.001 ! Complete description of rind by finding corners: CALL Corners (projpoint_uvec) ! Set up the (x,y) system in the projection plane: mp_user_xy = .FALSE. IF (PRESENT(x_projpoint_meters)) THEN mp_user_xy = .TRUE. mp_x_projpoint_meters = x_projpoint_meters ELSE mp_x_projpoint_meters = 0. END IF IF (PRESENT(y_projpoint_meters)) THEN mp_user_xy = .TRUE. mp_y_projpoint_meters = y_projpoint_meters ELSE mp_y_projpoint_meters = 0. END IF IF (PRESENT(y_azimuth_radians)) THEN mp_user_xy = .TRUE. mp_y_azimuth_radians = y_azimuth_radians ELSE mp_y_azimuth_radians = 0. END IF IF (mp_user_xy) THEN ! from "natural" (x,y) system of projection plane to user system: mp_pp_2_xy_matrix(1,1) = COS(mp_y_azimuth_radians) mp_pp_2_xy_matrix(2,2) = mp_pp_2_xy_matrix(1,1) mp_pp_2_xy_matrix(2,1) = SIN(mp_y_azimuth_radians) mp_pp_2_xy_matrix(1,2) = -mp_pp_2_xy_matrix(2,1) mp_pp_2_xy_vector(1) = mp_x_projpoint_meters mp_pp_2_xy_vector(2) = mp_y_projpoint_meters ! from user (x,y) system to "natural" system of projection plane: mp_xy_2_pp_matrix(1,1) = COS(mp_y_azimuth_radians) mp_xy_2_pp_matrix(2,2) = mp_xy_2_pp_matrix(1,1) mp_xy_2_pp_matrix(1,2) = SIN(mp_y_azimuth_radians) mp_xy_2_pp_matrix(2,1) = -mp_xy_2_pp_matrix(1,2) mp_xy_2_pp_vector(1) = -mp_xy_2_pp_matrix(1,1)*mp_x_projpoint_meters & & -mp_xy_2_pp_matrix(1,2)*mp_y_projpoint_meters mp_xy_2_pp_vector(2) = -mp_xy_2_pp_matrix(2,1)*mp_x_projpoint_meters & & -mp_xy_2_pp_matrix(2,2)*mp_y_projpoint_meters END IF END SUBROUTINE Set_Stereographic SUBROUTINE Set_Lambert_Azimuthal_EqualArea ( & & radius_meters, & & projpoint_uvec, & & x_projpoint_meters, y_projpoint_meters, & & y_azimuth_radians) ! Radius_meters is the radius of the planet. ! - - - - - - - - - - - - - - - - - - - - - - - - ! A "uvec" is a unit vector from the center of the ! sphere to a surface point in a right-handed Cartesian system ! with the first axis at (0 E, 0N), the second at (90 E, 0 N), ! and the third at the North pole (?? E, 90 N). ! (See notes on uvec's in the module data above.) ! - - - - - - - - - - - - - - - - - - - - - - - - ! Projpoint_uvec is the point of tangency ! between the sphere and the projection plane. ! Not necessarily on the equator!!! ! The last 3 arguments are provided in case an (x,y) system ! was already created, with no relation to the map projection: ! X_projpoint_meters is the x value desired at projpoint_uvec. ! Y_projpoint_meters is the y value desired at projpoint_uvec. ! Y_azimuth_radians is the azimuth of the +y axis desired at ! point projpoint_uvec. ! The default values for these last 3 arguments are zero. ! These will give the default definition of ! the 2-D Cartesian coordinates in the (unrolled) projection ! plane, which is (0., 0.) at projpoint_uvec, with +y Northward ! and +x Eastward at that point. IMPLICIT NONE REAL, INTENT(IN) :: radius_meters REAL, DIMENSION(3), INTENT(IN) :: projpoint_uvec REAL, INTENT(IN), OPTIONAL :: x_projpoint_meters, y_projpoint_meters, & & y_azimuth_radians REAL, DIMENSION(3) :: east_uvec, north_uvec mp_projection_defined = .TRUE. mp_projection = 'Lambert Azimuthal Equal-Area' mp_projection_number = 7 mp_radius_meters = radius_meters CALL Make_Uvec(projpoint_uvec, mp_projpoint_uvec) ! just for sure! ! find the 3 orthogonal uvec's used for Stereographic projection: mp_pole_1_uvec = mp_projpoint_uvec CALL NorthEast_Convention (mp_pole_1_uvec, north_uvec, east_uvec) mp_pole_2_uvec = east_uvec mp_pole_3_uvec = north_uvec ! Set up the single great-circle boundary of the projected area: mp_using_top = .FALSE. mp_using_left = .FALSE. mp_using_right = .FALSE. mp_using_bottom = .TRUE. mp_bottom_uvec = -mp_pole_1_uvec mp_bottom_applies_uvec = mp_bottom_uvec mp_bottom_maxdot = 0.001 ! Complete description of rind by finding corners: CALL Corners (projpoint_uvec) ! Set up the (x,y) system in the projection plane: mp_user_xy = .FALSE. IF (PRESENT(x_projpoint_meters)) THEN mp_user_xy = .TRUE. mp_x_projpoint_meters = x_projpoint_meters ELSE mp_x_projpoint_meters = 0. END IF IF (PRESENT(y_projpoint_meters)) THEN mp_user_xy = .TRUE. mp_y_projpoint_meters = y_projpoint_meters ELSE mp_y_projpoint_meters = 0. END IF IF (PRESENT(y_azimuth_radians)) THEN mp_user_xy = .TRUE. mp_y_azimuth_radians = y_azimuth_radians ELSE mp_y_azimuth_radians = 0. END IF IF (mp_user_xy) THEN ! from "natural" (x,y) system of projection plane to user system: mp_pp_2_xy_matrix(1,1) = COS(mp_y_azimuth_radians) mp_pp_2_xy_matrix(2,2) = mp_pp_2_xy_matrix(1,1) mp_pp_2_xy_matrix(2,1) = SIN(mp_y_azimuth_radians) mp_pp_2_xy_matrix(1,2) = -mp_pp_2_xy_matrix(2,1) mp_pp_2_xy_vector(1) = mp_x_projpoint_meters mp_pp_2_xy_vector(2) = mp_y_projpoint_meters ! from user (x,y) system to "natural" system of projection plane: mp_xy_2_pp_matrix(1,1) = COS(mp_y_azimuth_radians) mp_xy_2_pp_matrix(2,2) = mp_xy_2_pp_matrix(1,1) mp_xy_2_pp_matrix(1,2) = SIN(mp_y_azimuth_radians) mp_xy_2_pp_matrix(2,1) = -mp_xy_2_pp_matrix(1,2) mp_xy_2_pp_vector(1) = -mp_xy_2_pp_matrix(1,1)*mp_x_projpoint_meters & & -mp_xy_2_pp_matrix(1,2)*mp_y_projpoint_meters mp_xy_2_pp_vector(2) = -mp_xy_2_pp_matrix(2,1)*mp_x_projpoint_meters & & -mp_xy_2_pp_matrix(2,2)*mp_y_projpoint_meters END IF END SUBROUTINE Set_Lambert_Azimuthal_EqualArea SUBROUTINE Set_Azimuthal_Equidistant ( & & radius_meters, & & projpoint_uvec, & & x_projpoint_meters, y_projpoint_meters, & & y_azimuth_radians) ! Radius_meters is the radius of the planet. ! - - - - - - - - - - - - - - - - - - - - - - - - ! A "uvec" is a unit vector from the center of the ! sphere to a surface point in a right-handed Cartesian system ! with the first axis at (0 E, 0N), the second at (90 E, 0 N), ! and the third at the North pole (?? E, 90 N). ! (See notes on uvec's in the module data above.) ! - - - - - - - - - - - - - - - - - - - - - - - - ! Projpoint_uvec is the point of tangency ! between the sphere and the projection plane. ! Not necessarily on the equator!!! ! The last 3 arguments are provided in case an (x,y) system ! was already created, with no relation to the map projection: ! X_projpoint_meters is the x value desired at projpoint_uvec. ! Y_projpoint_meters is the y value desired at projpoint_uvec. ! Y_azimuth_radians is the azimuth of the +y axis desired at ! point projpoint_uvec. ! The default values for these last 3 arguments are zero. ! These will give the default definition of ! the 2-D Cartesian coordinates in the (unrolled) projection ! plane, which is (0., 0.) at projpoint_uvec, with +y Northward ! and +x Eastward at that point. IMPLICIT NONE REAL, INTENT(IN) :: radius_meters REAL, DIMENSION(3), INTENT(IN) :: projpoint_uvec REAL, INTENT(IN), OPTIONAL :: x_projpoint_meters, y_projpoint_meters, & & y_azimuth_radians REAL, DIMENSION(3) :: east_uvec, north_uvec mp_projection_defined = .TRUE. mp_projection = 'Azimuthal Equidistant' mp_projection_number = 8 mp_radius_meters = radius_meters CALL Make_Uvec(projpoint_uvec, mp_projpoint_uvec) ! just for sure! ! find the 3 orthogonal uvec's used for Stereographic projection: mp_pole_1_uvec = mp_projpoint_uvec CALL NorthEast_Convention (mp_pole_1_uvec, north_uvec, east_uvec) mp_pole_2_uvec = east_uvec mp_pole_3_uvec = north_uvec ! Set up the single great-circle boundary of the projected area: mp_using_top = .FALSE. mp_using_left = .FALSE. mp_using_right = .FALSE. mp_using_bottom = .TRUE. ! NOTE: With this projection, whole planet CAN be shown, but ! back hemisphere is terribly distorted! Thus, I only allow ! near hemisphere (as with Stereographic projection). mp_bottom_uvec = -mp_pole_1_uvec mp_bottom_applies_uvec = mp_bottom_uvec mp_bottom_maxdot = 0.001 ! Complete description of rind by finding corners: CALL Corners (projpoint_uvec) ! Set up the (x,y) system in the projection plane: mp_user_xy = .FALSE. IF (PRESENT(x_projpoint_meters)) THEN mp_user_xy = .TRUE. mp_x_projpoint_meters = x_projpoint_meters ELSE mp_x_projpoint_meters = 0. END IF IF (PRESENT(y_projpoint_meters)) THEN mp_user_xy = .TRUE. mp_y_projpoint_meters = y_projpoint_meters ELSE mp_y_projpoint_meters = 0. END IF IF (PRESENT(y_azimuth_radians)) THEN mp_user_xy = .TRUE. mp_y_azimuth_radians = y_azimuth_radians ELSE mp_y_azimuth_radians = 0. END IF IF (mp_user_xy) THEN ! from "natural" (x,y) system of projection plane to user system: mp_pp_2_xy_matrix(1,1) = COS(mp_y_azimuth_radians) mp_pp_2_xy_matrix(2,2) = mp_pp_2_xy_matrix(1,1) mp_pp_2_xy_matrix(2,1) = SIN(mp_y_azimuth_radians) mp_pp_2_xy_matrix(1,2) = -mp_pp_2_xy_matrix(2,1) mp_pp_2_xy_vector(1) = mp_x_projpoint_meters mp_pp_2_xy_vector(2) = mp_y_projpoint_meters ! from user (x,y) system to "natural" system of projection plane: mp_xy_2_pp_matrix(1,1) = COS(mp_y_azimuth_radians) mp_xy_2_pp_matrix(2,2) = mp_xy_2_pp_matrix(1,1) mp_xy_2_pp_matrix(1,2) = SIN(mp_y_azimuth_radians) mp_xy_2_pp_matrix(2,1) = -mp_xy_2_pp_matrix(1,2) mp_xy_2_pp_vector(1) = -mp_xy_2_pp_matrix(1,1)*mp_x_projpoint_meters & & -mp_xy_2_pp_matrix(1,2)*mp_y_projpoint_meters mp_xy_2_pp_vector(2) = -mp_xy_2_pp_matrix(2,1)*mp_x_projpoint_meters & & -mp_xy_2_pp_matrix(2,2)*mp_y_projpoint_meters END IF END SUBROUTINE Set_Azimuthal_Equidistant SUBROUTINE Set_Orthographic ( radius_meters, & & projpoint_uvec, & & x_projpoint_meters, y_projpoint_meters, & & y_azimuth_radians) ! Radius_meters is the radius of the planet. ! - - - - - - - - - - - - - - - - - - - - - - - - ! A "uvec" is a unit vector from the center of the ! sphere to a surface point in a right-handed Cartesian system ! with the first axis at (0 E, 0N), the second at (90 E, 0 N), ! and the third at the North pole (?? E, 90 N). ! (See notes on uvec's in the module data above.) ! - - - - - - - - - - - - - - - - - - - - - - - - ! Projpoint_uvec is the point of tangency ! between the sphere and the projection plane. ! Not necessarily on the equator!!! ! The last 3 arguments are provided in case an (x,y) system ! was already created, with no relation to the map projection: ! X_projpoint_meters is the x value desired at projpoint_uvec. ! Y_projpoint_meters is the y value desired at projpoint_uvec. ! Y_azimuth_radians is the azimuth of the +y axis desired at ! point projpoint_uvec. ! The default values for these last 3 arguments are zero. ! These will give the default definition of ! the 2-D Cartesian coordinates in the (unrolled) projection ! plane, which is (0., 0.) at projpoint_uvec, with +y Northward ! and +x Eastward at that point. IMPLICIT NONE REAL, INTENT(IN) :: radius_meters REAL, DIMENSION(3), INTENT(IN) :: projpoint_uvec REAL, INTENT(IN), OPTIONAL :: x_projpoint_meters, y_projpoint_meters, & & y_azimuth_radians REAL, DIMENSION(3) :: east_uvec, north_uvec mp_projection_defined = .TRUE. mp_projection = 'Orthographic' mp_projection_number = 9 mp_radius_meters = radius_meters CALL Make_Uvec(projpoint_uvec, mp_projpoint_uvec) ! just for sure! ! find the 3 orthogonal uvec's used for Stereographic projection: mp_pole_1_uvec = mp_projpoint_uvec CALL NorthEast_Convention (mp_pole_1_uvec, north_uvec, east_uvec) mp_pole_2_uvec = east_uvec mp_pole_3_uvec = north_uvec ! Set up the single great-circle boundary of the projected area: mp_using_top = .FALSE. mp_using_left = .FALSE. mp_using_right = .FALSE. mp_using_bottom = .TRUE. mp_bottom_uvec = -mp_pole_1_uvec mp_bottom_applies_uvec = mp_bottom_uvec mp_bottom_maxdot = 0.001 ! Complete description of rind by finding corners: CALL Corners (projpoint_uvec) ! Set up the (x,y) system in the projection plane: mp_user_xy = .FALSE. IF (PRESENT(x_projpoint_meters)) THEN mp_user_xy = .TRUE. mp_x_projpoint_meters = x_projpoint_meters ELSE mp_x_projpoint_meters = 0. END IF IF (PRESENT(y_projpoint_meters)) THEN mp_user_xy = .TRUE. mp_y_projpoint_meters = y_projpoint_meters ELSE mp_y_projpoint_meters = 0. END IF IF (PRESENT(y_azimuth_radians)) THEN mp_user_xy = .TRUE. mp_y_azimuth_radians = y_azimuth_radians ELSE mp_y_azimuth_radians = 0. END IF IF (mp_user_xy) THEN ! from "natural" (x,y) system of projection plane to user system: mp_pp_2_xy_matrix(1,1) = COS(mp_y_azimuth_radians) mp_pp_2_xy_matrix(2,2) = mp_pp_2_xy_matrix(1,1) mp_pp_2_xy_matrix(2,1) = SIN(mp_y_azimuth_radians) mp_pp_2_xy_matrix(1,2) = -mp_pp_2_xy_matrix(2,1) mp_pp_2_xy_vector(1) = mp_x_projpoint_meters mp_pp_2_xy_vector(2) = mp_y_projpoint_meters ! from user (x,y) system to "natural" system of projection plane: mp_xy_2_pp_matrix(1,1) = COS(mp_y_azimuth_radians) mp_xy_2_pp_matrix(2,2) = mp_xy_2_pp_matrix(1,1) mp_xy_2_pp_matrix(1,2) = SIN(mp_y_azimuth_radians) mp_xy_2_pp_matrix(2,1) = -mp_xy_2_pp_matrix(1,2) mp_xy_2_pp_vector(1) = -mp_xy_2_pp_matrix(1,1)*mp_x_projpoint_meters & & -mp_xy_2_pp_matrix(1,2)*mp_y_projpoint_meters mp_xy_2_pp_vector(2) = -mp_xy_2_pp_matrix(2,1)*mp_x_projpoint_meters & & -mp_xy_2_pp_matrix(2,2)*mp_y_projpoint_meters END IF END SUBROUTINE Set_Orthographic SUBROUTINE Set_Gnomonic ( radius_meters, & & projpoint_uvec, & & x_projpoint_meters, y_projpoint_meters, & & y_azimuth_radians) ! Radius_meters is the radius of the planet. ! - - - - - - - - - - - - - - - - - - - - - - - - ! A "uvec" is a unit vector from the center of the ! sphere to a surface point in a right-handed Cartesian system ! with the first axis at (0 E, 0N), the second at (90 E, 0 N), ! and the third at the North pole (?? E, 90 N). ! (See notes on uvec's in the module data above.) ! - - - - - - - - - - - - - - - - - - - - - - - - ! Projpoint_uvec is the point of tangency ! between the sphere and the projection plane. ! Not necessarily on the equator!!! ! The last 3 arguments are provided in case an (x,y) system ! was already created, with no relation to the map projection: ! X_projpoint_meters is the x value desired at projpoint_uvec. ! Y_projpoint_meters is the y value desired at projpoint_uvec. ! Y_azimuth_radians is the azimuth of the +y axis desired at ! point projpoint_uvec. ! The default values for these last 3 arguments are zero. ! These will give the default definition of ! the 2-D Cartesian coordinates in the (unrolled) projection ! plane, which is (0., 0.) at projpoint_uvec, with +y Northward ! and +x Eastward at that point. IMPLICIT NONE REAL, INTENT(IN) :: radius_meters REAL, DIMENSION(3), INTENT(IN) :: projpoint_uvec REAL, INTENT(IN), OPTIONAL :: x_projpoint_meters, y_projpoint_meters, & & y_azimuth_radians REAL, DIMENSION(3) :: east_uvec, north_uvec, t_vec mp_projection_defined = .TRUE. mp_projection = 'Gnomonic' mp_projection_number = 10 mp_radius_meters = radius_meters CALL Make_Uvec(projpoint_uvec, mp_projpoint_uvec) ! just for sure! ! find the 3 orthogonal uvec's used for Stereographic projection: mp_pole_1_uvec = mp_projpoint_uvec CALL NorthEast_Convention (mp_pole_1_uvec, north_uvec, east_uvec) mp_pole_2_uvec = east_uvec mp_pole_3_uvec = north_uvec ! Set up the single great-circle boundary of the projected area: ! The Gnomonic projection has a dangerous singularity at ! Pi_over_2 radians from projpoint_uvec. If paths are allowed ! to even close to this, then control points 1 and 2 of ! Bezier curves could hit or cross over the singularity. ! Thus, we use a combination of all 4 wall to prevent paths ! from going more than about 70 to 80 degrees from projpoint_uvec. ! Since Gnomonic converts all great circles to straight lines, ! the boundaries of the projected map will be straight ! (and, in fact, will form a square). mp_using_top = .TRUE. t_vec = -SIN(70.*radians_per_degree) * mp_pole_1_uvec & & +COS(70.*radians_per_degree) * mp_pole_3_uvec CALL Make_Uvec (t_vec, mp_top_uvec) mp_top_applies_uvec = mp_top_uvec mp_top_maxdot = 0.0 mp_using_left = .TRUE. t_vec = -SIN(70.*radians_per_degree) * mp_pole_1_uvec & & -COS(70.*radians_per_degree) * mp_pole_2_uvec CALL Make_Uvec (t_vec, mp_left_uvec) mp_left_applies_uvec = mp_left_uvec mp_left_maxdot = 0.0 mp_using_right = .TRUE. t_vec = -SIN(70.*radians_per_degree) * mp_pole_1_uvec & & +COS(70.*radians_per_degree) * mp_pole_2_uvec CALL Make_Uvec (t_vec, mp_right_uvec) mp_right_applies_uvec = mp_right_uvec mp_right_maxdot = 0.0 mp_using_bottom = .TRUE. t_vec = -SIN(70.*radians_per_degree) * mp_pole_1_uvec & & -COS(70.*radians_per_degree) * mp_pole_3_uvec CALL Make_Uvec (t_vec, mp_bottom_uvec) mp_bottom_applies_uvec = mp_bottom_uvec mp_bottom_maxdot = 0.0 ! Complete description of rind by finding corners: CALL Corners (projpoint_uvec) ! Set up the (x,y) system in the projection plane: mp_user_xy = .FALSE. IF (PRESENT(x_projpoint_meters)) THEN mp_user_xy = .TRUE. mp_x_projpoint_meters = x_projpoint_meters ELSE mp_x_projpoint_meters = 0. END IF IF (PRESENT(y_projpoint_meters)) THEN mp_user_xy = .TRUE. mp_y_projpoint_meters = y_projpoint_meters ELSE mp_y_projpoint_meters = 0. END IF IF (PRESENT(y_azimuth_radians)) THEN mp_user_xy = .TRUE. mp_y_azimuth_radians = y_azimuth_radians ELSE mp_y_azimuth_radians = 0. END IF IF (mp_user_xy) THEN ! from "natural" (x,y) system of projection plane to user system: mp_pp_2_xy_matrix(1,1) = COS(mp_y_azimuth_radians) mp_pp_2_xy_matrix(2,2) = mp_pp_2_xy_matrix(1,1) mp_pp_2_xy_matrix(2,1) = SIN(mp_y_azimuth_radians) mp_pp_2_xy_matrix(1,2) = -mp_pp_2_xy_matrix(2,1) mp_pp_2_xy_vector(1) = mp_x_projpoint_meters mp_pp_2_xy_vector(2) = mp_y_projpoint_meters ! from user (x,y) system to "natural" system of projection plane: mp_xy_2_pp_matrix(1,1) = COS(mp_y_azimuth_radians) mp_xy_2_pp_matrix(2,2) = mp_xy_2_pp_matrix(1,1) mp_xy_2_pp_matrix(1,2) = SIN(mp_y_azimuth_radians) mp_xy_2_pp_matrix(2,1) = -mp_xy_2_pp_matrix(1,2) mp_xy_2_pp_vector(1) = -mp_xy_2_pp_matrix(1,1)*mp_x_projpoint_meters & & -mp_xy_2_pp_matrix(1,2)*mp_y_projpoint_meters mp_xy_2_pp_vector(2) = -mp_xy_2_pp_matrix(2,1)*mp_x_projpoint_meters & & -mp_xy_2_pp_matrix(2,2)*mp_y_projpoint_meters END IF END SUBROUTINE Set_Gnomonic SUBROUTINE New_L45_Path (level, uvec) ! Level 5 accepts unit-vectors (uvec's) in a unit sphere from ! the user. ! Level 4 should only be called internally by other routines ! (especially Process_L5_Paths, as it builds L4 paths after ! windowing and duplicating the paths as needed). IMPLICIT NONE INTEGER, INTENT(IN) :: level REAL, DIMENSION(3), INTENT(IN) :: uvec INTEGER :: path IF ((level == 4).OR.(level == 5)) THEN IF (ai_page_open) THEN IF (mp_zoom_defined) THEN IF (mp_projection_defined) THEN IF (.NOT. ai_in_path) THEN ! USUAL CASE: ai_in_path = .TRUE. path = Next_Free_Path() ! Note: this function will check for full library. ai_current_path_index = path ai_current_path_level = level ai_path_level(path) = level ai_Ln_paths(level) = ai_Ln_paths(level) + 1 ai_total_paths = ai_total_paths + 1 ai_segments(path) = 0 ai_pathlib(1:6,0,path) = 0.0 ai_pathlib(1:3,0,path) = uvec mp_last_uvec = uvec ! (memory) ELSE ! ai_in_path already WRITE (*,"(' ERROR: Cannot New_L45_Path with a path open.')") CALL Traceback END IF ELSE ! projection not defined yet WRITE (*,"(' ERROR: Must set map projection before New_L45_Path.')") CALL Traceback END IF ELSE ! zoom from L3 <--> L2 not defined WRITE (*,"(' ERROR: Must Set_Zoom before New_L45_Path.')") CALL Traceback END IF ELSE ! no page open WRITE (*,"(' ERROR: Cannot New_L45_Path before Begin_Page.')") CALL Traceback END IF ELSE WRITE (*,"(' ERROR: Level ',I2,' not allowed in New_L45_Path.')") level CALL Traceback END IF END SUBROUTINE New_L45_Path SUBROUTINE Great_To_L45 (to_uvec) ! Adds a great-circle arc to the current path (on level 4 or 5), ! from the current point (mp_last_uvec) to the input "to_uvec". IMPLICIT NONE REAL, DIMENSION(3), INTENT(IN) :: to_uvec INTEGER :: next, path IF (ai_in_path) THEN IF ((ai_current_path_level == 4).OR.(ai_current_path_level == 5)) THEN path = ai_current_path_index IF (ai_segments(path) < ai_longest) THEN ! just ignore degenerate arcs (of zero length): IF (to_uvec(1) == mp_last_uvec(1)) THEN IF (to_uvec(2) == mp_last_uvec(2)) THEN IF (to_uvec(3) == mp_last_uvec(3)) RETURN END IF END IF IF ((to_uvec(1) == -mp_last_uvec(1)).AND. & &(to_uvec(2) == -mp_last_uvec(2)).AND. & &(to_uvec(3) == -mp_last_uvec(3))) THEN ! reject great-circle to antipode (nonunique) WRITE (*,"(' ERROR: Cannot Great_To_L45 to antipode.')") CALL Traceback ELSE ! USUAL CASE: next = ai_segments(path) + 1 ai_segments(path) = next ai_pathlib(4:6,next,path) = 0.0 ai_pathlib(1:3,next,path) = to_uvec ai_bent(next,path) = .FALSE. mp_last_uvec = to_uvec ! (memory) END IF ELSE WRITE (*,"(' ERROR: Increase ai_longest from current ',I6)") ai_longest CALL Traceback END IF ELSE WRITE (*,"(' ERROR: Cannot Great_To_L45 in path on level ',I2)") & & ai_current_path_level CALL Traceback END IF ELSE WRITE (*,"(' ERROR: Cannot Great_To_L45 before New_L45_Path.')") CALL Traceback END IF END SUBROUTINE Great_To_L45 SUBROUTINE Small_Through_L45 (through_uvec, to_uvec) ! Adds a small-circle arc to the current path (on level 4 or 5), ! from the current point (mp_last_uvec) to the point "to_uvec", ! passing through the intermediate point "through_uvec". ! Note that the final point to_uvec may NOT be the present point, ! because the result is undefined. ! If intermediate point through_uvec is equal to the initial ! or the final point, then the shorter great circle arc is used. ! If the start point is too close to the end point, then also ! a shorter great circle arc is substituted; this is not ! required by logic, but instead is to guard against artifacts ! from roundoff effects when very close uvecs are differenced. IMPLICIT NONE REAL, DIMENSION(3), INTENT(IN) :: through_uvec, to_uvec INTEGER :: next, path REAL, PARAMETER :: tolerance = 0.001 ! limit for small circle arc, in radians (= 6.4 km). REAL :: arc1 REAL, DIMENSION(3) :: from_uvec, & ! copy of mp_last_uvec, to see while debugging & pole_uvec, & ! local copy, computed here & vec, uvec1, uvec2 IF (ai_in_path) THEN IF ((ai_current_path_level == 4).OR.(ai_current_path_level == 5)) THEN from_uvec(1:3) = mp_last_uvec(1:3) ! otherwise not visible in debugger !test for illegal to_uvec == from_uvec: IF (from_uvec(1) == to_uvec(1)) THEN IF (from_uvec(2) == to_uvec(2)) THEN IF (from_uvec(3) == to_uvec(3)) THEN WRITE (*,"(' ERROR: Cannot Small_Through_L45 back to present location; nonunique.')") CALL Traceback END IF END IF END IF !test for to_uvec very close to from_uvec: vec = from_uvec - to_uvec arc1 = Magnitude (vec) IF (arc1 < tolerance) THEN CALL Great_to_L45 (to_uvec) RETURN END IF !test for simple case through_uvec == from_uvec: IF (through_uvec(1) == from_uvec(1)) THEN IF (through_uvec(2) == from_uvec(2)) THEN IF (through_uvec(3) == from_uvec(3)) THEN CALL Great_to_L45 (to_uvec) RETURN END IF END IF END IF !test for simple case through_uvec == to_uvec: IF (through_uvec(1) == to_uvec(1)) THEN IF (through_uvec(2) == to_uvec(2)) THEN IF (through_uvec(3) == to_uvec(3)) THEN CALL Great_to_L45 (to_uvec) RETURN END IF END IF END IF !test for space in path library: path = ai_current_path_index IF (ai_segments(path) < ai_longest) THEN !NORMAL CASE BEGINS HERE: !first, compute pole_uvec: vec(1:3) = through_uvec(1:3) - from_uvec(1:3) CALL Make_Uvec (vec, uvec1) vec(1:3) = to_uvec(1:3) - through_uvec(1:3) CALL Make_Uvec (vec, uvec2) CALL Cross (uvec1, uvec2, vec) CALL Make_Uvec (vec, pole_uvec) !then, write to path library: next = ai_segments(path) + 1 ai_segments(path) = next ai_pathlib(1:3,next,path) = pole_uvec ai_pathlib(4:6,next,path) = to_uvec ai_bent(next,path) = .TRUE. mp_last_uvec = to_uvec ! (memory) ELSE WRITE (*,"(' ERROR: Increase ai_longest from current ',I6)") ai_longest CALL Traceback END IF ELSE WRITE (*,"(' ERROR: Cannot Small_Through_L45 in path on level ',I2)") & & ai_current_path_level CALL Traceback END IF ELSE WRITE (*,"(' ERROR: Cannot Small_Through_L45 before New_L45_Path.')") CALL Traceback END IF END SUBROUTINE Small_Through_L45 SUBROUTINE Small_To_L45 (pole_uvec, to_uvec) ! Adds a small-circle arc to the current path (on level 4 or 5), ! from the current point (mp_last_uvec) to the input "to_uvec", ! swinging counterclockwise (as seen from outside the sphere) ! around "pole_uvec". (To swing clockwise, reverse each ccordinate ! of the pole to get a new pole at the antipode.) ! NOTE: Unlike Great_To_L45, this routine plots a complete ! circle when the goal point equals the current point! ! You can use it to plot a complete great circle by choosing ! an appropriate pole 90 degrees away. IMPLICIT NONE REAL, DIMENSION(3), INTENT(IN) :: pole_uvec, to_uvec INTEGER :: next, path REAL, PARAMETER :: tolerance = 0.04 ! precision for circle radius, radians (= 2.3 deg = 255 km). !Note that tolerances of less than 2.3 degrees cause FPS data from Harvard CMT to be !unacceptably inconsistent, due to rounding of trend and plunge to nearest degree. REAL :: arc1, arc2 REAL, DIMENSION(3) :: from_uvec ! copy of mp_last_uvec, to see while debugging IF (ai_in_path) THEN IF ((ai_current_path_level == 4).OR.(ai_current_path_level == 5)) THEN path = ai_current_path_index IF (ai_segments(path) < ai_longest) THEN from_uvec(1:3) = mp_last_uvec(1:3) ! otherwise not visible arc1 = Arc ( pole_uvec, from_uvec ) arc2 = Arc ( pole_uvec, to_uvec ) IF (ABS(arc1 - arc2) <= tolerance) THEN ! USUAL CASE: next = ai_segments(path) + 1 ai_segments(path) = next ai_pathlib(1:3,next,path) = pole_uvec ai_pathlib(4:6,next,path) = to_uvec ai_bent(next,path) = .TRUE. mp_last_uvec = to_uvec ! (memory) ELSE ! inconsistent data; not a small circle! WRITE (*,"(' ERROR: Initial and final radii of small circle in Small_To_L45' & &/' are ',F6.4,' and ',F6.4,', which disagree by more than' & &/' the permitted tolerance of ',F6.4,' radians.')") arc1, arc2, tolerance CALL Traceback END IF ELSE WRITE (*,"(' ERROR: Increase ai_longest from current ',I6)") ai_longest CALL Traceback END IF ELSE WRITE (*,"(' ERROR: Cannot Small_To_L45 in path on level ',I2)") & & ai_current_path_level CALL Traceback END IF ELSE WRITE (*,"(' ERROR: Cannot Small_To_L45 before New_L45_Path.')") CALL Traceback END IF END SUBROUTINE Small_To_L45 SUBROUTINE End_L45_Path (close, stroke, fill) IMPLICIT NONE LOGICAL, INTENT(IN) :: close, stroke, fill INTEGER :: level, path REAL, DIMENSION(3) :: t_uvec ! Note: The .AI macros which include "closepath" ! don't close it as one would like, with an extra ! "lineto". Instead, they do a curveto, which somehow ! misses the last point you specified. ! Therefore, it is not enough to specify close = T, ! you must also program so that the path actually ! comes back to its initial point! IF (ai_in_path) THEN path = ai_current_path_index level = ai_current_path_level IF ((level == 4).OR.(level == 5)) THEN IF (stroke .OR. fill) THEN ! USUAL ROUTE: IF (fill) THEN ! Logic of Process_L2_Paths requires these paths ! to end at their initial points: IF ((mp_last_uvec(1) /= ai_pathlib(1,0,path)).OR. & &(mp_last_uvec(2) /= ai_pathlib(2,0,path)).OR. & &(mp_last_uvec(3) /= ai_pathlib(3,0,path))) THEN t_uvec = ai_pathlib(1:3,0,path) CALL Great_To_L45(t_uvec) END IF ! completion of path needed END IF ! filled path ai_closed(path) = close ai_stroked(path) = stroke ai_filled(path) = fill ai_in_path = .FALSE. mp_last_uvec = (/0.,0.,0./) ! (undefined) IF (level == 4) THEN CALL Process_L4_Paths ELSE ! level 5 CALL Process_L5_Paths END IF ELSE !neither stroked nor filled WRITE (*,"(' ERROR: End_L45_path with neither fill nor stroke.')") CALL Traceback END IF ELSE WRITE (*,"(' ERROR: Cannot End_L45_path when other level open.')") CALL Traceback END IF ELSE WRITE (*,"(' ERROR: Cannot End_L45_Path when path not open.')") CALL Traceback END IF END SUBROUTINE End_L45_Path SUBROUTINE New_L67_Path (level, r1, r2) ! Accepts initial points at level 6 or 7 (on the spherical surface, ! with position expressed as either: ! Level 6: r1 = theta (colatitude from North, radians); ! r2 = phi (East longitude in radians); OR, ! Level 7: r1 = lon (East longitude in degrees); ! r2 = lat (North latitude in degrees). IMPLICIT NONE INTEGER, INTENT(IN) :: level REAL, INTENT(IN) :: r1, r2 INTEGER :: path REAL :: lat, lon, phi, theta IF ((level == 6).OR.(level == 7)) THEN IF (ai_page_open) THEN IF (mp_zoom_defined) THEN IF (mp_projection_defined) THEN IF (.NOT. ai_in_path) THEN ! USUAL CASE: ai_in_path = .TRUE. path = Next_Free_Path() ! Note: this function will check for full library. ai_current_path_index = path ai_current_path_level = level ai_path_level(path) = level ai_Ln_paths(level) = ai_Ln_paths(level) + 1 ai_total_paths = ai_total_paths + 1 ai_segments(path) = 0 ai_pathlib(3:6,0,path) = 0.0 IF (level == 6) THEN theta = r1 phi = r2 IF ((theta < 0.).OR.(theta > Pi)) THEN WRITE (*,"(' ERROR: Illegal colatitude theta for New_L67_Path.')") CALL Traceback ELSE ai_pathlib(1,0,path) = theta ai_pathlib(2,0,path) = phi mp_last_theta = theta mp_last_phi = phi CALL ThetaPhi_2_Uvec (theta, phi, mp_last_uvec) END IF ELSE ! level == 7 lon = r1 lat = r2 IF ((lat > 90.).OR.(lat < -90.)) THEN WRITE (*,"(' ERROR: Illegal latitude theta for New_L67_Path.')") CALL Traceback ELSE ai_pathlib(1,0,path) = lon ai_pathlib(2,0,path) = lat mp_last_lon = lon mp_last_lat = lat CALL LonLat_2_Uvec (lon, lat, mp_last_uvec) END IF END IF ELSE ! ai_in_path already WRITE (*,"(' ERROR: Cannot New_L67_Path with a path open.')") CALL Traceback END IF ELSE ! projection not defined yet WRITE (*,"(' ERROR: Must set map projection before New_L67_Path.')") CALL Traceback END IF ELSE ! zoom from L3 <--> L2 not defined WRITE (*,"(' ERROR: Must Set_Zoom before New_L67_Path.')") CALL Traceback END IF ELSE ! no page open WRITE (*,"(' ERROR: Cannot New_L67_Path before Begin_Page.')") CALL Traceback END IF ELSE WRITE (*,"(' ERROR: Level ',I2,' not allowed in New_L67_Path.')") level CALL Traceback END IF END SUBROUTINE New_L67_Path SUBROUTINE Great_To_L67 (r1, r2) ! Adds a great-circle arc to the current path (on level 6 or 7), ! from the current point to (r1, r2), which is either: ! (theta, phi) in radians on level 6, or: ! (lon, lat) in degrees on level 7. ! Theta is measured S from N pole. Phi is East from 0E. ! Lon is positive to E of 0E. Lat is + to N of 0N. IMPLICIT NONE REAL, INTENT(IN) :: r1, r2 INTEGER :: level, next, path REAL :: lat, lon, phi, theta REAL, DIMENSION(3) :: to_uvec IF (ai_in_path) THEN IF ((ai_current_path_level == 6).OR.(ai_current_path_level == 7)) THEN path = ai_current_path_index level = ai_current_path_level IF (ai_segments(path) < ai_longest) THEN IF (level == 6) THEN theta = r1 IF ((theta < 0.).OR.(theta > Pi)) THEN WRITE (*,"(' ERROR: Illegal colatitude to Great_To_L67.')") CALL Traceback END IF phi = r2 ! just ignore degenerate arcs (of zero length): IF (theta == mp_last_theta) THEN IF (phi == mp_last_phi) RETURN END IF ! reject great-circle to antipode (nonunique) CALL ThetaPhi_2_Uvec (theta, phi, to_uvec) IF ((to_uvec(1) == -mp_last_uvec(1)).AND. & &(to_uvec(2) == -mp_last_uvec(2)).AND. & &(to_uvec(3) == -mp_last_uvec(3))) THEN WRITE (*,"(' ERROR: Cannot Great_To_L67 to antipode.')") CALL Traceback ELSE ! USUAL CASE: next = ai_segments(path) + 1 ai_segments(path) = next ai_pathlib(3:6,next,path) = 0.0 ai_pathlib(1,next,path) = theta ai_pathlib(2,next,path) = phi ai_bent(next,path) = .FALSE. mp_last_uvec = to_uvec ! (memory) mp_last_theta = theta mp_last_phi = phi END IF ELSE ! level == 7 lon = r1 lat = r2 IF ((lat > 90.).OR.(lat < -90.)) THEN WRITE (*,"(' ERROR: Illegal latitude to Great_To_L67.')") CALL Traceback END IF ! just ignore degenerate arcs (of zero length): IF (lon == mp_last_lon) THEN IF (lat == mp_last_lat) RETURN END IF ! reject great-circle to antipode (nonunique) CALL LonLat_2_Uvec (lon, lat, to_uvec) IF ((to_uvec(1) == -mp_last_uvec(1)).AND. & &(to_uvec(2) == -mp_last_uvec(2)).AND. & &(to_uvec(3) == -mp_last_uvec(3))) THEN WRITE (*,"(' ERROR: Cannot Great_To_L67 to antipode.')") CALL Traceback ELSE ! USUAL CASE: next = ai_segments(path) + 1 ai_segments(path) = next ai_pathlib(3:6,next,path) = 0.0 ai_pathlib(1,next,path) = lon ai_pathlib(2,next,path) = lat ai_bent(next,path) = .FALSE. mp_last_uvec = to_uvec ! (memory) mp_last_lon = lon mp_last_lat = lat END IF END IF ELSE WRITE (*,"(' ERROR: Increase ai_longest from current ',I6)") ai_longest CALL Traceback END IF ELSE WRITE (*,"(' ERROR: Cannot Great_To_L67 in path on level ',I2)") & & ai_current_path_level CALL Traceback END IF ELSE WRITE (*,"(' ERROR: Cannot Great_To_L67 before New_L67_Path.')") CALL Traceback END IF END SUBROUTINE Great_To_L67 SUBROUTINE Small_To_L67 (r1, r2, r3, r4) ! Adds a small-circle arc to the current path (on level 6 or 7), ! from the current point to (r3, r4), which is either: ! (theta, phi) in radians on level 6, or: ! (lon, lat) in degrees on level 7. ! Theta is measured S from N pole. Phi is East from 0E. ! Lon is positive to E of 0E. Lat is + to N of 0N. ! The pole for the small circle is (r1, r2), which is either ! (theta, phi) or (lon, lat) depending on the level. ! Note that rotation is always counterclockwise about the pole. ! IF (r3, r4) are the same as the current point, then a ! complete circle is drawn. IMPLICIT NONE REAL, INTENT(IN) :: r1, r2, r3, r4 INTEGER :: level, next, path REAL :: lat, lon, phi, pole_lat, pole_lon, pole_phi, pole_theta, theta REAL, DIMENSION(3) :: to_uvec IF (ai_in_path) THEN IF ((ai_current_path_level == 6).OR.(ai_current_path_level == 7)) THEN path = ai_current_path_index level = ai_current_path_level IF (ai_segments(path) < ai_longest) THEN ! USUAL CASE: IF (level == 6) THEN pole_theta = r1 IF ((pole_theta < 0.).OR.(pole_theta > Pi)) THEN WRITE (*,"(' ERROR: Illegal pole colatitude to Small_To_L67.')") CALL Traceback END IF pole_phi = r2 theta = r3 IF ((theta < 0.).OR.(theta > Pi)) THEN WRITE (*,"(' ERROR: Illegal colatitude to Small_To_L67.')") CALL Traceback END IF phi = r4 next = ai_segments(path) + 1 ai_segments(path) = next ai_pathlib(5:6,next,path) = 0.0 ai_pathlib(1,next,path) = pole_theta ai_pathlib(2,next,path) = pole_phi ai_pathlib(3,next,path) = theta ai_pathlib(4,next,path) = phi ai_bent(next,path) = .TRUE. CALL ThetaPhi_2_Uvec (theta, phi, to_uvec) mp_last_uvec = to_uvec ! (memory) mp_last_theta = theta mp_last_phi = phi ELSE ! level == 7 pole_lon = r1 pole_lat = r2 IF ((pole_lat > 90.).OR.(pole_lat < -90.)) THEN WRITE (*,"(' ERROR: Illegal pole latitude to Small_To_L67.')") CALL Traceback END IF lon = r3 lat = r4 IF ((lat > 90.).OR.(lat < -90.)) THEN WRITE (*,"(' ERROR: Illegal latitude to Small_To_L67.')") CALL Traceback END IF next = ai_segments(path) + 1 ai_segments(path) = next ai_pathlib(5:6,next,path) = 0.0 ai_pathlib(1,next,path) = pole_lon ai_pathlib(2,next,path) = pole_lat ai_pathlib(3,next,path) = lon ai_pathlib(4,next,path) = lat ai_bent(next,path) = .TRUE. CALL LonLat_2_Uvec (lon, lat, to_uvec) mp_last_uvec = to_uvec ! (memory) mp_last_lon = lon mp_last_lat = lat END IF ELSE WRITE (*,"(' ERROR: Increase ai_longest from current ',I6)") ai_longest CALL Traceback END IF ELSE WRITE (*,"(' ERROR: Cannot Small_To_L67 in path on level ',I2)") & & ai_current_path_level CALL Traceback END IF ELSE WRITE (*,"(' ERROR: Cannot Small_To_L67 before New_L67_Path.')") CALL Traceback END IF END SUBROUTINE Small_To_L67 SUBROUTINE End_L67_Path (close, stroke, fill) IMPLICIT NONE LOGICAL, INTENT(IN) :: close, stroke, fill INTEGER :: level, path ! Note: The .AI macros which include "closepath" ! don't close it as one would like, with an extra ! "lineto". Instead, they do a curveto, which somehow ! misses the last point you specified. ! Therefore, it is not enough to specify close = T, ! you must also program so that the path actually ! comes back to its initial point! IF (ai_in_path) THEN path = ai_current_path_index level = ai_current_path_level IF ((level == 6).OR.(level == 7)) THEN IF (stroke .OR. fill) THEN ! USUAL ROUTE: IF (fill) THEN ! Logic of Process_L2_Paths requires these paths ! to end at their initial points: IF (level == 6) THEN IF ((mp_last_theta /= ai_pathlib(1,0,path)).OR. & &(mp_last_phi /= ai_pathlib(2,0,path))) THEN CALL Great_To_L67(ai_pathlib(1,0,path),ai_pathlib(2,0,path)) END IF ! completion of path needed ELSE ! level == 7 THEN IF ((mp_last_lon /= ai_pathlib(1,0,path)).OR. & &(mp_last_lat /= ai_pathlib(2,0,path))) THEN CALL Great_To_L67(ai_pathlib(1,0,path),ai_pathlib(2,0,path)) END IF ! completion of path needed END IF ! L6 or L7 END IF ! filled path ai_closed(path) = close ai_stroked(path) = stroke ai_filled(path) = fill ai_in_path = .FALSE. IF (level == 6) THEN CALL Process_L6_Paths ELSE ! level 7 CALL Process_L7_Paths END IF ELSE !neither stroked nor filled WRITE (*,"(' ERROR: End_L67_path with neither fill nor stroke.')") CALL Traceback END IF ELSE WRITE (*,"(' ERROR: Cannot End_L67_path when other level open.')") CALL Traceback END IF ELSE WRITE (*,"(' ERROR: Cannot End_L67_Path when path not open.')") CALL Traceback END IF END SUBROUTINE End_L67_Path ! =============================================== ! | INTERNAL UTILITY ROUTINES | ! | (in alphabetical order ) | ! =============================================== 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 LOGICAL FUNCTION Clockways (level_45_path) ! Determines whether a given level-5 or level-4 path on a sphere ! (expressed as a series of great- and small-circle arcs, with ! control points expressed as unit vectors) travels clockwise or ! not (i.e., counterclockwise). ! It does this by adding up internal and segment-boundary bends ! (to the right is +, to left is -) around the path. ! CAUTION: The property of clockwise/counterclockwise is only ! DEFINED for paths on a sphere which: ! -close by returning to their initial point; ! -divide the sphere into two unequal areas; ! -cross themselves either nowhere (0 points) or at an ! even number of points. ! If a non-conforming path (e.g., open path, great circle, or ! figure-8 path) is processed, results are meaningless! ! The single argument "L45_path" is the integer identifier of ! one path among those in ai_pathlib (module's global data). IMPLICIT NONE INTEGER, INTENT(IN) :: level_45_path INTEGER :: i, j, path, segments REAL :: azim_before, azim_end, azim_start, internal_bend, & & last_lat, last_lon, pazim0, pazim1, size, start_bend, total_bend REAL, DIMENSION(3) :: last_uvec, pole_uvec, t_vec, t0_uvec, t1_uvec path = level_45_path segments = ai_segments(path) IF ((ai_path_level(path) == 5).OR. & &(ai_path_level(path) == 4)) THEN IF (segments > 0) THEN total_bend = 0. last_uvec = ai_pathlib(1:3,0,path) DO j = 1, segments+1 ! one extra time to pick up start_bend i = 1 + MOD(j-1, segments) ! = j, except last time: 1 ! Compute bending (to left) within the arc itself: IF (ai_bent(i,path)) THEN ! small circle arc pole_uvec = ai_pathlib(1:3,i,path) ! pole of small t0_uvec = last_uvec ! initial point of small circle arc t1_uvec = ai_pathlib(4:6,i,path) ! end point size = Arc(pole_uvec, t1_uvec) ! radius in arc radians IF ((t0_uvec(1) == t1_uvec(1)).AND. & &(t0_uvec(2) == t1_uvec(2)).AND. & &(t0_uvec(3) == t1_uvec(3))) THEN ! complete circle pazim0 = Two_Pi pazim1 = 0. ! one counterclockwise (negative; left) circuit azim_start = Pi_over_2 + Relative_Compass (t0_uvec, pole_uvec) azim_end = azim_start ! same trend as start ELSE ! incomplete circle pazim0 = Relative_Compass (pole_uvec, t0_uvec) pazim1 = Relative_Compass (pole_uvec, t1_uvec) IF (pazim1 >= pazim0) pazim1 = pazim1 - Two_Pi ! because, by definition, small circle arcs ! are ALWAYS counterclockwise about their pole, so ! azimuth (measured at pole) should decrease. azim_start = Pi_over_2 + Relative_Compass (t0_uvec, pole_uvec) azim_end = Pi_over_2 + Relative_Compass (t1_uvec, pole_uvec) END IF internal_bend = (pazim1 - pazim0) * COS(size) ELSE ! great circle arc; COS(size) = 0. internal_bend = 0. t0_uvec = last_uvec t1_uvec = ai_pathlib(1:3,i,path) CALL Cross (t0_uvec, t1_uvec, t_vec) CALL Make_Uvec (t_vec, pole_uvec) azim_start = Pi_over_2 + Relative_Compass (t0_uvec, pole_uvec) azim_end = Pi_over_2 + Relative_Compass (t1_uvec, pole_uvec) END IF IF (j <= segments) THEN total_bend = total_bend + internal_bend END IF IF (j > 1) THEN ! azim_before is known start_bend = azim_start - azim_before ! keep in range of +Pi to -Pi (near 0.) IF (start_bend > Pi) start_bend = start_bend - Two_Pi IF (start_bend < -Pi) start_bend = start_bend + Two_Pi total_bend = total_bend + start_bend END IF ! set memory before processing next segment: azim_before = azim_end IF (ai_bent(i,path)) THEN last_uvec = ai_pathlib(4:6,i,path) ELSE last_uvec = ai_pathlib(1:3,i,path) END IF END DO ! i = 1, segments Clockways = (total_bend > 0.) ELSE WRITE (*,"(' ERROR: Single-point paths cannot have the Clockways property.')") CALL Traceback END IF ELSE ! wrong path level WRITE (*,"(' ERROR: Wrong path level for argument of Clockways: ',I2)") ai_path_level(path) CALL Traceback END IF END FUNCTION Clockways 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 REAL FUNCTION Conformal_Deflation (uvec) ! Gives the local "deflation factor" which CAN be used to scale ! down the size of an object at "uvec" (on the sphere) BEFORE it is ! plotted, so that AFTER it is plotted its size is not affected by ! the map projection. ! For example, using a standard Mercator projection, Deflation ! is 1.00 along the equator, but shrinks to = COS(latitude) ! at other latitudes. ! Note that this correction only works for the conformal ! map projections: Mercator, Lambert Conformal Conic, and ! Stereographic! In other projections, we return ! 1.00 (no correction) because this is an unsolvable problem. IMPLICIT NONE REAL, DIMENSION(3), INTENT(IN) :: uvec REAL :: colatitude, latitude SELECT CASE(mp_projection_number) CASE (1) ! Mercator colatitude = Arc (uvec, mp_pole_3_uvec) Conformal_Deflation = SIN(colatitude) CASE (2) ! Lambert Conformal Conic ! mp_proj2_F = COS(latitude_1_radians) * TAN(Pi_over_4 + latitude_1_radians/2.)**mp_proj2_n / mp_proj2_n latitude = Pi_over_2 - Arc(uvec, mp_pole_3_uvec) Conformal_Deflation = COS(latitude)*TAN(0.5*(Pi_over_2 + latitude))**mp_proj2_n / & & (mp_proj2_F * mp_proj2_n) CASE (6) ! Stereographic colatitude = Arc (uvec, mp_projpoint_uvec) Conformal_Deflation = (1. + COS(colatitude)) * 0.5 CASE DEFAULT Conformal_Deflation = 1.00 END SELECT END FUNCTION Conformal_Deflation SUBROUTINE Corners (projpoint_uvec) ! After Set_[projection] has defined the (up to) 4 walls of ! the projectable "rind" on the sphere, and before Trace_Rind() ! is asked to draw around it, this routine defines the topology. ! It fills in two set of global variables: ! INTEGER, DIMENSION(0:1, 4) :: mp_next_wall ! has 4 columns for the four walls of the rind, in the ! order: 1 = bottom, 2 = right, 3 = top, 4 = left. ! It has two rows (0:1,) for the low-t_wall and the ! high t-wall ends of each wall. Low-t_wall occurs at ! the left end of the bottom, the bottom end of the right, ! the right end of the top, and the top end of the right ! wall. (Note that these are in reverse order from ! the vectors "first_x_uvec" and "last_x_uvec" below!) ! A value of 0 indicates that the wall is a complete ! circle and connects only to itself at that end. ! A value of 1,2,3,4 indicates that it connects to the ! next wall of that number (by the system above). ! REAL, DIMENSION(3) :: mp_first_top_uvec, mp_last_top_uvec, & ! & mp_first_left_uvec, mp_last_left_uvec, & ! & mp_first_right_uvec, mp_last_right_uvec, & ! & mp_first_bottom_uvec, mp_last_bottom_uvec ! These 8 unit vectors give "first" and "last" points on ! each of the 4 walls of the rind, meaning corner points ! where one wall ceases to be relevant and another takes over. ! (The parts of the wall outside these limits are not relevant ! either to finding intersections with paths, or to tracing ! around the rind.) First and last refer to a counterclockwise ! circling about the corresponding mp_x_uvec pole. Since ! these poles all point OUTWARD from the rind, the "first" ! end corresponds to the "high-t_wall end" above, and the ! "last" end corresponds to the "low-t_wall end". IMPLICIT NONE REAL, DIMENSION(3), INTENT(IN) :: projpoint_uvec INTEGER :: i, number LOGICAL :: overlap, using_wall REAL, PARAMETER :: tolerance = 0.00001 REAL :: close1, close2, dot1a, dot1b, dot2a, dot2b, error, error1, & & error2, maxdot, sidestep REAL, DIMENSION(3), PARAMETER :: north_uvec = (/ 1., 0., 0. /) REAL, DIMENSION(3), PARAMETER :: zero_vec = (/ 0., 0., 0. /) ! (This is a useful flag in "CALL Circles_Intersect" that we ! want to consider the whole circle, not an arc.) REAL, DIMENSION(3) :: first_uvec, perp1_uvec, & & perp2_uvec, point1_uvec, point2_uvec, pole_uvec, & & pole_applies_uvec, right_uvec, t_vec, trial_uvec ! Decide which way is "right" from projpoint; "East" of NorthEast_Convention: CALL NorthEast_Convention (projpoint_uvec, t_vec, right_uvec) ! Initialize intersections to "none": mp_next_wall = 0 ! whole array; default is no connections ! Initialize "first" and "last" points on each circle as equal. !(Also, check wall definitions for some easy-to-spot errors!) DO i = 1, 4 SELECT CASE(i) CASE(1) using_wall = mp_using_bottom pole_uvec = mp_bottom_uvec pole_applies_uvec = mp_bottom_applies_uvec maxdot = mp_bottom_maxdot CASE(2) using_wall = mp_using_right pole_uvec = mp_right_uvec pole_applies_uvec = mp_right_applies_uvec maxdot = mp_right_maxdot CASE(3) using_wall = mp_using_top pole_uvec = mp_top_uvec pole_applies_uvec = mp_top_applies_uvec maxdot = mp_top_maxdot CASE(4) using_wall = mp_using_left pole_uvec = mp_left_uvec pole_applies_uvec = mp_left_applies_uvec maxdot = mp_left_maxdot END SELECT IF (using_wall) THEN ! verify maxdot is in range: 0. <= maxdot < 1. IF ((maxdot < 0.).OR.(maxdot >= 1.)) THEN WRITE (*,"(' ERROR: Violation of limits 0. <= maxdot < 1. on rind side ',I1)") i CALL Traceback END IF ! check that this wall does not exclude projection point! IF (Dot(pole_applies_uvec,projpoint_uvec) > 0.) THEN IF (Dot(projpoint_uvec,pole_uvec) > maxdot) THEN WRITE (*,"(' ERROR: Projection point excluded by rind wall ',I1)") i CALL Traceback END IF END IF ! find "any old" uvec in the wall for a default first/last: sidestep = SQRT(1. - maxdot**2) CALL Cross (projpoint_uvec, pole_uvec, t_vec) IF (Length(t_vec) == 0.) CALL Cross (north_uvec, pole_uvec, t_vec) IF (Length(t_vec) == 0.) CALL Cross (right_uvec, pole_uvec, t_vec) CALL Make_Uvec(t_vec, perp1_uvec) CALL Cross (perp1_uvec, pole_uvec, perp2_uvec) ! Now, perp2_uvec is perp. to pole, about opposite projpoint. first_uvec = maxdot*pole_uvec + sidestep*perp2_uvec ELSE first_uvec = zero_vec END IF SELECT CASE(i) CASE(1) mp_first_bottom_uvec = first_uvec mp_last_bottom_uvec = mp_first_bottom_uvec CASE(2) mp_first_right_uvec = first_uvec mp_last_right_uvec = mp_first_right_uvec CASE(3) mp_first_top_uvec = first_uvec mp_last_top_uvec = mp_first_top_uvec CASE(4) mp_first_left_uvec = first_uvec mp_last_left_uvec = mp_first_left_uvec END SELECT END DO ! Survey all four sides (if used) for intersections with the others: IF (mp_using_top) THEN IF (mp_using_left) THEN ! look for left\top corner CALL Circles_Intersect (mp_top_uvec, mp_top_maxdot, zero_vec, zero_vec, & & mp_left_uvec, mp_left_maxdot, zero_vec, zero_vec, & ! input & overlap, number, point1_uvec, point2_uvec) ! output IF (overlap) THEN WRITE (*,"(' ERROR: Walls of rind may not overlap.')") CALL Traceback END IF IF (number > 0) THEN ! found some possible points ! choose which answer is most reasonable (least error) dot1a = Dot(point1_uvec, mp_top_applies_uvec) dot1b = Dot(point1_uvec, mp_left_applies_uvec) error1 = -MIN(dot1a, dot1b) ! problem if either is <0. error = error1 trial_uvec = point1_uvec IF (number == 2) THEN dot2a = Dot(point2_uvec, mp_top_applies_uvec) dot2b = Dot(point2_uvec, mp_left_applies_uvec) error2 = -MIN(dot2a, dot2b) ! problem if either is <0. IF (ABS(error2 - error1) <= tolerance) THEN ! essentially equal ! choose the one closer to projpoint: close1 = Dot(point1_uvec, projpoint_uvec) close2 = Dot(point2_uvec, projpoint_uvec) IF (close2 > close1) THEN error = error2 trial_uvec = point2_uvec END IF ELSE IF (error2 < error1) THEN ! second point is a better solution error = error2 trial_uvec = point2_uvec END IF IF (error <= tolerance) THEN ! ***** Found the left \ top corner! ******* mp_next_wall(1,3) = 4 mp_first_top_uvec = trial_uvec mp_next_wall(0,4) = 3 mp_last_left_uvec = trial_uvec ! ********************************************** END IF END IF ! 2 possible intersection points END IF ! got candidate points END IF ! mp_using_left IF (mp_using_right) THEN ! look for top/right corner CALL Circles_Intersect (mp_top_uvec, mp_top_maxdot, zero_vec, zero_vec, & & mp_right_uvec, mp_right_maxdot, zero_vec, zero_vec, & ! input & overlap, number, point1_uvec, point2_uvec) ! output IF (overlap) THEN WRITE (*,"(' ERROR: Walls of rind may not overlap.')") CALL Traceback END IF IF (number > 0) THEN ! found some possible points ! choose which answer is most reasonable (least error) dot1a = Dot(point1_uvec, mp_top_applies_uvec) dot1b = Dot(point1_uvec, mp_right_applies_uvec) error1 = -MIN(dot1a, dot1b) ! problem if either is <0. error = error1 trial_uvec = point1_uvec IF (number == 2) THEN dot2a = Dot(point2_uvec, mp_top_applies_uvec) dot2b = Dot(point2_uvec, mp_right_applies_uvec) error2 = -MIN(dot2a, dot2b) ! problem if either is <0. IF (ABS(error2 - error1) <= tolerance) THEN ! essentially equal ! choose the one closer to projpoint: close1 = Dot(point1_uvec, projpoint_uvec) close2 = Dot(point2_uvec, projpoint_uvec) IF (close2 > close1) THEN error = error2 trial_uvec = point2_uvec END IF ELSE IF (error2 < error1) THEN ! second point is a better solution error = error2 trial_uvec = point2_uvec END IF IF (error <= tolerance) THEN ! ***** Found the top / right corner! ******* mp_next_wall(0,3) = 2 mp_last_top_uvec = trial_uvec mp_next_wall(1,2) = 3 mp_first_right_uvec = trial_uvec ! ********************************************** END IF END IF ! 2 possible intersection points END IF ! got candidate points END IF ! mp_using_right IF (mp_using_bottom) THEN ! If no connections yet, use any top/bottom intersections ! at the left and right ends of a "lens" shape: CALL Circles_Intersect (mp_top_uvec, mp_top_maxdot, zero_vec, zero_vec, & & mp_bottom_uvec, mp_bottom_maxdot, zero_vec, zero_vec, & ! input & overlap, number, point1_uvec, point2_uvec) ! output IF (overlap) THEN WRITE (*,"(' ERROR: Walls of rind may not overlap.')") CALL Traceback END IF IF (number == 2) THEN ! found two points ! decide which is right, which is left IF (Dot(point1_uvec,right_uvec) > Dot(point2_uvec,right_uvec)) THEN ! swap t_vec = point1_uvec point1_uvec = point2_uvec point2_uvec = t_vec END IF !Now, point1 is to left, point2 to right. IF (mp_next_wall(1,3) == 0) THEN ! ********************************************** ! no left \ top; need left corner: mp_next_wall(1,3) = 1 mp_first_top_uvec = point1_uvec mp_next_wall(0,1) = 3 mp_last_bottom_uvec = point1_uvec ! ********************************************** END IF IF (mp_next_wall(0,3) == 0) THEN ! ********************************************** ! no top / right; need right corner: mp_next_wall(0,3) = 1 mp_last_top_uvec = point2_uvec mp_next_wall(1,1) = 3 mp_first_bottom_uvec = point2_uvec ! ********************************************** END IF END IF ! found 2 top/bottom intersections; used if needed END IF ! mp_using_bottom END IF ! mp_using_top ! Have now found all possible intersections involving the top. ! From now on, only consider those which don't. IF (mp_using_bottom) THEN IF (mp_using_left) THEN ! look for left / bottom corner CALL Circles_Intersect (mp_bottom_uvec, mp_bottom_maxdot, zero_vec, zero_vec, & & mp_left_uvec, mp_left_maxdot, zero_vec, zero_vec, & ! input & overlap, number, point1_uvec, point2_uvec) ! output IF (overlap) THEN WRITE (*,"(' ERROR: Walls of rind may not overlap.')") CALL Traceback END IF IF (number > 0) THEN ! found some possible points ! choose which answer is most reasonable (least error) dot1a = Dot(point1_uvec, mp_bottom_applies_uvec) dot1b = Dot(point1_uvec, mp_left_applies_uvec) error1 = -MIN(dot1a, dot1b) ! problem if either is <0. error = error1 trial_uvec = point1_uvec IF (number == 2) THEN dot2a = Dot(point2_uvec, mp_bottom_applies_uvec) dot2b = Dot(point2_uvec, mp_left_applies_uvec) error2 = -MIN(dot2a, dot2b) ! problem if either is <0. IF (ABS(error2 - error1) <= tolerance) THEN ! essentially equal ! choose the one closer to projpoint: close1 = Dot(point1_uvec, projpoint_uvec) close2 = Dot(point2_uvec, projpoint_uvec) IF (close2 > close1) THEN error = error2 trial_uvec = point2_uvec END IF ELSE IF (error2 < error1) THEN ! second point is a better solution error = error2 trial_uvec = point2_uvec END IF IF (error <= tolerance) THEN ! ***** Found the left / bottom corner! ******* mp_next_wall(0,1) = 4 mp_last_bottom_uvec = trial_uvec mp_next_wall(1,4) = 1 mp_first_left_uvec = trial_uvec ! ********************************************** END IF END IF ! 2 possible intersection points END IF ! got candidate points END IF ! mp_using_left IF (mp_using_right) THEN ! look for bottom \ right corner CALL Circles_Intersect (mp_bottom_uvec, mp_bottom_maxdot, zero_vec, zero_vec, & & mp_right_uvec, mp_right_maxdot, zero_vec, zero_vec, & ! input & overlap, number, point1_uvec, point2_uvec) ! output IF (overlap) THEN WRITE (*,"(' ERROR: Walls of rind may not overlap.')") CALL Traceback END IF IF (number > 0) THEN ! found some possible points ! choose which answer is most reasonable (least error) dot1a = Dot(point1_uvec, mp_bottom_applies_uvec) dot1b = Dot(point1_uvec, mp_right_applies_uvec) error1 = -MIN(dot1a, dot1b) ! problem if either is <0. error = error1 trial_uvec = point1_uvec IF (number == 2) THEN dot2a = Dot(point2_uvec, mp_bottom_applies_uvec) dot2b = Dot(point2_uvec, mp_right_applies_uvec) error2 = -MIN(dot2a, dot2b) ! problem if either is <0. IF (ABS(error2 - error1) <= tolerance) THEN ! essentially equal ! choose the one closer to projpoint: close1 = Dot(point1_uvec, projpoint_uvec) close2 = Dot(point2_uvec, projpoint_uvec) IF (close2 > close1) THEN error = error2 trial_uvec = point2_uvec END IF ELSE IF (error2 < error1) THEN ! second point is a better solution error = error2 trial_uvec = point2_uvec END IF IF (error <= tolerance) THEN ! ***** Found the bottom \ right corner! ******* mp_next_wall(1,1) = 2 mp_first_bottom_uvec = trial_uvec mp_next_wall(0,2) = 1 mp_last_right_uvec = trial_uvec ! ********************************************** END IF END IF ! 2 possible intersection points END IF ! got candidate points END IF ! mp_using_right END IF ! mp_using_bottom !Finally, check for zenith and nadir corners, if needed IF (mp_using_left) THEN IF (mp_using_right) THEN CALL Circles_Intersect (mp_left_uvec, mp_left_maxdot, zero_vec, zero_vec, & & mp_right_uvec, mp_right_maxdot, zero_vec, zero_vec, & ! input & overlap, number, point1_uvec, point2_uvec) ! output IF (overlap) THEN WRITE (*,"(' ERROR: Walls of rind may not overlap.')") CALL Traceback END IF IF (number == 2) THEN ! found two points ! decide which is zenith (up), which nadir (down) IF (Dot(point1_uvec,north_uvec) > Dot(point2_uvec,north_uvec)) THEN ! swap t_vec = point1_uvec point1_uvec = point2_uvec point2_uvec = t_vec END IF !Now, point1 is to nadir (down), point2 to zenth (up). IF (mp_next_wall(1,4) == 0) THEN ! ********************************************** ! no left \ bottom; need nadir corner mp_next_wall(1,4) = 2 mp_first_left_uvec = point1_uvec mp_next_wall(0,2) = 4 mp_last_right_uvec = point1_uvec ! ********************************************** END IF IF (mp_next_wall(0,4) == 0) THEN ! ********************************************** ! no left \ top; need zenith corner: mp_next_wall(0,4) = 2 mp_last_left_uvec = point2_uvec mp_next_wall(1,2) = 4 mp_first_right_uvec = point2_uvec ! ********************************************** END IF END IF ! found 2 zenith/nadir intersections; used if needed END IF ! mp_using_right END IF ! mp_using_left END SUBROUTINE Corners 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 CHARACTER*1 FUNCTION In_Rind (uvec) ! The "Rind" is the part of the planet surface that is ! projectable (within any limits and outside any cut wedges). ! Global variables at module level contain these limits. ! "Uvec" is the input unit-vector representing postion ! on the surface of a unit sphere representing the planet. ! Returns 'I' for inside, 'O' for outside, 'B' for on-boundary. IMPLICIT NONE REAL, DIMENSION(3), INTENT(IN) :: uvec INTEGER :: i REAL :: upper_limit, value In_Rind = 'I' ! but, may be changed below if any test fails DO i = 1, 4 SELECT CASE(i) CASE(1) ! top boundary IF (mp_using_top.AND.(Dot(uvec, mp_top_applies_uvec) > 0.)) THEN value = Dot (uvec, mp_top_uvec) upper_limit = mp_top_maxdot ELSE ! top limit doesn't apply; give a "free pass" value = 0. upper_limit = 1. END IF CASE(2) ! left boundary IF (mp_using_left.AND.(Dot(uvec, mp_left_applies_uvec) > 0.)) THEN value = Dot (uvec, mp_left_uvec) upper_limit = mp_left_maxdot ELSE ! left limit doesn't apply; give a "free pass" value = 0. upper_limit = 1. END IF CASE(3) ! right boundary IF (mp_using_right.AND.(Dot(uvec, mp_right_applies_uvec) > 0.)) THEN value = Dot (uvec, mp_right_uvec) upper_limit = mp_right_maxdot ELSE ! right limit doesn't apply; give a "free pass" value = 0. upper_limit = 1. END IF CASE(4) ! bottom boundary IF (mp_using_bottom.AND.(Dot(uvec, mp_bottom_applies_uvec) > 0.)) THEN value = Dot (uvec, mp_bottom_uvec) upper_limit = mp_bottom_maxdot ELSE ! bottom limit doesn't apply; give a "free pass" value = 0. upper_limit = 1. END IF END SELECT IF (value > upper_limit) THEN In_Rind = 'O' ! outside RETURN ELSE IF (value == upper_limit) THEN In_Rind = 'B' ! ELSE ! value < upper_limit ! In_Rind remains as 'I' or 'B' (whichever it was). END IF END DO END FUNCTION In_Rind LOGICAL FUNCTION L5_In_Window (uvec) !Tests whether a point on the planet, !represented by level-5 Cartesian unit vector "uvec", !will be in sight in the final map window. !Note that use of this function is strictly OPTIONAL; !it is normal to just plot features all over the planet and !to let the automatic processing of paths and text determine !what is visible. Also, calling this routine for points !that are likely to be visible wastes time, !because the map projection and windowing will be done all !over again when something is actually plotted. !So, calling this function in advance is ONLY useful if !it will take very extensive computation to create the thing !being plotted, and the programmer wants to avoid doing that !computation for invisible objects, which are expected to !be very numerous. IMPLICIT NONE REAL, DIMENSION(3), INTENT(IN) :: uvec CHARACTER*1 :: c1 REAL :: x_meters, x_points, y_meters, y_points !------------------------- c1 = In_Rind(uvec) IF (c1 == 'O') THEN ! outside rind L5_In_Window = .FALSE. ELSE CALL Project (uvec = uvec, x = x_meters, y = y_meters) CALL Meters_2_Points (x_meters,y_meters, x_points,y_points) c1 = In_Window (x_points, y_points) IF (c1 == 'O') THEN ! outside window L5_In_Window = .FALSE. ELSE L5_In_Window = .TRUE. END IF END IF END FUNCTION L5_In_Window 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 Meters_2_Points (x_meters,y_meters, x_points,y_points) ! Takes one point from level 3 (projection plane, in meters) ! to level 1 or 2 (page, in points) IMPLICIT NONE REAL, INTENT(IN) :: x_meters, y_meters REAL, INTENT(OUT) :: x_points, y_points x_points = mp_m_2_pts_vector(1) + & & mp_m_2_pts_matrix(1,1) * x_meters + & & mp_m_2_pts_matrix(1,2) * y_meters y_points = mp_m_2_pts_vector(2) + & & mp_m_2_pts_matrix(2,1) * x_meters + & & mp_m_2_pts_matrix(2,2) * y_meters END SUBROUTINE Meters_2_Points 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 SUBROUTINE Points_2_LonLat (x_points,y_points, success, lon,lat) ! For a given position on the page (measured from lower left), ! attempts to find a corresponding (longitude, latitude) ! according to the current zoom and map projection. ! Results "lon" and "lat" are only useful when "success" = T. ! This routine is often used to locate degree ticks along the ! edge of the map window. Therefore, points are allowed ! to be outside the window. IMPLICIT NONE REAL, INTENT(IN) :: x_points, y_points LOGICAL, INTENT(OUT) :: success REAL, INTENT(OUT) :: lon, lat REAL :: x_meters, y_meters REAL, DIMENSION(3) :: uvec CALL Points_2_Meters (x_points,y_points, x_meters,y_meters) CALL Reject (x_meters,y_meters, success, uvec) IF (success) CALL Uvec_2_LonLat (uvec, lon, lat) END SUBROUTINE Points_2_LonLat SUBROUTINE Points_2_Meters (x_points,y_points, x_meters,y_meters) ! Takes one point from level 1 or 2 (page, in points) ! to level 3 (projection plane, in meters) IMPLICIT NONE REAL, INTENT(IN) :: x_points, y_points REAL, INTENT(OUT) :: x_meters, y_meters x_meters = mp_pts_2_m_vector(1) + & & mp_pts_2_m_matrix(1,1) * x_points + & & mp_pts_2_m_matrix(1,2) * y_points y_meters = mp_pts_2_m_vector(2) + & & mp_pts_2_m_matrix(2,1) * x_points + & & mp_pts_2_m_matrix(2,2) * y_points END SUBROUTINE Points_2_Meters SUBROUTINE Process_L3_Paths () ! Perform an in-place translation from meters to points, ! and pass on to the next level. ! Another action performed at this point is to replace all ! curved segments which are exactly or nearly straight ! with straight-line segments. ! This action will catch all segments from higher levels, too. ! (This is necessary because in Module Adobe_Illustrator, ! the test was applied in Curve_To_L12, and these segments ! do not pass throught that user interface.) IMPLICIT NONE INTEGER :: i, path LOGICAL inline1, inline2 REAL :: crossx, crossy, dot, & & last_x_points, last_y_points, length, offline, & & ux, uy, vx, vy, & & x0_points, x1_points, x2_points, x3_points, & & y0_points, y1_points, y2_points, y3_points DO WHILE (ai_Ln_paths(3) > 0) path = Next_Path(level = 3) CALL Meters_2_Points(ai_pathlib(1,0,path),ai_pathlib(2,0,path),x0_points,y0_points) ai_pathlib(1,0,path) = x0_points ai_pathlib(2,0,path) = y0_points last_x_points = x0_points last_y_points = y0_points DO i = 1, ai_segments(path) CALL Meters_2_Points(ai_pathlib(1,i,path),ai_pathlib(2,i,path),x1_points,y1_points) ai_pathlib(1,i,path) = x1_points ai_pathlib(2,i,path) = y1_points IF (ai_bent(i,path)) THEN CALL Meters_2_Points(ai_pathlib(3,i,path),ai_pathlib(4,i,path),x2_points,y2_points) ai_pathlib(3,i,path) = x2_points ai_pathlib(4,i,path) = y2_points CALL Meters_2_Points(ai_pathlib(5,i,path),ai_pathlib(6,i,path),x3_points,y3_points) ai_pathlib(5,i,path) = x3_points ai_pathlib(6,i,path) = y3_points ! Simplify degenerate curveto's !(those which are straight within a tolerance ! of 1 point) by converting them to lineto's. !(This will also cover the case of handle points ! that erroneously coincide with end points.) IF ((x3_points /= last_x_points).OR. & &(y3_points /= last_y_points)) THEN ! straight line is defined; test can proceed: ux = x3_points - last_x_points uy = y3_points - last_y_points length = SQRT(ux**2 + uy**2) ! > 0. ux = ux / length ! unit vector of baseline uy = uy / length vx = x1_points - last_x_points ! lever1 vector vy = y1_points - last_y_points dot = vx * ux + vy * uy ! length of // component crossx = vx - dot * ux ! perpendicular component crossy = vy - dot * uy offline = SQRT(crossx**2 + crossy**2) inline1 = (offline <= ai_min_amplitude_points) vx = x3_points - x2_points ! -lever2 vector vy = y3_points - y2_points dot = vx * ux + vy * uy ! length of // component crossx = vx - dot * ux ! perpendicular component crossy = vy - dot * uy offline = SQRT(crossx**2 + crossy**2) inline2 = (offline <= ai_min_amplitude_points) IF (inline1.AND.inline2) THEN ! replace curved segment with straight ai_bent(i,path) = .FALSE. ai_pathlib(1,i,path) = x3_points ai_pathlib(2,i,path) = y3_points END IF END IF ! in-line test was possible last_x_points = x3_points last_y_points = y3_points ELSE ! segment presented as straight last_x_points = x1_points last_y_points = y1_points END IF ! segment presented as bent / straight END DO ! segments in path ai_Ln_paths(3) = ai_Ln_paths(3) - 1 ai_path_level(path) = 2 ai_Ln_paths(2) = ai_Ln_paths(2) + 1 CALL Process_L2_Paths END DO END SUBROUTINE Process_L3_Paths SUBROUTINE Process_L3_Text (x_meters, y_meters, page_angle_radians, & & font_points, lr_fraction, ud_fraction, & & text) ! Apply a simple meters -> points position conversion and pass on. IMPLICIT NONE INTEGER, INTENT(IN) :: font_points REAL, INTENT(IN) :: x_meters, y_meters, page_angle_radians, & & lr_fraction, ud_fraction CHARACTER*(*),INTENT(IN) :: text REAL :: x_points, y_points CALL Meters_2_Points (x_meters,y_meters, x_points,y_points) CALL Process_L2_Text (x_points, y_points, page_angle_radians, & & font_points, lr_fraction, ud_fraction, & & text) END SUBROUTINE Process_L3_Text SUBROUTINE Process_L4_Paths () ! Takes path(s) on the sphere which has (have) already been cut/limited ! and projects it (them) to the [x,y] system of the (planet-sized) ! projection plane. Then, pass it (them) on to Process_L3_Paths. ! There are 3 steps: ! (1) Limit the curvature in any one segment to 90 degrees left/right ! (like circles in Adobe Illustrator) and limit arc-length to ! about 10 degrees (because projection distortion is non-uniform). ! (2) Find equivalents of Bezier handle points on sphere. ! (3) Project Bezier handle points and end-points to plane. IMPLICIT NONE INTEGER :: i, j, path, pieces REAL, PARAMETER :: piece_size = 0.174533 ! 10 degrees, in radians REAL :: azim_0, azim_1, azim_2, azim_3, azim_a, azim_b, azimuth, & & arc_length, crossed, curviness, dotted, forward, & & radians, small_size_radians, x0_meters, y0_meters, & & x1_meters, y1_meters, x2_meters, y2_meters, & & x3_meters, y3_meters REAL, DIMENSION(3) :: ccw_pole_uvec, goal_uvec, guide_uvec, last_piece_uvec, & & last_uvec, omega_uvec, result_uvec, t_vec DO WHILE (ai_Ln_paths(4) > 0) path = Next_Path(level = 4) last_uvec = ai_pathlib(1:3,0,path) ! initialize memory CALL Project(uvec = last_uvec, x = x0_meters, y = y0_meters) CALL New_L3_Path(x0_meters, y0_meters) ! start L3 path DO i = 1, ai_segments(path) !------------------------------------------------ ! decide how many pieces (at least 1) IF (ai_bent(i,path)) THEN ! small-circle arc ccw_pole_uvec = ai_pathlib(1:3,i,path) goal_uvec = ai_pathlib(4:6,i,path) small_size_radians = Arc(ccw_pole_uvec, goal_uvec) IF ((goal_uvec(1) == last_uvec(1)).AND. & & (goal_uvec(2) == last_uvec(2)).AND. & & (goal_uvec(3) == last_uvec(3))) THEN ! End-point = start-point; means complete circle radians = Two_Pi pieces = 4 ! no more than 90 deg. turning each azim_a = Relative_Compass(ccw_pole_uvec, last_uvec) azim_b = azim_a - Two_Pi ELSE ! less than a complete circle azim_a = Relative_Compass(ccw_pole_uvec, last_uvec) azim_b = Relative_Compass(ccw_pole_uvec, goal_uvec) IF (azim_b > azim_a) azim_b = azim_b - Two_Pi ! ccw radians = azim_a - azim_b ! > 0. curviness = ABS( COS(small_size_radians) ) pieces = MAX (1, IAbove(curviness*radians/Pi_over_2) ) END IF ! complete small circle ! Also limit to no more than piece_size arclength. arc_length = radians * SIN(small_size_radians) pieces = MAX(pieces, IAbove(arc_length/piece_size) ) ELSE ! great-circle arc small_size_radians = Pi_over_2 ! find pole of great circle, and rotation angle goal_uvec = ai_pathlib(1:3,i,path) CALL Cross(last_uvec, goal_uvec, t_vec) CALL Make_Uvec(t_vec, ccw_pole_uvec) crossed = Length(t_vec) ! > 0. dotted = Dot(last_uvec, goal_uvec) radians = ATAN2(crossed, dotted) ! > 0. azim_a= Relative_Compass (ccw_pole_uvec, last_uvec) azim_b = azim_a - radians arc_length = radians pieces = IAbove(radians/piece_size) END IF ! small-circle / great-circle segment !------------------------------------------------ ! find equivalent Bezier control points on sphere, ! project, and write immediately to L3 path azim_0 = azim_a last_piece_uvec = last_uvec DO j = 1, pieces IF (j == pieces) THEN azim_3 = azim_b result_uvec = goal_uvec ELSE ! one of several pieces, going ccw azim_3 = azim_a - (j*radians)/pieces CALL Turn_To (azim_3, ccw_pole_uvec, & & small_size_radians, & ! inputs & omega_uvec, result_uvec) ! results END IF ! only piece, or one-of-many pieces CALL Project(uvec = result_uvec, x = x3_meters, y = y3_meters) guide_uvec = result_uvec ! enforces correct side of cut! IF (ai_bent(i,path)) THEN ! piece is a small-circle arc forward = 0.3516 * arc_length / pieces ! Bezier handle on initial end of piece: azimuth = Pi_over_2 + Relative_Compass (last_piece_uvec, ccw_pole_uvec) CALL Turn_To (azimuth, last_piece_uvec, & & forward, & ! inputs & omega_uvec, result_uvec) ! outputs CALL Project (uvec = result_uvec, guide = guide_uvec, & & x = x1_meters, y = y1_meters) ! Bezier handle on final end of piece: azimuth = Relative_Compass (guide_uvec, ccw_pole_uvec) - Pi_over_2 CALL Turn_To (azimuth, guide_uvec, & & forward, & ! inputs & omega_uvec, result_uvec) ! outputs CALL Project (uvec = result_uvec, guide = guide_uvec, & & x = x2_meters, y = y2_meters) ELSE ! piece is a great circle arc ! Bezier handle on initial end of piece: azim_1 = 0.6484*azim_0 + 0.3516*azim_3 CALL Turn_To (azim_1, ccw_pole_uvec, & & small_size_radians, & ! inputs & omega_uvec, result_uvec) ! results CALL Project (uvec = result_uvec, guide = guide_uvec, & & x = x1_meters, y = y1_meters) ! Bezier handle on final end of piece: azim_2 = 0.3516*azim_0 + 0.6484*azim_3 CALL Turn_To (azim_2, ccw_pole_uvec, & & small_size_radians, & ! inputs & omega_uvec, result_uvec) ! results CALL Project (uvec = result_uvec, guide = guide_uvec, & & x = x2_meters, y = y2_meters) END IF CALL Curve_To_L3 (x1_meters,y1_meters, & & x2_meters,y2_meters, & & x3_meters,y3_meters) azim_0 = azim_3 ! set memory for next j last_piece_uvec = guide_uvec ! remember old final point END DO ! j = 1, pieces last_uvec = goal_uvec ! update memory for loop i END DO ! i = 1, segments (of old L4 path) CALL End_L3_Path(close = ai_closed(path), & & stroke = ai_stroked(path), fill = ai_filled(path)) !- - - - - - - - - - - - - - - - - - - - - - - - ! free original L4 path ai_path_level(path) = 0 ai_Ln_paths(4) = ai_Ln_paths(4) - 1 ai_total_paths = ai_total_paths - 1 !- - - - - - - - - - - - - - - - - - - - - - - - CALL Process_L3_Paths !- - - - - - - - - - - - - - - - - - - - - - - - END DO ! on waiting L4 paths END SUBROUTINE Process_L4_Paths SUBROUTINE Process_L4_Text (uvec, angle_radians, from_east, & & font_points, lr_fraction, ud_fraction, & & text) ! Project position from sphere (already checked as In_Rind) to ! projection plane (x,y) in meters, and pass on. IMPLICIT NONE REAL, DIMENSION(3), INTENT(IN) :: uvec REAL, INTENT(IN) :: angle_radians, & & lr_fraction, ud_fraction LOGICAL, INTENT(IN) :: from_east INTEGER, INTENT(IN) :: font_points CHARACTER*(*), INTENT(IN) :: text REAL :: equat, lon_radians, page_angle_radians, x_meters, x_points, xe_meters, & & xe_points, y_meters, y_points, ye_meters, ye_points REAL, DIMENSION(3) :: east_post CALL Project (uvec = uvec, x = x_meters, y = y_meters) IF (from_east) THEN IF (ABS(uvec(3)) >= 1.) THEN ! Cannot measure from_east at N or S pole. Pretend not requested. page_angle_radians = angle_radians ELSE ! normal case: ! create a temporary vector position "east_post" east of uvec: equat = SQRT(uvec(1)**2 + uvec(2)**2) lon_radians = ATAN2(uvec(2),uvec(1)) east_post(1) = equat * COS(lon_radians + radians_per_degree) east_post(2) = equat * SIN(lon_radians + radians_per_degree) east_post(3) = uvec(3) ! project both points to the (x,y) system: CALL Project (uvec = east_post, guide = uvec, x = xe_meters, y = ye_meters) ! map both points onto the page: CALL Meters_2_Points (x_meters,y_meters, x_points,y_points) CALL Meters_2_Points (xe_meters,ye_meters, xe_points,ye_points) page_angle_radians = ATAN2((ye_points - y_points), (xe_points - x_points)) & & + angle_radians END IF ELSE ! from_east feature not requested; measure from horizontal on page. page_angle_radians = angle_radians END IF CALL Process_L3_Text (x_meters,y_meters, page_angle_radians, & & font_points, lr_fraction, ud_fraction, & & text) END SUBROUTINE Process_L4_Text SUBROUTINE Process_L5_Paths () ! Checks Level-5 paths (on the unit sphere, composed of arcs ! of great and small circles, and expressed in terms of ! Cartesian unit vectors (uvec's) from the center of the sphere. ! If they go outside any of the 4 predefined limits ! (all of which are great or small circles or half-circles), ! they are truncated. Complex logic for filled paths adds ! segments as needed along the boundaries to keep them closed. ! At the end, the path(s) are passed to Process_L4_Paths for ! subdivision into small pieces, conversion to Bezier curves, ! and projection onto the projection plane. IMPLICIT NONE CHARACTER*1 :: abyte CHARACTER*1, DIMENSION(1:ai_longest) :: seg_memo ! O, B, I CHARACTER*1, DIMENSION(0:ai_longest) :: point_memo ! O, B, I DOUBLE PRECISION :: dot_a, dot_b INTEGER, PARAMETER :: inout_list_limit = 40 INTEGER :: first, i, intersections, inout_list_count, & & ip1, j, k, & & kp1, last, n_area, new_path, & & new_segments, number, path, trial_first, trial_last LOGICAL :: all_segs_in, all_segs_out, area_done, & & clockwise, complete_circle, connects, curve_shows, & & filled, found_in_list, & & last_shown, next_is_out, null_arc, overlap, & & point0_shows, point1_shows, seg_shows, this_is_out, this_shown, unchanged LOGICAL, DIMENSION(ai_longest) :: new_segment_used REAL :: azimuth, base_of_t, far_radians, from_t_wall, phi, & & phi_end, phi_start, phi0, phi1, t, t_wall, to_t_wall, zeta, zeta0, zeta1 REAL, DIMENSION(3) :: closure_uvec, first_a_uvec, first_b_uvec, from_point_uvec, & & goal_uvec, last_uvec, last_a_uvec, last_b_uvec, omega_uvec, & & point1_uvec, point2_uvec, pole_uvec, pole_a_uvec, pole_b_uvec, & & start_uvec, t_vec, t_uvec, to_point_uvec REAL, DIMENSION(16) :: t_edge, t15_edge, a_edge, b_edge, g_edge REAL, DIMENSION(inout_list_limit) :: t_1_to_5, & & a_inout_list, b_inout_list, g_inout_list REAL, DIMENSION(1:ai_longest) :: t15_memo DO WHILE (ai_Ln_paths(5) > 0) path = Next_Path(level = 5) filled = ai_filled(path) ! Copy path to a new slot at level -1 (because it may expand), ! adding new control points WITHIN SEGMENTS where they cross ! the window boundary. (If a segment coincides with the ! boundary for some distance, both initial and final points ! will be in the new list, either as original or added points.) ! While working, record whether each point (either an original ! control point, or an added control point) is outside ('O'), ! on the boundary ('B'), or inside ('I'). new_path = Next_Free_Path() ai_path_level(new_path) = -1 ! temporary working array, not to be "Process"ed! ai_total_paths = ai_total_paths + 1 ai_closed(new_path) = ai_closed(path) ai_stroked(new_path) = ai_stroked(path) ai_filled(new_path) = ai_filled(path) ai_pathlib(1:3,0,new_path) = ai_pathlib(1:3,0,path) last_uvec = ai_pathlib(1:3,0,path) point_memo(0) = In_Rind (last_uvec) unchanged = (point_memo(0) == 'I').OR.(point_memo(0) == 'B') new_segments = 0 DO i = 1, ai_segments(path) IF (ai_bent(i,path)) THEN goal_uvec = ai_pathlib(4:6,i,path) pole_uvec = ai_pathlib(1:3,i,path) far_radians = Arc (pole_uvec, goal_uvec) complete_circle = (goal_uvec(1) == last_uvec(1)).AND. & & (goal_uvec(2) == last_uvec(2)).AND. & & (goal_uvec(3) == last_uvec(3)) null_arc = .FALSE. ELSE ! great circle arc goal_uvec = ai_pathlib(1:3,i,path) CALL Cross (last_uvec, goal_uvec, t_vec) null_arc = (MAX(ABS(t_vec(1)),ABS(t_vec(2)),ABS(t_vec(3))) == 0.0) IF (null_arc) THEN ! pole is undefined IF (ABS(goal_uvec(3)) > 0.999) THEN t_uvec = (/ 1.0, 0.0, 0.0 /) ELSE t_uvec = (/ 0.0, 0.0, 1.0 /) END IF CALL Cross (t_uvec, goal_uvec, t_vec) !result will be perpendicular to goal_uvec END IF CALL Make_Uvec (t_vec, pole_uvec) far_radians = Pi_over_2 complete_circle = .FALSE. END IF abyte = In_Rind (goal_uvec) ! preserve until end!!! IF ((point_memo(i-1) == 'O').OR.(abyte == 'O')) unchanged = .FALSE. ! check for any intersections (0...2 each) with each side of rind pole_a_uvec = pole_uvec dot_a = COS(1.0D0 * Arc(pole_uvec, goal_uvec)) !Note: This is logically equivalent to: ! dot_a = Dot (pole_uvec, goal_uvec) !but is more accurate when the two vectors are not PRECISELY unit-vectors. first_a_uvec = last_uvec last_a_uvec = goal_uvec ! set up internal coordinate phi = phi0 to phi1 for segment: phi0 = -Relative_Compass (pole_uvec, last_uvec) phi1 = -Relative_Compass (pole_uvec, goal_uvec) IF (phi1 < phi0) THEN phi1 = phi1 + Two_Pi ELSE IF (phi1 == phi0) THEN IF (complete_circle) phi1 = phi1+ Two_Pi END IF intersections = 0 four_walls: DO j = 1, 4 ! 4 sides of rind, going around counterclockwise SELECT CASE(j) CASE(1) ! bottom IF (.NOT.mp_using_bottom) CYCLE four_walls pole_b_uvec = mp_bottom_uvec dot_b = mp_bottom_maxdot first_b_uvec = mp_first_bottom_uvec last_b_uvec = mp_last_bottom_uvec base_of_t = 1. CASE(2) ! right IF (.NOT.mp_using_right) CYCLE four_walls pole_b_uvec = mp_right_uvec dot_b = mp_right_maxdot first_b_uvec = mp_first_right_uvec last_b_uvec = mp_last_right_uvec base_of_t = 2. CASE(3) ! top IF (.NOT.mp_using_top) CYCLE four_walls pole_b_uvec = mp_top_uvec dot_b = mp_top_maxdot first_b_uvec = mp_first_top_uvec last_b_uvec = mp_last_top_uvec base_of_t = 3. CASE(4) ! left IF (.NOT.mp_using_left) CYCLE four_walls pole_b_uvec = mp_left_uvec dot_b = mp_left_maxdot first_b_uvec = mp_first_left_uvec last_b_uvec = mp_last_left_uvec base_of_t = 4. END SELECT CALL 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 IF (number >= 1) THEN unchanged = .FALSE. intersections = intersections + 1 phi = -Relative_Compass (pole_uvec, point1_uvec) IF (phi < phi0) phi = phi + Two_Pi t = (phi - phi0) / (phi1 - phi0) ! 0. to 1. t_edge(intersections) = t ! find internal coordinate of intersection point in wall ! (considering that internal coordinate of rind boundary ! winds in the opposite way from first->last in each wall): zeta0 = Relative_Compass (pole_b_uvec, last_b_uvec) zeta1 = Relative_Compass (pole_b_uvec, first_b_uvec) IF (zeta1 <= zeta0) zeta1 = zeta1 + Two_Pi zeta = Relative_Compass (pole_b_uvec, point1_uvec) IF (zeta < zeta0) zeta = zeta + Two_Pi t_wall = (zeta - zeta0) / (zeta1 - zeta0) t_wall = MAX(0.,MIN(t_wall,0.99999)) t15_edge(intersections) = base_of_t + t_wall a_edge(intersections) = point1_uvec(1) b_edge(intersections) = point1_uvec(2) g_edge(intersections) = point1_uvec(3) IF (number == 2) THEN intersections = intersections + 1 ! internal coordinate of intersection point in segment phi = -Relative_Compass (pole_uvec, point2_uvec) IF (phi < phi0) phi = phi + Two_Pi t = (phi - phi0) / (phi1 - phi0) t_edge(intersections) = t ! find internal coordinate of intersection point in wall: zeta = Relative_Compass (pole_b_uvec, point2_uvec) IF (zeta < zeta0) zeta = zeta + Two_Pi t_wall = (zeta - zeta0) / (zeta1 - zeta0) t_wall = MAX(0.,MIN(t_wall,0.99999)) t15_edge(intersections) = base_of_t + t_wall a_edge(intersections) = point2_uvec(1) b_edge(intersections) = point2_uvec(2) g_edge(intersections) = point2_uvec(3) END IF ! 2 intersections found with this side END IF ! any intersections found with this side END DO four_walls ! check all four sides of rind for intersections !sort the list of intersections (up to 8) by t (of segment) IF (intersections > 1) CALL Sort_Lists(intersections, t_edge, t15_edge, a_edge, b_edge, g_edge) IF (intersections > 0) THEN phi_start = phi0 DO j = 1, intersections new_segments = new_segments + 1 IF (new_segments >= ai_longest) THEN WRITE (*,"(' ERROR: Path length > ai_longest =',I6)") ai_longest CALL Traceback END IF IF (ai_bent(i,path)) THEN ai_bent(new_segments,new_path) = .TRUE. ai_pathlib(1:3,new_segments,new_path) = pole_uvec ai_pathlib(4,new_segments,new_path) = a_edge(j) ai_pathlib(5,new_segments,new_path) = b_edge(j) ai_pathlib(6,new_segments,new_path) = g_edge(j) ELSE ai_bent(new_segments,new_path) = .FALSE. ai_pathlib(1,new_segments,new_path) = a_edge(j) ai_pathlib(2,new_segments,new_path) = b_edge(j) ai_pathlib(3,new_segments,new_path) = g_edge(j) END IF ! find midpoint of segment to classify it: IF (phi1 == phi0) THEN ! avoid: NaN = Infinity * 0 phi_end = phi0 ELSE phi_end = phi0 + t_edge(j)*(phi1 - phi0) END IF azimuth = -0.5*(phi_start + phi_end) CALL Turn_To (azimuth, pole_uvec, far_radians, & ! inputs & omega_uvec, t_uvec) seg_memo(new_segments) = In_Rind (t_uvec) point_memo(new_segments) = 'B' ! border point t15_memo(new_segments) = t15_edge(j) phi_start = phi_end ! memory, for next j (or completion) END DO ! j = 1, intersections IF (t_edge(intersections) < 1.0) THEN ! complete the segment new_segments = new_segments + 1 IF (new_segments >= ai_longest) THEN WRITE (*,"(' ERROR: Path length > ai_longest =',I6)") ai_longest CALL Traceback END IF IF (ai_bent(i,path)) THEN ai_bent(new_segments,new_path) = .TRUE. ai_pathlib(1:3,new_segments,new_path) = pole_uvec ai_pathlib(4:6,new_segments,new_path) = goal_uvec ELSE ai_bent(new_segments,new_path) = .FALSE. ai_pathlib(1:3,new_segments,new_path) = goal_uvec END IF ! find midpoint of segment to classify it: azimuth = -0.5*(phi_start + phi1) CALL Turn_To (azimuth, pole_uvec, far_radians, & ! inputs & omega_uvec, t_uvec) seg_memo(new_segments) = In_Rind (t_uvec) point_memo(new_segments) = abyte END IF ! divided segment needs completion ELSE ! no intersections found for this segment new_segments = new_segments + 1 IF (new_segments >= ai_longest) THEN WRITE (*,"(' ERROR: Path length > ai_longest =',I6)") ai_longest CALL Traceback END IF IF (ai_bent(i,path)) THEN ai_bent(new_segments,new_path) = .TRUE. ai_pathlib(1:3,new_segments,new_path) = pole_uvec ai_pathlib(4:6,new_segments,new_path) = goal_uvec ELSE ai_bent(new_segments,new_path) = .FALSE. ai_pathlib(1:3,new_segments,new_path) = goal_uvec END IF ! find midpoint of segment to classify it: azimuth = -0.5*(phi0 + phi1) CALL Turn_To (azimuth, pole_uvec, far_radians, & ! inputs & omega_uvec, t_uvec) seg_memo(new_segments) = In_Rind (t_uvec) point_memo(new_segments) = abyte END IF ! intersections found / not-found for this segment ! set memory with initial point for next segment: last_uvec = goal_uvec END DO ! i = 1, ai_segments(path) IF (filled) THEN ! Filled (and maybe stroked?) paths may require ! completions along the window boundary!!! ! First, decide whether path is all-in, all-out, ! or some complex mixture: all_segs_in = .TRUE. all_segs_out = .TRUE. DO i = 1, new_segments IF (seg_memo(i) == 'O') all_segs_in = .FALSE. IF ((seg_memo(i) == 'I').OR.(seg_memo(i) == 'B')) all_segs_out = .FALSE. END DO IF (all_segs_in) THEN ! No need to modify path at all! Just change level. ai_Ln_paths(5) = ai_Ln_paths(5) - 1 ai_path_level(path) = 4 ai_Ln_paths(4) = ai_Ln_paths(4) + 1 CALL Process_L4_Paths ELSE IF (all_segs_out) THEN ! No action. While it is technically possible to ! have a path that lies completely outside the rind ! but "surrounds" it, we ignore this possibility ! in practice because the rind is so big ! (almost the whole planet). ELSE ! filled path crosses boundary; complex logic! ! Does path circulate clockwise or counterclockwise??? clockwise = Clockways(path) ! Build a short list of boundary-crossing points, and ! then sort them into order that will be needed for interpolation: inout_list_count = 0 DO i = 1, new_segments ip1 = MOD(i,new_segments) + 1 this_is_out = (seg_memo(i) == 'O') next_is_out = (seg_memo(ip1) == 'O') IF (this_is_out .NEQV. next_is_out) THEN inout_list_count = inout_list_count + 1 IF (inout_list_count > inout_list_limit) THEN WRITE (*,"(' ERROR: Increase parameter inout_list_limit.')") CALL Traceback END IF t_1_to_5(inout_list_count) = t15_memo(i) IF (ai_bent(i,new_path)) THEN ! small-circle segment a_inout_list(inout_list_count) = ai_pathlib(4,i,new_path) b_inout_list(inout_list_count) = ai_pathlib(5,i,new_path) g_inout_list(inout_list_count) = ai_pathlib(6,i,new_path) ELSE ! great-circle segment a_inout_list(inout_list_count) = ai_pathlib(1,i,new_path) b_inout_list(inout_list_count) = ai_pathlib(2,i,new_path) g_inout_list(inout_list_count) = ai_pathlib(3,i,new_path) END IF END IF ! got another inout point at end of segment END DO ! i = 1, new_segments IF (clockwise) THEN ! reverse definition of t to get points in clockwise order: DO j = 1, inout_list_count t_1_to_5(j) = -t_1_to_5(j) END DO END IF IF (inout_list_count > 2) CALL Sort_Lists(inout_list_count, & & t_1_to_5, a_inout_list, b_inout_list, g_inout_list) IF (clockwise) THEN ! undo the redefinition of t; get it back! DO j = 1, inout_list_count t_1_to_5(j) = -t_1_to_5(j) END DO END IF DO i = 1, new_segments new_segment_used(i) = (seg_memo(i) == 'O') END DO all_areas: DO n_area = 1, 4 ! although probably there are only 1 or 2 ! Find a unused starting segment which is in, but follows ! a segment which is out. Thus, its initial point will ! be a loose end on the boundary that ultimately needs to ! be the end point of the last segment in a path. first = 0 ! if not changed, will show failure find_first: DO i = 1, new_segments trial_first = i IF (.NOT.new_segment_used(trial_first)) THEN trial_last = MOD((i - 2 + new_segments),new_segments) + 1 IF (((seg_memo(trial_first)=='I').OR.(seg_memo(trial_first)=='B')) & & .AND.(seg_memo(trial_last)=='O')) THEN first = trial_first last = trial_last EXIT find_first END IF ! this is visible, previous was not END IF ! trial_first is still available END DO find_first IF (first == 0) EXIT all_areas ! all visible segments used last_shown = .FALSE. DO i = 1, new_segments ! but j will be true loop variable j = MOD((first + i - 2), new_segments) + 1 ! j = first...last seg_shows = ((seg_memo(j) == 'I').OR.(seg_memo(j) == 'B')) IF (seg_shows) THEN IF (ai_in_path) THEN ! Does this segment continue from where we are? IF (j == 1) THEN start_uvec(1:3) = ai_pathlib(1:3,0,new_path) ELSE ! j > 1 IF (ai_bent(j-1,new_path)) THEN ! previous was small circle start_uvec(1:3) = ai_pathlib(4:6,j-1,new_path) ELSE ! previous segment was a great circle arc start_uvec(1:3) = ai_pathlib(1:3,j-1,new_path) END IF ! previous seg bent / straight END IF ! j == 1 or greater connects = last_shown .OR. & ! simple sequence or good match & ((start_uvec(1) == mp_last_uvec(1)).AND. & & (start_uvec(2) == mp_last_uvec(2)).AND. & & (start_uvec(3) == mp_last_uvec(3))) this_shown = connects ! set switch for code below ELSE ! not ai_in_path; start one this_shown = .TRUE. ! Need to start new path (and remember the closure point, ! which should always be a boundary point) IF (j == 1) THEN closure_uvec(1:3) = ai_pathlib(1:3,0,new_path) ELSE ! j > 1 IF (ai_bent(j-1,new_path)) THEN ! previous was small circle closure_uvec(1:3) = ai_pathlib(4:6,j-1,new_path) ELSE ! previous segment was a great circle arc closure_uvec(1:3) = ai_pathlib(1:3,j-1,new_path) END IF ! previous seg bent / straight END IF ! j == 1 or greater CALL New_L45_Path(4, closure_uvec) END IF ! ai_in_path, or not IF (this_shown) THEN IF (ai_bent(j,new_path)) THEN ! small pole_uvec = ai_pathlib(1:3,j,new_path) goal_uvec = ai_pathlib(4:6,j,new_path) CALL Small_To_L45 (pole_uvec, goal_uvec) new_segment_used(j) = .TRUE. ELSE ! segment j is great goal_uvec = ai_pathlib(1:3,j,new_path) CALL Great_To_L45 (goal_uvec) new_segment_used(j) = .TRUE. END IF ! segment j is small/great circle arc END IF ! this_shown ELSE ! segment j is not shown this_shown = .FALSE. IF (last_shown) THEN ! Need to complete (along boundary) and close segment!!! ! Starting point of boundary loop is from_point_uvec: IF (j == 1) THEN from_point_uvec = ai_pathlib(1:3,0,new_path) ELSE ! j > 1 IF (ai_bent(j-1,new_path)) THEN from_point_uvec = ai_pathlib(4:6,j-1,new_path) ELSE ! previous segment great circle arc from_point_uvec = ai_pathlib(1:3,j-1,new_path) END IF ! previous seg small / great END IF ! j == 1 or greater ! Ending point of boundary loop (to_point_uvec) may be ! either (closure_uvec) where this path started, ! or else a different boundary point (where path re-enters): ! Use the ordered list of boundary points to find the next ! re-entry in the "clockwise"(?) direction: found_in_list = .FALSE. lookup: DO k = 1, inout_list_count IF ((from_point_uvec(1) == a_inout_list(k)).AND. & &(from_point_uvec(2) == b_inout_list(k)).AND. & &(from_point_uvec(3) == g_inout_list(k))) THEN ! found from_point_uvec in list found_in_list = .TRUE. ! next point (cyclically) in list: kp1 = MOD(k,inout_list_count) + 1 to_point_uvec(1) = a_inout_list(kp1) to_point_uvec(2) = b_inout_list(kp1) to_point_uvec(3) = g_inout_list(kp1) ! t15 ccordinates of each point: from_t_wall = t_1_to_5(k) to_t_wall = t_1_to_5(kp1) EXIT lookup END IF ! found current from_point_uvec in list END DO lookup ! k = 1, inout_list_count IF (.NOT.found_in_list) THEN WRITE (*,"(' ERROR: Boundary point not found in ordered list.')") CALL Traceback END IF CALL Trace_Rind(from_t_wall, clockwise,to_point_uvec,to_t_wall) area_done = (to_point_uvec(1) == closure_uvec(1)).AND. & &(to_point_uvec(2) == closure_uvec(2)).AND. & &(to_point_uvec(3) == closure_uvec(3)) IF (area_done) THEN CALL End_L45_Path(close = .TRUE., stroke = ai_stroked(path), fill = .TRUE.) CYCLE all_areas END IF ! area_done END IF ! last_shown END IF ! segment j is shown / not shown last_shown = this_shown END DO ! on all segments in filled, boundary-cutting path END DO all_areas ! multiple passes, in case of disconnected areas !All areas should now be closed, but check this! IF (ai_in_path) THEN WRITE (*,"(' ERROR: Filled path did not close as expected.')") CALL Traceback END IF END IF ! path is all-in / all-out / crosses rind-wall boundary ELSE ! path is only stroked last_uvec = ai_pathlib(1:3,0,new_path) point0_shows = (point_memo(0) == 'I').OR.(point_memo(0) == 'B') IF (point0_shows) CALL New_L45_Path(4,last_uvec) DO i = 1, new_segments point1_shows = (point_memo(i) == 'I').OR.(point_memo(i) == 'B') IF ((point_memo(i-1) == 'I').OR.(point_memo(i) == 'I')) THEN curve_shows = .TRUE. ELSE IF ((point_memo(i-1) == 'O').OR.(point_memo(i) == 'O')) THEN curve_shows = .FALSE. ELSE ! both control points are on the boundary; ! which way does the line bulge? curve_shows = (seg_memo(i) == 'I').OR.(seg_memo(i) == 'B') END IF IF (curve_shows) THEN IF (ai_in_path) THEN IF (ai_bent(i,new_path)) THEN pole_uvec = ai_pathlib(1:3,i,new_path) t_uvec = ai_pathlib(4:6,i,new_path) CALL Small_To_L45 (pole_uvec, t_uvec) ELSE ! great circle arc t_uvec = ai_pathlib(1:3,i,new_path) CALL Great_To_L45 (t_uvec) END IF ELSE ! not in path CALL New_L45_Path (4,last_uvec) IF (ai_bent(i,new_path)) THEN pole_uvec = ai_pathlib(1:3,i,new_path) t_uvec = ai_pathlib(4:6,i,new_path) CALL Small_To_L45 (pole_uvec, t_uvec) ELSE t_uvec = ai_pathlib(1:3,i,new_path) CALL Great_To_L45 (t_uvec) END IF END IF ! in path or not ELSE ! curve doesn't show IF (ai_in_path) THEN CALL End_L45_Path (close = .FALSE., stroke = .TRUE., fill = .FALSE.) END IF ! in path END IF ! curve shows or not !memory IF (ai_bent(i,new_path)) THEN last_uvec = ai_pathlib(4:6,i,new_path) ELSE last_uvec = ai_pathlib(1:3,i,new_path) END IF point0_shows = point1_shows END DO ! i = 1, new_segments IF (ai_in_path) CALL End_L45_Path( close = (ai_closed(path).AND.unchanged), & & stroke = .TRUE., fill = .FALSE.) END IF ! filled? or only stroked? ! Free up original L5 path (if still L5): IF (ai_path_level(path) == 5) THEN ai_path_level(path) = 0 ai_Ln_paths(5) = ai_Ln_paths(5) - 1 ai_total_paths = ai_total_paths - 1 END IF ! Free up temporary L(-1) path: ai_path_level(new_path) = 0 ai_total_paths = ai_total_paths - 1 CALL Process_L4_Paths END DO ! on all L5 paths presented END SUBROUTINE Process_L5_Paths SUBROUTINE Process_L5_Text (uvec, angle_radians, from_east, & & font_points, lr_fraction, ud_fraction, & & text) ! Check whether uvec is In_Rind, and if so, pass on to L4. IMPLICIT NONE REAL, DIMENSION(3), INTENT(IN) :: uvec REAL, INTENT(IN) :: angle_radians, & & lr_fraction, ud_fraction LOGICAL, INTENT(IN) :: from_east INTEGER, INTENT(IN) :: font_points CHARACTER*(*),INTENT(IN) :: text CHARACTER*1 :: result result = In_Rind (uvec) IF ((result == 'I').OR.(result == 'B')) THEN CALL Process_L4_Text (uvec, angle_radians, from_east, & & font_points, lr_fraction, ud_fraction, & & text) END IF ! no ELSE; text outside Rind is discarded. END SUBROUTINE Process_L5_Text SUBROUTINE Process_L6_Paths () ! Apply simple in-place (theta,phi) -> uvec transformation, ! and pass on to next level. IMPLICIT NONE INTEGER :: i, path REAL :: phi, phi2, theta, theta2 REAL, DIMENSION(3) :: uvec DO WHILE (ai_Ln_paths(6) > 0) path = Next_Path(level = 6) theta = ai_pathlib(1,0,path) phi = ai_pathlib(2,0,path) CALL ThetaPhi_2_Uvec (theta, phi, uvec) ai_pathlib(1:3,0,path) = uvec DO i = 1, ai_segments(path) theta = ai_pathlib(1,i,path) phi = ai_pathlib(2,i,path) IF (ai_bent(i,path)) THEN ! extract 3rd number before it gets overwritten: theta2 = ai_pathlib(3,i,path) phi2 = ai_pathlib(4,i,path) CALL ThetaPhi_2_Uvec (theta2,phi2, uvec) ai_pathlib(4:6,i,path) = uvec END IF CALL ThetaPhi_2_Uvec (theta,phi, uvec) ai_pathlib(1:3,i,path) = uvec END DO ! i = 1, segments ai_path_level(path) = 5 ai_Ln_paths(6) = ai_Ln_paths(6) - 1 ai_Ln_paths(5) = ai_Ln_paths(5) + 1 CALL Process_L5_Paths END DO ! while L6 paths remain END SUBROUTINE Process_L6_Paths SUBROUTINE Process_L6_Text (theta, phi, angle_radians, from_east, & & font_points, lr_fraction, ud_fraction, & & text) ! Transform position coordinates, and pass on. IMPLICIT NONE REAL, INTENT(IN) :: theta, phi, angle_radians, & & lr_fraction, ud_fraction LOGICAL, INTENT(IN) :: from_east INTEGER, INTENT(IN) :: font_points CHARACTER*(*),INTENT(IN) :: text REAL, DIMENSION(3) :: uvec CALL ThetaPhi_2_Uvec (theta, phi, uvec) CALL Process_L5_Text (uvec, angle_radians, from_east, & & font_points, lr_fraction, ud_fraction, & & text) END SUBROUTINE Process_L6_Text SUBROUTINE Process_L7_Paths () ! Apply simple in-place (lon,lat) -> (theta,phi) transformation, ! and pass on to next level. IMPLICIT NONE INTEGER :: i, path REAL :: lat, lon, phi, theta DO WHILE (ai_Ln_paths(7) > 0) path = Next_Path(level = 7) lon = ai_pathlib(1,0,path) lat = ai_pathlib(2,0,path) CALL LonLat_2_ThetaPhi (lon, lat, theta,phi) ai_pathlib(1,0,path) = theta ai_pathlib(2,0,path) = phi DO i = 1, ai_segments(path) lon = ai_pathlib(1,i,path) lat = ai_pathlib(2,i,path) CALL LonLat_2_ThetaPhi (lon,lat, theta,phi) ai_pathlib(1,i,path) = theta ai_pathlib(2,i,path) = phi IF (ai_bent(i,path)) THEN lon = ai_pathlib(3,i,path) lat = ai_pathlib(4,i,path) CALL LonLat_2_ThetaPhi (lon,lat, theta,phi) ai_pathlib(3,i,path) = theta ai_pathlib(4,i,path) = phi END IF END DO ! i = 1, segments ai_path_level(path) = 6 ai_Ln_paths(7) = ai_Ln_paths(7) - 1 ai_Ln_paths(6) = ai_Ln_paths(6) + 1 CALL Process_L6_Paths END DO ! while L7 paths remain END SUBROUTINE Process_L7_Paths SUBROUTINE Process_L7_Text (lon, lat, angle_radians, from_east, & & font_points, lr_fraction, ud_fraction, & & text) ! Transform position coordinates and pass on. IMPLICIT NONE REAL, INTENT(IN) :: lon, lat, angle_radians, & & lr_fraction, ud_fraction LOGICAL, INTENT(IN) :: from_east INTEGER, INTENT(IN) :: font_points CHARACTER*(*),INTENT(IN) :: text REAL :: theta, phi CALL LonLat_2_ThetaPhi (lon, lat, theta, phi) CALL Process_L6_Text (theta, phi, angle_radians, from_east, & & font_points, lr_fraction, ud_fraction, & & text) END SUBROUTINE Process_L7_Text SUBROUTINE Project (uvec, guide, x, y) ! Converts a position on the surface of a planet ! (uvec, a unit vector from the center of the sphere, ! in an orthogonal Cartesian system with the 1st axis ! at (0E,0N), 2nd axis at (90E,0N), and 3rd axis at 90N), ! into an (x,y) position (in meters) on the planet-sized ! projection plane. ! Uses parameters in global data of module Map_Projections. ! If "guide" is present, this is an additional surface ! position on the planet; the projection of "uvec" to ! (x,y) will take place AS IF uvec were on the same side of ! the cut as "guide", even if it really isn't! !========================================================== ! NOTE: It is NOT the job of Project to test inputs, ! give warning messages, create abends for projections ! to infinity or the complex plane, et cetera! ! Other programs (especially Process_L5_Paths) will ! guard against ever calling Project for troublesome ! points outside the walls of the "rind" of projectable ! surface points on the sphere. Therefore, Project ! can concentrate on simplicity and speed. !================================================================ IMPLICIT NONE REAL, DIMENSION(3), INTENT(IN) :: uvec REAL, DIMENSION(3), INTENT(IN), OPTIONAL :: guide REAL, INTENT(OUT) :: x, y REAL :: alpha, arc_radians, azimuth, beta, dot1, dot2, dot3, & & E, east_radians, equat, equat2, guide_east, & & latitude_radians, longitude_radians, & & radius_meters, rho, theta, & & up_radians, xpp, ypp dot1 = Dot (uvec, mp_pole_1_uvec) dot2 = Dot (uvec, mp_pole_2_uvec) dot3 = Dot (uvec, mp_pole_3_uvec) SELECT CASE (mp_projection_number) CASE (1) ! Mercator east_radians = ATAN2(dot2, dot1) equat = SQRT(dot1**2 + dot2**2) IF (PRESENT(guide)) THEN dot1 = Dot (guide, mp_pole_1_uvec) dot2 = Dot (guide, mp_pole_2_uvec) guide_east = ATAN2(dot2, dot1) IF (east_radians < (guide_east - Pi)) east_radians = east_radians + Two_Pi IF (east_radians > (guide_east + Pi)) east_radians = east_radians - Two_Pi END IF alpha = mp_radius_meters * east_radians up_radians = ATAN2(dot3, equat) beta = mp_radius_meters * LOG(TAN((up_radians * 0.5) + 0.78539816)) ! rotate from (alpha, beta) to "natural" (x,y): xpp = mp_upright(1,1)*alpha + mp_upright(1,2)*beta ypp = mp_upright(2,1)*alpha + mp_upright(2,2)*beta CASE (2) ! Lambert Conformal Conic ! Note: following are in a possibly-rotated system! equat = SQRT(dot1**2 + dot2**2) latitude_radians = ATAN2(dot3, equat) IF (equat > 0.) THEN longitude_radians = ATAN2(dot2, dot1) ELSE longitude_radians = mp_proj2_lambda0 END IF IF (PRESENT(guide)) THEN dot1 = Dot (guide, mp_pole_1_uvec) dot2 = Dot (guide, mp_pole_2_uvec) equat2 = dot1**2 + dot2**2 IF (equat2 > 0.) THEN guide_east = ATAN2(dot2, dot1) IF (longitude_radians < (guide_east - Pi)) longitude_radians = longitude_radians + Two_Pi IF (longitude_radians > (guide_east + Pi)) longitude_radians = longitude_radians - Two_Pi END IF END IF theta = mp_proj2_n * (longitude_radians - mp_proj2_lambda0) IF (ABS(latitude_radians - Pi_over_2) > 0.001) THEN rho = mp_radius_meters * mp_proj2_F / & & TAN((Pi_over_2 + latitude_radians)/2.)**mp_proj2_n ELSE ! avoid TAN(Pi_over_2) ---> Domain Error (crash) rho = 0. END IF alpha = rho * SIN(theta) beta = mp_proj2_rho0 - rho * COS(theta) ! rotate from (alpha, beta) to "natural" (x,y): xpp = mp_upright(1,1)*alpha + mp_upright(1,2)*beta ypp = mp_upright(2,1)*alpha + mp_upright(2,2)*beta CASE (3) ! Albers Equal-Area Conic ! Note: following are in a possibly-rotated system! equat = SQRT(dot1**2 + dot2**2) latitude_radians = ATAN2(dot3, equat) IF (equat > 0.) THEN longitude_radians = ATAN2(dot2, dot1) ELSE longitude_radians = mp_proj3_lambda0 END IF IF (PRESENT(guide)) THEN dot1 = Dot (guide, mp_pole_1_uvec) dot2 = Dot (guide, mp_pole_2_uvec) equat2 = dot1**2 + dot2**2 IF (equat2 > 0.) THEN guide_east = ATAN2(dot2, dot1) IF (longitude_radians < (guide_east - Pi)) longitude_radians = longitude_radians + Two_Pi IF (longitude_radians > (guide_east + Pi)) longitude_radians = longitude_radians - Two_Pi END IF END IF theta = mp_proj3_n * (longitude_radians - mp_proj3_lambda0) rho = mp_radius_meters * SQRT(MAX(0., mp_proj3_C - 2. * mp_proj3_n * SIN(latitude_radians))) / mp_proj3_n alpha = rho * SIN(theta) beta = mp_proj3_rho0 - rho * COS(theta) ! rotate from (alpha, beta) to "natural" (x,y): xpp = mp_upright(1,1)*alpha + mp_upright(1,2)*beta ypp = mp_upright(2,1)*alpha + mp_upright(2,2)*beta CASE (4) ! Polyconic ! Note: following are in a possibly-rotated system! equat = SQRT(dot1**2 + dot2**2) latitude_radians = ATAN2(dot3, equat) IF (equat > 0.) THEN longitude_radians = ATAN2(dot2, dot1) ELSE longitude_radians = mp_proj4_lon0_radians END IF IF (PRESENT(guide)) THEN dot1 = Dot (guide, mp_pole_1_uvec) dot2 = Dot (guide, mp_pole_2_uvec) equat2 = dot1**2 + dot2**2 IF (equat2 > 0.) THEN guide_east = ATAN2(dot2, dot1) IF (longitude_radians < (guide_east - Pi)) longitude_radians = longitude_radians + Two_Pi IF (longitude_radians > (guide_east + Pi)) longitude_radians = longitude_radians - Two_Pi END IF END IF E = (longitude_radians - mp_proj4_lon0_radians) * SIN(latitude_radians) IF (ABS(latitude_radians - Pi_over_2) > 0.0001) THEN alpha = mp_radius_meters * SIN(E) / TAN(latitude_radians) ELSE ! avoid TAN(Pi_over_2) ---> Domain error (crash)! alpha = 0. END IF beta = mp_radius_meters * (latitude_radians - mp_proj4_lat0_radians & & + (1. - COS(E)) / TAN(latitude_radians)) ! rotate from (alpha, beta) to "natural" (x,y): xpp = mp_upright(1,1)*alpha + mp_upright(1,2)*beta ypp = mp_upright(2,1)*alpha + mp_upright(2,2)*beta CASE (5) ! Geometric Conic ! Note: following are in a possibly-rotated system! equat = SQRT(dot1**2 + dot2**2) latitude_radians = ATAN2(dot3, equat) IF (equat > 0.) THEN longitude_radians = ATAN2(dot2, dot1) ELSE longitude_radians = 0. END IF IF (PRESENT(guide)) THEN dot1 = Dot (guide, mp_pole_1_uvec) dot2 = Dot (guide, mp_pole_2_uvec) equat2 = dot1**2 + dot2**2 IF (equat2 > 0.) THEN guide_east = ATAN2(dot2, dot1) IF (longitude_radians < (guide_east - Pi)) longitude_radians = longitude_radians + Two_Pi IF (longitude_radians > (guide_east + Pi)) longitude_radians = longitude_radians - Two_Pi END IF END IF theta = mp_proj5_n * (longitude_radians - mp_proj5_lambda0) IF (ABS((latitude_radians - mp_proj5_lat0_radians) - Pi_over_2) > 0.0001) THEN rho = mp_proj5_rtan - mp_radius_meters * TAN(latitude_radians - mp_proj5_lat0_radians) ELSE ! avoid TAN(Pi_over_2) ---> Domain error (crash)! rho = 0. END IF alpha = rho * SIN(theta) beta = mp_proj5_ypole - rho * COS(theta) ! rotate from (alpha, beta) to "natural" (x,y): xpp = mp_upright(1,1)*alpha + mp_upright(1,2)*beta ypp = mp_upright(2,1)*alpha + mp_upright(2,2)*beta CASE (6) ! Stereographic equat2 = dot2**2 + dot3**2 IF (equat2 > 0.) THEN azimuth = ATAN2 (dot2, dot3) arc_radians = Arc (mp_pole_1_uvec, uvec) radius_meters = 2. * mp_radius_meters * TAN(0.5 * arc_radians) xpp = radius_meters * SIN(azimuth) ypp = radius_meters * COS(azimuth) ELSE xpp = 0.0 ypp = 0.0 END IF CASE (7) ! Lambert Azimuthal Equal-Area equat2 = dot2**2 + dot3**2 IF (equat2 > 0.) THEN azimuth = ATAN2 (dot2, dot3) arc_radians = Arc (mp_pole_1_uvec, uvec) radius_meters = 2. * mp_radius_meters * SIN(0.5 * arc_radians) xpp = radius_meters * SIN(azimuth) ypp = radius_meters * COS(azimuth) ELSE xpp = 0.0 ypp = 0.0 END IF CASE (8) ! Azimuthal Equidistant equat2 = dot2**2 + dot3**2 IF (equat2 > 0.) THEN azimuth = ATAN2 (dot2, dot3) arc_radians = Arc (mp_pole_1_uvec, uvec) radius_meters = mp_radius_meters * arc_radians xpp = radius_meters * SIN(azimuth) ypp = radius_meters * COS(azimuth) ELSE xpp = 0.0 ypp = 0.0 END IF CASE (9) ! Orthographic xpp = mp_radius_meters * dot2 ypp = mp_radius_meters * dot3 CASE (10) ! Gnomonic equat2 = dot2**2 + dot3**2 IF (equat2 > 0.) THEN azimuth = ATAN2 (dot2, dot3) arc_radians = Arc (mp_pole_1_uvec, uvec) !WARNING: Paths must be edited before projection ! to avoid hitting, stepping over, or even approaching ! the dangerous singularity at arc_radians = Pi_over_2 !!! radius_meters = mp_radius_meters * TAN(arc_radians) xpp = radius_meters * SIN(azimuth) ypp = radius_meters * COS(azimuth) ELSE xpp = 0.0 ypp = 0.0 END IF CASE DEFAULT WRITE (*,"(' ERROR: Map projection [',A,'] is not supported.')") mp_projection CALL Traceback END SELECT IF (mp_user_xy) THEN ! correct from "natural", default (x,y) to user's (x,y): x = mp_pp_2_xy_matrix(1,1)*xpp + mp_pp_2_xy_matrix(1,2)*ypp + mp_pp_2_xy_vector(1) y = mp_pp_2_xy_matrix(2,1)*xpp + mp_pp_2_xy_matrix(2,2)*ypp + mp_pp_2_xy_vector(2) ELSE x = xpp y = ypp END IF END SUBROUTINE Project SUBROUTINE Reject (x_meters, y_meters, success, uvec) ! Opposite of "SUBROUTINE Project": ! Converts a position on the surface of a planet-sized ! projection plane, with (x_meters,y_meters) coordinates ! to "uvec"(a unit vector from the center of the sphere, ! in an orthogonal Cartesian system with the 1st axis ! at (0E,0N), 2nd axis at (90E,0N), and 3rd axis at 90N). ! Uses parameters in global data of module Map_Projections. ! If the (x,y) point requested has no inverse projection, ! then success = .FALSE. on return. IMPLICIT NONE REAL, INTENT(IN) :: x_meters, y_meters LOGICAL, INTENT(OUT) :: success REAL, DIMENSION(3), INTENT(OUT) :: uvec INTEGER :: i LOGICAL :: converged, converging REAL :: A, alpha, arc_radians, argument, azimuth, & & B, beta, beta_max, beta_min, & & change, dot1, dot2, dot3, east_radians, last_change, & & latitude_radians, longitude_radians, new_lat_radians, & & radius_meters, rho, tanlat, theta, up_radians, xpp, ypp REAL, DIMENSION(3) :: omega_uvec IF (mp_user_xy) THEN ! correct from user's (x,y) to "natural", default (x,y): xpp = mp_xy_2_pp_matrix(1,1)*x_meters + mp_xy_2_pp_matrix(1,2)*y_meters + mp_xy_2_pp_vector(1) ypp = mp_xy_2_pp_matrix(2,1)*x_meters + mp_xy_2_pp_matrix(2,2)*y_meters + mp_xy_2_pp_vector(2) ELSE xpp = x_meters ypp = y_meters END IF SELECT CASE (mp_projection_number) CASE (1) ! Mercator ! rotate from "natural" (x,y) to (alpha, beta): alpha = mp_downright(1,1)*xpp + mp_downright(1,2)*ypp beta = mp_downright(2,1)*xpp + mp_downright(2,2)*ypp ! undo the Mercator projection: east_radians = alpha / mp_radius_meters IF (ABS(east_radians) > Pi) THEN success = .FALSE. RETURN END IF up_radians = Pi_over_2 - 2. * ATAN(EXP(-beta/mp_radius_meters)) IF (ABS(up_radians) > 1.221381) THEN ! 69.98 degrees success = .FALSE. ELSE success = .TRUE. dot1 = COS(up_radians) * COS(east_radians) dot2 = COS(up_radians) * SIN(east_radians) dot3 = SIN(up_radians) uvec(1) = dot1*mp_pole_1_uvec(1) + dot2*mp_pole_2_uvec(1) + dot3*mp_pole_3_uvec(1) uvec(2) = dot1*mp_pole_1_uvec(2) + dot2*mp_pole_2_uvec(2) + dot3*mp_pole_3_uvec(2) uvec(3) = dot1*mp_pole_1_uvec(3) + dot2*mp_pole_2_uvec(3) + dot3*mp_pole_3_uvec(3) END IF CASE (2) ! Lambert Conformal Conic ! rotate from "natural" (x,y) to (alpha, beta) [cone pole up]: alpha = mp_downright(1,1)*xpp + mp_downright(1,2)*ypp beta = mp_downright(2,1)*xpp + mp_downright(2,2)*ypp ! UNDO the Lambert Conformal Conic projection: rho = SIGN(1.,mp_proj2_n) * SQRT(alpha**2 + (mp_proj2_rho0 - beta)**2) IF (mp_proj2_n > 0.) THEN theta = ATAN2(alpha, (mp_proj2_rho0 - beta)) ELSE theta = ATAN2(-alpha, (-mp_proj2_rho0 + beta)) END IF longitude_radians = (theta / mp_proj2_n) + mp_proj2_lambda0 IF (rho == 0.) THEN latitude_radians = SIGN(Pi_over_2, mp_proj2_n) ELSE latitude_radians = 2. * ATAN((mp_radius_meters * mp_proj2_f / rho)**(1./mp_proj2_n)) - Pi_over_2 END IF ! Note: both of the the above may be in a rotated system: uvec = mp_pole_1_uvec*COS(latitude_radians)*COS(longitude_radians) + & & mp_pole_2_uvec*COS(latitude_radians)*SIN(longitude_radians) + & & mp_pole_3_uvec*SIN(latitude_radians) success = .TRUE. ! provisionally, unless changed below: IF (mp_using_top) THEN IF (Dot(uvec,mp_top_uvec) > mp_top_maxdot) success = .FALSE. END IF IF (mp_using_bottom) THEN IF (Dot(uvec,mp_bottom_uvec) > mp_bottom_maxdot) success = .FALSE. END IF CASE (3) ! Albers Equal-Area Conic ! rotate from "natural" (x,y) to (alpha, beta) [cone pole up]: alpha = mp_downright(1,1)*xpp + mp_downright(1,2)*ypp beta = mp_downright(2,1)*xpp + mp_downright(2,2)*ypp ! UNDO the Albers Equal-Area Conic projection: rho = SQRT(alpha**2 + (mp_proj3_rho0 - beta)**2) IF (mp_proj3_n > 0.) THEN theta = ATAN2(alpha, (mp_proj3_rho0 - beta)) ELSE theta = ATAN2(-alpha, (-mp_proj3_rho0 + beta)) END IF longitude_radians = (theta / mp_proj3_n) + mp_proj3_lambda0 argument = (mp_proj3_C - (rho * mp_proj3_n / mp_radius_meters)**2) / (2. * mp_proj3_n) IF (ABS(argument) > 1.) THEN success = .FALSE. ELSE ! usual case latitude_radians = ASIN(argument) ! Note: both of the the above may be in a rotated system: uvec = mp_pole_1_uvec*COS(latitude_radians)*COS(longitude_radians) + & & mp_pole_2_uvec*COS(latitude_radians)*SIN(longitude_radians) + & & mp_pole_3_uvec*SIN(latitude_radians) success = .TRUE. ! unless changed below: IF (mp_using_top) THEN IF (Dot(uvec,mp_top_uvec) > mp_top_maxdot) success = .FALSE. END IF IF (mp_using_bottom) THEN IF (Dot(uvec,mp_bottom_uvec) > mp_bottom_maxdot) success = .FALSE. END IF END IF ! ABS(argument) <= 1. CASE (4) ! Polyconic ! rotate from "natural" (x,y) to (alpha, beta) [cone pole up]: alpha = mp_downright(1,1)*xpp + mp_downright(1,2)*ypp beta = mp_downright(2,1)*xpp + mp_downright(2,2)*ypp ! UNDO the Polyconic projection, IF POSSIBLE. beta_min = -mp_radius_meters * mp_proj4_lat0_radians beta_max = mp_radius_meters * (Pi_over_2 - mp_proj4_lat0_radians) IF (beta < beta_min) THEN ! wrong hemisphere success = .FALSE. ELSE IF (beta == beta_min) THEN ! equator success = .TRUE. latitude_radians = 0. longitude_radians = (alpha/mp_radius_meters) + mp_proj4_lon0_radians ELSE IF (beta > beta_max) THEN ! above pole success = .FALSE. ELSE IF (beta == beta_max) THEN ! at pole latitude_radians = Pi_over_2 longitude_radians = mp_proj4_lon0_radians ELSE ! general case; ! Only an iterated formula is available. It will fail for ! "latitudes" (with respect to the cone pole) of 0 or less, ! and also for "longitudes" (with respect to the projection point) ! 90 or more degrees away from the prime meridian. A = mp_proj4_lat0_radians + beta / mp_radius_meters B = (alpha / mp_radius_meters)**2 + A**2 latitude_radians = A converging = .FALSE. last_change = 1.0 iterate: DO i = 1, 10 tanlat = TAN(latitude_radians) new_lat_radians = latitude_radians - (A*(latitude_radians*tanlat + 1.) & & - latitude_radians - 0.5*(latitude_radians**2 + B)*tanlat) & & / (((latitude_radians - A)/tanlat) - 1.) change = ABS(new_lat_radians - latitude_radians) IF (change > 1.) EXIT iterate ! probably diverging converging = (change < last_change) IF (converging.AND.(change < 1.E-5)) EXIT iterate latitude_radians = new_lat_radians ! memory for next iteration last_change = change END DO iterate converged = converging .AND. (change < 1.E-4) IF (converged) THEN IF (latitude_radians < 0.) THEN ! wrong hemisphere success = .FALSE. ELSE IF (latitude_radians == 0.) THEN success = .TRUE. longitude_radians = (alpha/mp_radius_meters) + mp_proj4_lon0_radians ELSE IF (latitude_radians > Pi_over_2) THEN ! above pole??? success = .FALSE. ELSE IF (latitude_radians == Pi_over_2) THEN success = .TRUE. longitude_radians = mp_proj4_lon0_radians ELSE ! general case success = .TRUE. longitude_radians = ((ASIN(alpha*tanlat/mp_radius_meters)) & & / SIN(latitude_radians)) + mp_proj4_lon0_radians END IF ! at pole, or general case ELSE ! not converging success = .FALSE. END IF ! converging or not END IF ! choice of solution method for lon, lat in rotated system IF (success) THEN ! check that point is in bounds uvec = mp_pole_1_uvec*COS(latitude_radians)*COS(longitude_radians) + & & mp_pole_2_uvec*COS(latitude_radians)*SIN(longitude_radians) + & & mp_pole_3_uvec*SIN(latitude_radians) IF (mp_using_left) THEN IF (Dot(uvec,mp_left_uvec) > mp_left_maxdot) success = .FALSE. END IF IF (mp_using_right) THEN IF (Dot(uvec,mp_right_uvec) > mp_right_maxdot) success = .FALSE. END IF IF (mp_using_bottom) THEN IF (Dot(uvec,mp_bottom_uvec) > mp_bottom_maxdot) success = .FALSE. END IF END IF ! apparent success CASE (5) ! Geometric Conic ! rotate from "natural" (x,y) to (alpha, beta) [cone pole up]: alpha = mp_downright(1,1)*xpp + mp_downright(1,2)*ypp beta = mp_downright(2,1)*xpp + mp_downright(2,2)*ypp ! undo the geometric conic projection rho = SQRT(alpha**2 + (mp_proj5_ypole - beta)**2) theta = ATAN2(alpha, (mp_proj5_ypole - beta)) longitude_radians = (theta / mp_proj5_n) + mp_proj5_lambda0 IF (rho == 0.) THEN latitude_radians = SIGN(Pi_over_2, mp_proj5_n) ELSE latitude_radians = mp_proj5_lat0_radians + ATAN2(mp_proj5_rtan - rho, mp_radius_meters) END IF ! Note: both of the the above may be in a rotated system: uvec = mp_pole_1_uvec*COS(latitude_radians)*COS(longitude_radians) + & & mp_pole_2_uvec*COS(latitude_radians)*SIN(longitude_radians) + & & mp_pole_3_uvec*SIN(latitude_radians) success = .TRUE. IF (mp_using_top) THEN IF (Dot(uvec,mp_top_uvec) > mp_top_maxdot) success = .FALSE. END IF IF (mp_using_bottom) THEN IF (Dot(uvec,mp_bottom_uvec) > mp_bottom_maxdot) success = .FALSE. END IF CASE (6) ! Stereographic radius_meters = SQRT (xpp**2 + ypp**2) IF (radius_meters > 0.) THEN azimuth = ATAN2(xpp, ypp) arc_radians = 2.0 * ATAN(0.5 * radius_meters / mp_radius_meters) CALL Turn_To (azimuth, mp_pole_1_uvec, arc_radians, & ! inputs & omega_uvec, uvec) ELSE azimuth = 0. arc_radians = 0. uvec = mp_pole_1_uvec END IF IF (arc_radians <= Pi_over_2) THEN ! NOTE: This projection CAN show the whole planet, ! but the far hemisphere is terribly distorted. ! Thus, we map only the near hemisphere. success = .TRUE. ELSE success = .FALSE. END IF CASE (7) ! Lambert Azimuthal Equal-Area radius_meters = SQRT (xpp**2 + ypp**2) IF (radius_meters > 0.) THEN azimuth = ATAN2(xpp, ypp) arc_radians = 2.0 * ASIN(0.5 * radius_meters / mp_radius_meters) CALL Turn_To (azimuth, mp_pole_1_uvec, arc_radians, & ! inputs & omega_uvec, uvec) ELSE azimuth = 0. arc_radians = 0. uvec = mp_pole_1_uvec END IF IF (arc_radians <= Pi_over_2) THEN success = .TRUE. ELSE success = .FALSE. END IF CASE (8) ! Azimuthal Equidistant radius_meters = SQRT (xpp**2 + ypp**2) IF (radius_meters > 0.) THEN azimuth = ATAN2(xpp, ypp) arc_radians = radius_meters / mp_radius_meters CALL Turn_To (azimuth, mp_pole_1_uvec, arc_radians, & ! inputs & omega_uvec, uvec) ELSE azimuth = 0. arc_radians = 0. uvec = mp_pole_1_uvec END IF IF (arc_radians <= Pi_over_2) THEN ! NOTE: This projection CAN show the whole planet, ! but the far hemisphere is terribly distorted. ! Thus, we map only the near hemisphere. success = .TRUE. ELSE success = .FALSE. END IF CASE (9) ! Orthographic radius_meters = SQRT (xpp**2 + ypp**2) IF (radius_meters <= mp_radius_meters) THEN success = .TRUE. dot2 = xpp / mp_radius_meters dot3 = ypp / mp_radius_meters dot1 = SQRT (1.00 - dot2**2 - dot3**2) uvec = dot1*mp_pole_1_uvec + dot2*mp_pole_2_uvec + dot3*mp_pole_3_uvec ELSE success = .FALSE. END IF CASE (10) ! Gnomonic radius_meters = SQRT (xpp**2 + ypp**2) IF (radius_meters > 0.) THEN azimuth = ATAN2(xpp, ypp) arc_radians = ATAN(radius_meters / mp_radius_meters) CALL Turn_To (azimuth, mp_pole_1_uvec, arc_radians, & ! inputs & omega_uvec, uvec) ELSE uvec = mp_pole_1_uvec END IF success = .TRUE. CASE DEFAULT WRITE (*,"(' ERROR: Map projection [',A,'] is not supported.')") TRIM(mp_projection) CALL Traceback END SELECT END SUBROUTINE Reject 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 Restore_mp_State () IMPLICIT NONE IF (mp_state_saved) THEN mp_zoom_defined = sp_zoom_defined mp_scale_denominator = sp_scale_denominator mp_x_center_meters = sp_x_center_meters mp_y_center_meters = sp_y_center_meters mp_xy_wrt_page_radians = sp_xy_wrt_page_radians mp_meters_per_point = sp_meters_per_point mp_m_2_pts_matrix = sp_m_2_pts_matrix mp_pts_2_m_matrix = sp_pts_2_m_matrix mp_m_2_pts_vector = sp_m_2_pts_vector mp_pts_2_m_vector = sp_pts_2_m_vector mp_last_L3_x_meters = sp_last_L3_x_meters mp_last_L3_y_meters = sp_last_L3_y_meters mp_last_uvec = sp_last_uvec mp_last_theta = sp_last_theta mp_last_phi = sp_last_phi mp_last_lon = sp_last_lon mp_last_lat = sp_last_lat mp_using_top = sp_using_top mp_using_left = sp_using_left mp_using_right = sp_using_right mp_using_bottom = sp_using_bottom mp_bottom_uvec = sp_bottom_uvec mp_right_uvec = sp_right_uvec mp_top_uvec = sp_top_uvec mp_left_uvec = sp_left_uvec mp_bottom_maxdot = sp_bottom_maxdot mp_right_maxdot = sp_right_maxdot mp_top_maxdot = sp_top_maxdot mp_left_maxdot = sp_left_maxdot mp_bottom_applies_uvec = sp_bottom_applies_uvec mp_right_applies_uvec = sp_right_applies_uvec mp_top_applies_uvec = sp_top_applies_uvec mp_left_applies_uvec = sp_left_applies_uvec mp_first_top_uvec = sp_first_top_uvec mp_last_top_uvec = sp_last_top_uvec mp_first_left_uvec = sp_first_left_uvec mp_last_left_uvec = sp_last_left_uvec mp_first_right_uvec = sp_first_right_uvec mp_last_right_uvec = sp_last_right_uvec mp_first_bottom_uvec = sp_first_bottom_uvec mp_last_bottom_uvec = sp_last_bottom_uvec mp_next_wall = sp_next_wall mp_projection_defined = sp_projection_defined mp_projection = sp_projection mp_projection_number = sp_projection_number mp_radius_meters = sp_radius_meters mp_pole_1_uvec = sp_pole_1_uvec mp_pole_2_uvec = sp_pole_2_uvec mp_pole_3_uvec = sp_pole_3_uvec mp_projpoint_uvec = sp_projpoint_uvec mp_proj2_F = sp_proj2_F mp_proj2_lambda0 = sp_proj2_lambda0 mp_proj2_n = sp_proj2_n mp_proj2_rho0 = sp_proj2_rho0 mp_proj3_C = sp_proj3_C mp_proj3_lambda0 = sp_proj3_lambda0 mp_proj3_n = sp_proj3_n mp_proj3_rho0 = sp_proj3_rho0 mp_proj4_lat0_radians = sp_proj4_lat0_radians mp_proj4_lon0_radians = sp_proj4_lon0_radians mp_proj5_lambda0 = sp_proj5_lambda0 mp_proj5_lat0_radians = sp_proj5_lat0_radians mp_proj5_n = sp_proj5_n mp_proj5_rtan = sp_proj5_rtan mp_proj5_ypole = sp_proj5_ypole mp_user_xy = sp_user_xy mp_x_projpoint_meters = sp_x_projpoint_meters mp_y_projpoint_meters = sp_y_projpoint_meters mp_y_azimuth_radians = sp_y_azimuth_radians mp_pp_2_xy_matrix = sp_pp_2_xy_matrix mp_xy_2_pp_matrix = sp_xy_2_pp_matrix mp_pp_2_xy_vector = sp_pp_2_xy_vector mp_xy_2_pp_vector = sp_xy_2_pp_vector mp_upright = sp_upright mp_downright = sp_downright ELSE ! no state available to restore WRITE (*,"(' ERROR: Must CALL Save_mp_State() before CALL Restore_mp_State()')") CALL Traceback END IF ! state was previously saved, or not END SUBROUTINE Restore_mp_State SUBROUTINE Save_mp_State () IMPLICIT NONE IF (mp_zoom_defined) THEN mp_state_saved = .TRUE. !(statement above is placed here for greater prominence) sp_zoom_defined = mp_zoom_defined sp_scale_denominator = mp_scale_denominator sp_x_center_meters = mp_x_center_meters sp_y_center_meters = mp_y_center_meters sp_xy_wrt_page_radians = mp_xy_wrt_page_radians sp_meters_per_point = mp_meters_per_point sp_m_2_pts_matrix = mp_m_2_pts_matrix sp_pts_2_m_matrix = mp_pts_2_m_matrix sp_m_2_pts_vector = mp_m_2_pts_vector sp_pts_2_m_vector = mp_pts_2_m_vector sp_last_L3_x_meters = mp_last_L3_x_meters sp_last_L3_y_meters = mp_last_L3_y_meters sp_last_uvec = mp_last_uvec sp_last_theta = mp_last_theta sp_last_phi = mp_last_phi sp_last_lon = mp_last_lon sp_last_lat = mp_last_lat sp_using_top = mp_using_top sp_using_left = mp_using_left sp_using_right = mp_using_right sp_using_bottom = mp_using_bottom sp_bottom_uvec = mp_bottom_uvec sp_right_uvec = mp_right_uvec sp_top_uvec = mp_top_uvec sp_left_uvec = mp_left_uvec sp_bottom_maxdot = mp_bottom_maxdot sp_right_maxdot = mp_right_maxdot sp_top_maxdot = mp_top_maxdot sp_left_maxdot = mp_left_maxdot sp_bottom_applies_uvec = mp_bottom_applies_uvec sp_right_applies_uvec = mp_right_applies_uvec sp_top_applies_uvec = mp_top_applies_uvec sp_left_applies_uvec = mp_left_applies_uvec sp_first_top_uvec = mp_first_top_uvec sp_last_top_uvec = mp_last_top_uvec sp_first_left_uvec = mp_first_left_uvec sp_last_left_uvec = mp_last_left_uvec sp_first_right_uvec = mp_first_right_uvec sp_last_right_uvec = mp_last_right_uvec sp_first_bottom_uvec = mp_first_bottom_uvec sp_last_bottom_uvec = mp_last_bottom_uvec sp_next_wall = mp_next_wall sp_projection_defined = mp_projection_defined sp_projection = mp_projection sp_projection_number = mp_projection_number sp_radius_meters = mp_radius_meters sp_pole_1_uvec = mp_pole_1_uvec sp_pole_2_uvec = mp_pole_2_uvec sp_pole_3_uvec = mp_pole_3_uvec sp_projpoint_uvec = mp_projpoint_uvec sp_proj2_F = mp_proj2_F sp_proj2_lambda0 = mp_proj2_lambda0 sp_proj2_n = mp_proj2_n sp_proj2_rho0 = mp_proj2_rho0 sp_proj3_C = mp_proj3_C sp_proj3_lambda0 = mp_proj3_lambda0 sp_proj3_n = mp_proj3_n sp_proj3_rho0 = mp_proj3_rho0 sp_proj4_lat0_radians = mp_proj4_lat0_radians sp_proj4_lon0_radians = mp_proj4_lon0_radians sp_proj5_lambda0 = mp_proj5_lambda0 sp_proj5_lat0_radians = mp_proj5_lat0_radians sp_proj5_n = mp_proj5_n sp_proj5_rtan = mp_proj5_rtan sp_proj5_ypole = mp_proj5_ypole sp_user_xy = mp_user_xy sp_x_projpoint_meters = mp_x_projpoint_meters sp_y_projpoint_meters = mp_y_projpoint_meters sp_y_azimuth_radians = mp_y_azimuth_radians sp_pp_2_xy_matrix = mp_pp_2_xy_matrix sp_xy_2_pp_matrix = mp_xy_2_pp_matrix sp_pp_2_xy_vector = mp_pp_2_xy_vector sp_xy_2_pp_vector = mp_xy_2_pp_vector sp_upright = mp_upright sp_downright = mp_downright ELSE ! no state defined yet WRITE (*,"(' ERROR: Must Set_Zoom before Save_mp_State')") CALL Traceback END IF ! a state of mp is defined, or not END SUBROUTINE Save_mp_State 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 Trace_Rind (from_t_wall,clockwise, & & to_point_uvec,to_t_wall) ! Continues an open path on level 4 by adding boundary segments ! along the edge of the rind (defined by Set_[projection] and Corners) ! until reaching "to_point_uvec", in direction "clockwise"(?). ! Parameters "from_t_wall" and "to_t_wall" give the intial and ! final points in a dimensionless coordinate system where the ! whole-number part of "t_wall" is: ! 1 = bottom; 2 = right; 3 = top; 4 = left wall of rind, and ! fractional part gives relative position along this wall ! in the counterclockwise direction (e.g., 1.01 is near left ! end of bottom wall; 3.99 is near left end of top wall). IMPLICIT NONE LOGICAL, INTENT(IN) :: clockwise REAL, INTENT(IN) :: from_t_wall, to_t_wall REAL, DIMENSION(3), INTENT(IN) :: to_point_uvec INTEGER :: current_wall, end_wall, i, start_wall LOGICAL first_pass, one_step REAL :: end_f, start_f REAL, DIMENSION(3) :: first_uvec, last_uvec, pole_uvec ! check for "just touching" boundary: IF (from_t_wall == to_t_wall) RETURN start_wall = MAX(1,MIN(4,INT(from_t_wall))) end_wall = MAX(1,MIN(4,INT (to_t_wall))) start_f = from_t_wall - start_wall end_f = to_t_wall - end_wall !Decide whether this can be done with a single small-circle: one_step = (start_wall == end_wall) .AND. & & ((mp_next_wall(0,start_wall)*mp_next_wall(1,start_wall) == 0) & ! no escape & .OR. (clockwise.AND.(end_f < start_f)) & & .OR. ((.NOT.clockwise).AND.(end_f > start_f))) IF (one_step) THEN SELECT CASE(start_wall) CASE(1) pole_uvec = mp_bottom_uvec CASE(2) pole_uvec = mp_right_uvec CASE(3) pole_uvec = mp_top_uvec CASE(4) pole_uvec = mp_left_uvec END SELECT IF (clockwise) THEN ! clockwise around the boundary is ccw around the pole CALL Small_To_L45 (pole_uvec, to_point_uvec) ELSE ! ccw around the boundary is ccw around the anti-pole pole_uvec = -pole_uvec CALL Small_To_L45 (pole_uvec, to_point_uvec) END IF ELSE ! complex circuit needed: current_wall = start_wall first_pass = .TRUE. DO i = 1, 5 ! enough to cover all cases SELECT CASE(current_wall) CASE(1) pole_uvec = mp_bottom_uvec first_uvec = mp_first_bottom_uvec last_uvec = mp_last_bottom_uvec CASE(2) pole_uvec = mp_right_uvec first_uvec = mp_first_right_uvec last_uvec = mp_last_right_uvec CASE(3) pole_uvec = mp_top_uvec first_uvec = mp_first_top_uvec last_uvec = mp_last_top_uvec CASE(4) pole_uvec = mp_left_uvec first_uvec = mp_first_left_uvec last_uvec = mp_last_left_uvec END SELECT IF (.NOT.clockwise) THEN ! swap pole and last point pole_uvec = -pole_uvec last_uvec = first_uvec END IF IF ((.NOT.first_pass).AND.(current_wall == end_wall)) THEN ! go to end-point and exit CALL Small_To_L45 (pole_uvec, to_point_uvec) RETURN ELSE ! go to corner point CALL Small_To_L45 (pole_uvec, last_uvec) END IF !Prepare for next loop: first_pass = .FALSE. IF (clockwise) THEN current_wall = mp_next_wall(0,current_wall) ELSE current_wall = mp_next_wall(1,current_wall) END IF END DO END IF END SUBROUTINE Trace_Rind 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 REAL FUNCTION Vector_Azimuth (site_uvec, vector) !Returns azimuth (in radians, clockwise from North) !of a vector (which does not need to be horizontal, !but which should not be vertical! IMPLICIT NONE REAL, DIMENSION(3), INTENT(IN) :: site_uvec, vector REAL, DIMENSION(3) :: theta_uvec, phi_uvec REAL :: N_part, E_part CALL Local_Theta(site_uvec, theta_uvec) CALL Local_Phi (site_uvec, phi_uvec) N_part = -Dot(theta_uvec, vector) E_part = Dot(phi_uvec, vector) IF ((N_part /= 0.0).OR.(E_part /= 0.0)) THEN Vector_Azimuth = ATAN2(E_part, N_part) ELSE Vector_Azimuth = 0.0 END IF END FUNCTION Vector_Azimuth END MODULE Map_Projections !===============================================