MODULE DMap_Projections ! ! Initial letter "D" indicates a DOUBLE PRECISION (REAL*8) version, ! created 2015.02 as part of a systematic upgrade of several of my codes ! (Shells, NeoKinema, FiniteMap, & NeoKineMap) to 64-bit precision. ! It is intended that the "D" FUNCTIONs and SUBROUTINEs here should be ! logically equivalent to those in MODULE Map_Projections, ! just more precise. Naturally, they now all take REAL*8 arguments ! in place of the old REAL arguments. ! ! Updated version copyright 2015 by Peter Bird and the ! Regents of the University of California. ! !---------------------------------------------------------------------- ! ! 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 DAdobe_Illustrator to produce graphics files readable by ! Adobe Illustrator for Windows, versions 4, 7, ...? USE DAdobe_Illustrator ! provided as file DAdobe_Illustrator.f90 ! ! Original code of MODULE Map_Projections ! by Peter Bird, UCLA, May 1997 - September 1999. ! Upgraded to "D" precision in 2015. ! Copyright (c) 1997, 1998, 1999, 2015 by ! Peter Bird and the Regents of the University of California. !----------------------------------------------------------------- ! ! USAGE NOTES FOR THIS MODULE: ! ------------------------------ ! 1. Always begin by initializing module DAdobe_Illustrator, with ! CALL DSelect_Paper ! CALL DSet_Background ! CALL DDefine_Margins ! CALL DBegin_Page ! CALL DSet_Window [ optional; DBegin_Page computes best window ] ! ! 2. Before plotting any map data, you MUST: ! CALL DSet_Zoom (scale_denominator, x_center_meters, & ! & y_center_meters, xy_wrt_page_radians) ! [ OR ] ! CALL DSet_Zoom (scale_denominator) ! This call must occur after CALL DSet_Window; if you re-call ! DSet_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.D6 to 1.D9 ?) ! 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 DAdobe_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 DNew_L3_Path [ begins a path; must be first ] ! CALL DLine_To_L3 [ optional; may be repeated ] ! CALL DCurve_To_L3 [ optional; may be repeated ] ! CALL DEnd_L3_Path [ ends the path; must be last ] ! possibly in combination with: ! CALL DCircle_on_L3 [ either before or after drawing path(s) ] ! possibly in combination with: ! CALL DL3_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 DSet_Mercator (...) ! [ OR ] ! CALL DSet_Conic (...) ! CALL DSet_Stereographic (...) ! CALL DSet_Orthographic (...) ! ! 6. Once the map projection is defined (above), you may draw ! on the spherical planet surface with a combination of: ! CALL DNew_L567_Path [ begins a path; must be first ] ! CALL DGreat_To_L567 [ optional; may be repeated ] ! CALL DSmall_To_L567 [ optional; may be repeated ] ! CALL DEnd_L567_Path [ ends the path; must be last ] ! possibly in combination with: ! CALL DL567_Text [ either before or after drawing path(s) ] ! Note that Great_To_L567 adds a great-circle arc to the path, ! while DSmall_To_L567 adds a small-circle arc to the path. ! ! 7. Always end by closing the DAdobe_Illustrator output file: ! CALL DEnd_Page ! ! CONTENTS OF THIS MODULE: ! -------------------------- ! ! USER SUBROUTINES (in a possible order-of-use): ! SUBROUTINE DSet_Zoom [ Always required, first! ] ! SUBROUTINE DNew_L3_Path ! SUBROUTINE DLine_To_L3 ! SUBROUTINE DCurve_To_L3 ! SUBROUTINE DEnd_L3_Path [ and/or ] ! SUBROUTINE DCircle_on_L3 [ and/or ] ! SUBROUTINE DL3_Text ! ---- transition from projection plane (level 3) to ! surface of spherical planet (levels 5,6,7) ------ ! SUBROUTINE DSet_Mercator ! [ OR ] SUBROUTINE DSet_Lambert_Conformal_Conic ! [ OR ] SUBROUTINE DSet_Albers_Equal_Area_Conic ! [ OR ] SUBROUTINE DSet_Polyconic ! [ OR ] SUBROUTINE DSet_Geometric_Conic ! [ OR ] SUBROUTINE DSet_Stereographic ! [ OR ] SUBROUTINE DSet_Lambert_Azimuthal_EqualArea ! [ OR ] SUBROUTINE DSet_Azimuthal_Equidistant ! [ OR ] SUBROUTINE DSet_Orthographic ! [ OR ] SUBROUTINE DSet_Gnomonic ! [ use ONE of the above ] ! SUBROUTINE DNew_L45_Path [(Level 4 for internal use) ] ! SUBROUTINE DGreat_To_L45 [ Level 5 accepts vectors ] ! SUBROUTINE DSafe_Small_To_L45 [ See also DSmall_to_L45. ] ! SUBROUTINE DSmall_Through_L45 ! SUBROUTINE DSmall_To_L45 [ See also DSafe_Small_to_L45. ] ! SUBROUTINE DEnd_L45_Path ! SUBROUTINE DL5_Text [(There is no L4_Text.) ] ! SUBROUTINE DNew_L67_Path [ Level 6 accepts theta/phi;] ! SUBROUTINE DGreat_To_L67 [ Level 7 accepts lon/lat. ] ! SUBROUTINE DSmall_To_L67 ! SUBROUTINE DEnd_L67_Path ! SUBROUTINE DL67_Text ! ! UTILITY SUBROUTINES AND FUNCTIONS: ! REAL*8 FUNCTION DArc ! SUBROUTINE DCircles_Intersect ! LOGICAL FUNCTION DClockways ! REAL*8 FUNCTION DCompass ! REAL*8 FUNCTION DConformal_Deflation ! SUBROUTINE DCross ! REAL*8 FUNCTION DDot ! REAL*8 FUNCTION DEasting ! CHARACTER*1 FUNCTION DIn_Rind ! LOGICAL FUNCTION DL5_In_Window ! REAL*8 FUNCTION DLength ! SUBROUTINE DLocal_Phi ! SUBROUTINE DLocal_Theta ! SUBROUTINE DLonLat_2_ThetaPhi ! SUBROUTINE DLonLat_2_Uvec ! SUBROUTINE DMagnitude ! SUBROUTINE DMake_Uvec ! SUBROUTINE DMeters_2_Points ! SUBROUTINE DNorthEast_Convention ! SUBROUTINE DOld_Complex_Process_L5_Paths ! SUBROUTINE DPoints_2_LonLat ! SUBROUTINE DPoints_2_Meters ! SUBROUTINE DProcess_L3_Paths ! SUBROUTINE DProcess_L3_Text ! SUBROUTINE DProcess_L4_Paths ! SUBROUTINE DProcess_L4_Text ! SUBROUTINE DProcess_L5_Paths ! SUBROUTINE DProcess_L5_Text ! SUBROUTINE DProcess_L6_Paths ! SUBROUTINE DProcess_L6_Text ! SUBROUTINE DProcess_L7_Paths ! SUBROUTINE DProcess_L7_Text ! SUBROUTINE DProject ! SUBROUTINE DReject ! REAL*8 FUNCTION DRelative_Compass ! SUBROUTINE DRestore_mp_State ! SUBROUTINE DSave_mp_State ! SUBROUTINE DThetaPhi_2_LonLat ! SUBROUTINE DThetaPhi_2_Uvec ! SUBROUTINE DTrace_Rind ! SUBROUTINE DTurn_To ! SUBROUTINE DUvec_2_LonLat ! SUBROUTINE DUvec_2_ThetaPhi ! REAL*8 FUNCTION DVector_Azimuth !----------------------------------------------------------------- !declares IMPLICIT NONE DOUBLE PRECISION, PARAMETER :: Two_Pi = 6.28318530717959D0 DOUBLE PRECISION, PARAMETER :: Pi = 3.14159265358979D0 ! These values were confirmed DOUBLE PRECISION, PARAMETER :: Pi_over_2 = 1.57079632679490D0 ! with PROGRAM Check_Pi. DOUBLE PRECISION, PARAMETER :: Pi_over_4 = 0.785398163397448D0 DOUBLE PRECISION, PARAMETER :: degrees_per_radian = 57.2957795130823D0 DOUBLE PRECISION, PARAMETER :: radians_per_degree = 0.0174532925199433D0 ! Relation of map-projection plane (x,y in meters) to the ! map window on the graphics page: LOGICAL :: mp_zoom_defined = .FALSE. REAL*8 :: mp_scale_denominator, & ! denominator of map scale (8.D8 ?). & mp_x_center_meters, & ! x of center of map, meters (0.D0 ?). & mp_y_center_meters, & ! y of center of map, meters (0.D0 ?). & 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.D0 ?). ! and convenience values, arrays and factors derived from these: REAL*8 :: mp_meters_per_point REAL*8, DIMENSION(2,2) :: mp_m_2_pts_matrix, mp_pts_2_m_matrix REAL*8, DIMENSION(2) :: mp_m_2_pts_vector, mp_pts_2_m_vector ! Memory of pen position (to avoid null segments in paths): REAL*8 :: mp_last_L3_x_meters, mp_last_L3_y_meters REAL*8, DIMENSION(3) :: mp_last_uvec REAL*8 :: 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 INTEGER :: mp_walls_in_use ! count of .TRUE. values in list above ! 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*8, DIMENSION(3) :: mp_bottom_uvec, mp_right_uvec, & & mp_top_uvec, mp_left_uvec LOGICAL :: mp_using_cut = .FALSE. ! Initialized F in case ! no map projection is ever defined, and we work ! only with projection-plane (x, y). Specific ! map projections will set to T later, and they ! must also then define values for the next 2: REAL*8, DIMENSION(3) :: mp_cutPole_uvec REAL*8 :: mp_intoCut_RC_azimuth_radians ! --------------------------------------------------------- ! | 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): !{N.B. Following DOUBLE PRECISION was present in original MODULE Map_Projections:} 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*8, 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*8, 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*8 :: mp_radius_meters = 6371000.D0 ! 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*8, 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*8 :: mp_proj2_F, mp_proj2_lambda0, mp_proj2_n, mp_proj2_rho0 REAL*8 :: mp_proj3_C, mp_proj3_lambda0, mp_proj3_n, mp_proj3_rho0 REAL*8 :: mp_proj4_lat0_radians, mp_proj4_lon0_radians REAL*8 :: 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*8 :: 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.D0,0.D0) 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*8, DIMENSION(2,2) :: mp_pp_2_xy_matrix, mp_xy_2_pp_matrix REAL*8, 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*8, DIMENSION(2,2) :: mp_upright, mp_downright !============= variables for saving state of module Map_Projections (mp_): LOGICAL :: mp_state_saved = .FALSE. ! until CALL DSave_mp_State() ! (The following are listed in same order as above, but without the comments) LOGICAL :: sp_zoom_defined REAL*8 :: sp_scale_denominator, & & sp_x_center_meters, & & sp_y_center_meters, & & sp_xy_wrt_page_radians REAL*8 :: sp_meters_per_point REAL*8, DIMENSION(2,2) :: sp_m_2_pts_matrix, sp_pts_2_m_matrix REAL*8, DIMENSION(2) :: sp_m_2_pts_vector, sp_pts_2_m_vector REAL*8 :: sp_last_L3_x_meters, sp_last_L3_y_meters REAL*8, DIMENSION(3) :: sp_last_uvec REAL*8 :: 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*8, DIMENSION(3) :: sp_bottom_uvec, sp_right_uvec, & & sp_top_uvec, sp_left_uvec REAL*8 :: sp_bottom_maxdot, sp_right_maxdot, sp_top_maxdot, sp_left_maxdot REAL*8, DIMENSION(3) :: sp_bottom_applies_uvec, sp_right_applies_uvec, & & sp_top_applies_uvec, sp_left_applies_uvec REAL*8, 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*8 :: sp_radius_meters = 6371000. REAL*8, DIMENSION(3) :: sp_pole_1_uvec, sp_pole_2_uvec, sp_pole_3_uvec, & & sp_projpoint_uvec REAL*8 :: sp_proj2_F, sp_proj2_lambda0, sp_proj2_n, sp_proj2_rho0 REAL*8 :: sp_proj3_C, sp_proj3_lambda0, sp_proj3_n, sp_proj3_rho0 REAL*8 :: sp_proj4_lat0_radians, sp_proj4_lon0_radians REAL*8 :: sp_proj5_lambda0, sp_proj5_lat0_radians, sp_proj5_n, sp_proj5_rtan, sp_proj5_ypole LOGICAL :: sp_user_xy REAL*8 :: sp_x_projpoint_meters, sp_y_projpoint_meters, sp_y_azimuth_radians REAL*8, DIMENSION(2,2) :: sp_pp_2_xy_matrix, sp_xy_2_pp_matrix REAL*8, DIMENSION(2) :: sp_pp_2_xy_vector, sp_xy_2_pp_vector REAL*8, DIMENSION(2,2) :: sp_upright, sp_downright !============= end of mp_ - saving variables ========================= CONTAINS ! =============================================== ! | USER ROUTINES | ! | (in the typical order of use) | ! =============================================== SUBROUTINE DSet_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.D6 ~ 1.D9). ! 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*8, INTENT(IN) :: scale_denominator REAL*8, INTENT(IN), OPTIONAL :: x_center_meters, & & y_center_meters, xy_wrt_page_radians REAL*8 :: factor ! Store values in module-level variables: mp_zoom_defined = .TRUE. mp_scale_denominator = scale_denominator mp_meters_per_point = scale_denominator / (72.D0 * 39.370079D0) IF (PRESENT(x_center_meters)) THEN mp_x_center_meters = x_center_meters ELSE mp_x_center_meters = 0.0D0 END IF IF (PRESENT(y_center_meters)) THEN mp_y_center_meters = y_center_meters ELSE mp_y_center_meters = 0.D0 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.D0 END IF ! Pre-compute meters -> points conversion matrix and vector: factor = 72.D0 * 39.370079D0 / mp_scale_denominator mp_m_2_pts_matrix(1,1) = factor * DCOS(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 * DSIN(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.D0 * 39.370079D0) mp_pts_2_m_matrix(1,1) = factor * DCOS(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 * DSIN(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 DSet_Zoom SUBROUTINE DNew_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*8, 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 = DNext_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.0D0 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 DNew_L3_Path with a path open.')") CALL DTraceback END IF ELSE ! zoom not defined WRITE (*,"(' ERROR: Must DSet_Zoom before DNew_L3_Path.')") CALL DTraceback END IF ELSE ! no page open WRITE (*,"(' ERROR: Cannot DNew_L3_Path before DBegin_Page.')") CALL DTraceback END IF END SUBROUTINE DNew_L3_Path SUBROUTINE DLine_To_L3 (x_meters, y_meters) IMPLICIT NONE REAL*8, 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.0D0 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 DTraceback END IF ELSE WRITE (*,"(' ERROR: Cannot DLine_To_L3 in path on level ',I2)") & & ai_current_path_level CALL DTraceback END IF ELSE WRITE (*,"(' ERROR: Cannot DLine_To_L3 before DNew_L3_Path.')") CALL DTraceback END IF END SUBROUTINE DLine_To_L3 SUBROUTINE DCurve_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*8, 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 DTraceback END IF ELSE WRITE (*,"(' ERROR: Cannot DCurve_To_L3 in path on level ',I2)") & & ai_current_path_level CALL DTraceback END IF ELSE WRITE (*,"(' ERROR: Cannot DCurve_To_L3 before DNew_L3_Path.')") CALL DTraceback END IF END SUBROUTINE DCurve_To_L3 SUBROUTINE DEnd_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 DLine_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.D25 ! (undefined) mp_last_L3_y_meters = -9.D25 CALL DProcess_L3_Paths ELSE !neither stroked nor filled WRITE (*,"(' ERROR: DEnd_L3_path with neither fill nor stroke.')") CALL DTraceback END IF ELSE WRITE (*,"(' ERROR: Cannot DEnd_L3_path when other level open:',I3)") level CALL DTraceback END IF ELSE WRITE (*,"(' ERROR: Cannot DEnd_L3_Path when path not open.')") CALL DTraceback END IF END SUBROUTINE DEnd_L3_Path SUBROUTINE DCircle_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*8, INTENT(IN) :: x_meters, y_meters, radius_meters LOGICAL, INTENT(IN) :: stroke, fill REAL*8 :: radius_points, x_points, y_points IF (.NOT.ai_page_open) THEN WRITE (*,"(' ERROR: Cannot DCircle_on_L3 before DBegin_Page.')") CALL DTraceback END IF IF (.NOT.mp_zoom_defined) THEN WRITE (*,"(' ERROR: Cannot DCircle_on_L3 before DSet_Zoom.')") CALL DTraceback END IF IF (ai_in_path) THEN WRITE (*,"(' ERROR: Cannot DCircle_on_L3 with another path already open.')") CALL DTraceback END IF IF (radius_meters < 0.0D0) THEN WRITE (*,"(' ERROR: Negative radius to DCircle_on_L3: ',1P,E9.2)") radius_meters CALL DTraceback END IF CALL DMeters_2_Points (x_meters,y_meters, x_points,y_points) radius_points = 2834.6D0 * (radius_meters / mp_scale_denominator) CALL DCircle_on_L12 (2, x_points,y_points, radius_points, stroke, fill) END SUBROUTINE DCircle_on_L3 SUBROUTINE DL3_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.D0 for left, 0.5D0 ! for middle, or 1.D0 for right-end. Values outside this range ! will also work, within reason. Note that alignment is ! most accurate at values of 0.D0, 0.5D0, or 1.D0, 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.D0 gives alignment with base of characters (normal); ! -0.4D0 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.0D0 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*8, INTENT(IN) :: x_meters, y_meters, angle_radians, & & lr_fraction, ud_fraction LOGICAL, INTENT(IN) :: from_x CHARACTER*(*),INTENT(IN) :: text INTEGER :: bytes REAL*8 :: 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 DProcess_L3_Text (x_meters, y_meters, page_angle_radians, & & font_points, lr_fraction, ud_fraction, & & text) END SUBROUTINE DL3_Text SUBROUTINE DL5_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.D0 for left, 0.5D0 ! for middle, or 1.D0 for right-end. Values outside this range ! will also work, within reason. Note that alignment is ! most accurate at values of 0.D0, 0.5D0, or 1.D0, 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.D0 gives alignment with base of characters (normal); ! -0.4D0 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.0D0 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*8, DIMENSION(3), INTENT(IN) :: uvec ! unit vector = position INTEGER, INTENT(IN) :: font_points REAL*8, 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 DProcess_L5_Text (uvec, angle_radians, from_east, & & font_points, lr_fraction, ud_fraction, & & text) END SUBROUTINE DL5_Text SUBROUTINE DL67_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.D0 for left, 0.5D0 ! for middle, or 1.D0 for right-end. Values outside this range ! will also work, within reason. Note that alignment is ! most accurate at values of 0.D0, 0.5D0, or 1.D0, 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.D0 gives alignment with base of characters (normal); ! -0.4D0 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.0D0 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*8, INTENT(IN) :: r1, r2 REAL*8, INTENT(IN) :: angle_radians, & & lr_fraction, ud_fraction LOGICAL, INTENT(IN) :: from_east CHARACTER*(*), INTENT(IN) :: text INTEGER :: bytes REAL*8 :: lat, lon, phi, theta bytes = LEN_TRIM(text) IF (bytes < 1) RETURN IF (level == 6) THEN theta = r1 phi = r2 IF ((theta < 0.D0).OR.(theta > Pi)) THEN WRITE (*,"(' ERROR: Illegal colatitude theta for DL67_Text.')") CALL DTraceback ELSE CALL DProcess_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.D0).OR.(lat < -90.D0)) THEN WRITE (*,"(' ERROR: Illegal latitude for DL67_Text.')") CALL DTraceback ELSE CALL DProcess_L7_Text (lon, lat, angle_radians, from_east, & & font_points, lr_fraction, ud_fraction, & & text) END IF ELSE WRITE (*,"(' ERROR: DL67_Text cannot accept level ',I2,'.')") level CALL DTraceback END IF END SUBROUTINE DL67_Text SUBROUTINE DSet_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.D0, 0.D0) at projpoint_uvec, with +y Northward ! and +x Eastward at that point. IMPLICIT NONE REAL*8, INTENT(IN) :: radius_meters REAL*8, DIMENSION(3), INTENT(IN) :: projpoint_uvec REAL*8, INTENT(IN) :: belt_azimuth_radians REAL*8, INTENT(IN), OPTIONAL :: x_projpoint_meters, y_projpoint_meters, & & y_azimuth_radians REAL*8 :: using_belt_azimuth_radians REAL*8, DIMENSION(3) :: omega_uvec, t_vec mp_projection_defined = .TRUE. mp_projection = 'Mercator' mp_projection_number = 1 mp_radius_meters = radius_meters CALL DMake_Uvec(projpoint_uvec, mp_projpoint_uvec) ! just for sure! ! find the 3 orthogonal uvec's used for Mercator projection: CALL DMake_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 (DSIN(belt_azimuth_radians) == 0.0D0) THEN using_belt_azimuth_radians = 0.0D0 ELSE ! use Eastward direction on great circle IF (DSIN(belt_azimuth_radians) > 0.0D0) THEN using_belt_azimuth_radians = belt_azimuth_radians ELSE using_belt_azimuth_radians = belt_azimuth_radians + Pi END IF END IF CALL DTurn_To (using_belt_azimuth_radians, mp_pole_1_uvec, Pi_over_2, & ! inputs & omega_uvec, mp_pole_2_uvec) CALL DMake_Uvec (mp_pole_2_uvec, mp_pole_2_uvec) ! just to be sure CALL DCross (mp_pole_1_uvec, mp_pole_2_uvec, mp_pole_3_uvec) CALL DMake_Uvec (mp_pole_3_uvec, mp_pole_3_uvec) ! just to be sure ! 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) = DSIN(using_belt_azimuth_radians) mp_upright(2,2) = mp_upright(1,1) mp_upright(2,1) = DCOS(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: mp_walls_in_use = 4 ! 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 = DSIN(69.98D0 / 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 DProcess_L4_Paths below 0.174533D0 (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.1D0 * radians_per_degree * mp_pole_1_uvec CALL DMake_Uvec(t_vec, mp_left_uvec) ! Notes: Because of wrap-around, it points "right". ! It is 0.2D0 degrees of "longitude" away ! from the corresponding "right" limit. t_vec = -mp_pole_2_uvec - 0.05D0 * radians_per_degree * mp_pole_1_uvec CALL DMake_Uvec(t_vec, mp_left_applies_uvec) mp_left_maxdot = 0.D0 ! cut is a great-circle arc ! Right ("East-longitude") side of the cut: mp_using_right = .TRUE. t_vec = -mp_pole_2_uvec - 0.1D0 * radians_per_degree * mp_pole_1_uvec CALL DMake_Uvec(t_vec, mp_right_uvec) ! Notes: Because of wrap-around, it points "right". ! It is 0.2D0 degrees of "longitude" away ! from the corresponding "right" limit. t_vec = mp_pole_2_uvec - 0.05D0 * radians_per_degree * mp_pole_1_uvec CALL DMake_Uvec(t_vec, mp_right_applies_uvec) mp_right_maxdot = 0.D0 ! cut is a great-circle arc ! Describe the "cut" edges of the rind (on the sphere), ! to help detect cut-spanning line segments later, in DProcess_L5_Paths: mp_using_cut = .TRUE. mp_cutPole_uvec = mp_pole_3_uvec mp_intoCut_RC_azimuth_radians = Pi + DRelative_Compass (from_uvec = mp_cutPole_uvec, & & to_uvec = projpoint_uvec) !Note that usage will be via DCOS and DSIN, so +- Two_Pi makes no difference. ! Complete description of rind by finding corners: CALL DCorners (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.D0 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.D0 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.D0 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) = DCOS(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) = DSIN(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) = DCOS(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) = DSIN(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 DSet_Mercator SUBROUTINE DSet_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.D0, 0.D0) at projpoint_uvec, with +y Northward ! and +x Eastward at that point. IMPLICIT NONE REAL*8, INTENT(IN) :: radius_meters, standard_parallel_gap_radians REAL*8, DIMENSION(3), INTENT(IN) :: projpoint_uvec, cone_pole_uvec REAL*8, INTENT(IN), OPTIONAL :: x_projpoint_meters, y_projpoint_meters, & & y_azimuth_radians REAL*8 :: latitude_0_radians, latitude_1_radians, latitude_2_radians, & & Pi_over_4, toconepole_azimuth_radians REAL*8, DIMENSION(3) :: t_vec mp_projection_defined = .TRUE. mp_projection = 'Lambert Conformal Conic' mp_projection_number = 2 mp_radius_meters = radius_meters CALL DMake_Uvec(projpoint_uvec, mp_projpoint_uvec) ! just for sure! IF (DArc(mp_projpoint_uvec, cone_pole_uvec) > 1.396D0) THEN WRITE (*,"(' ERROR: Distance from projection point to pole of cone exceeds 80 degrees.')") CALL DTraceback END IF ! find the 3 orthogonal uvec's used for all conic projections: CALL DMake_Uvec (cone_pole_uvec, mp_pole_3_uvec) ! just for insurance!!! CALL DCross (mp_pole_3_uvec, mp_projpoint_uvec, t_vec) CALL DMake_Uvec (t_vec, mp_pole_2_uvec) CALL DCross (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 - DArc(mp_projpoint_uvec, cone_pole_uvec) ! Note: "Latitude" is now with respect to conic_pole_uvec as "North pole"! Pi_over_4 = 0.5D0 * Pi_over_2 IF (standard_parallel_gap_radians > 0.0D0) THEN latitude_1_radians = latitude_0_radians - 0.5D0 * standard_parallel_gap_radians latitude_2_radians = latitude_0_radians + 0.5D0 * standard_parallel_gap_radians IF (latitude_2_radians >= Pi_over_2) THEN WRITE (*,"(' ERROR: Excessive gap between standard parallels.')") CALL DTraceback END IF mp_proj2_n = DLOG(DCOS(latitude_1_radians)/DCOS(latitude_2_radians)) / & & DLOG(DTAN(Pi_over_4 + latitude_2_radians/2.D0)/DTAN(Pi_over_4 + latitude_1_radians/2.D0)) mp_proj2_F = DCOS(latitude_1_radians) * DTAN(Pi_over_4 + latitude_1_radians/2.D0)**mp_proj2_n / mp_proj2_n ELSE ! only one standard parallel, assigned to be latitude_1 = latitude_0. mp_proj2_n = DSIN(latitude_0_radians) mp_proj2_F = DCOS(latitude_0_radians) * DTAN(Pi_over_4 + latitude_0_radians/2.D0)**mp_proj2_n / mp_proj2_n END IF mp_proj2_rho0 = mp_radius_meters * mp_proj2_F / & & DTAN(Pi_over_4 + latitude_0_radians/2.D0)**mp_proj2_n mp_proj2_lambda0 = 0.D0 ! 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 = DRelative_Compass(mp_projpoint_uvec, mp_pole_3_uvec) mp_upright(1,1) = DCOS(toconepole_azimuth_radians) mp_upright(2,2) = mp_upright(1,1) mp_upright(1,2) = DSIN(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) mp_walls_in_use = 0 ! 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_walls_in_use = mp_walls_in_use + 1 mp_top_uvec = mp_pole_3_uvec mp_top_applies_uvec = mp_top_uvec mp_top_maxdot = DSIN(latitude_0_radians + Pi_over_4) ! < 1.D0 END IF ! Low "latitude" (in the possibly-rotated coordinate system) limit: mp_using_bottom = .TRUE. mp_walls_in_use = mp_walls_in_use + 1 mp_bottom_uvec = -mp_pole_3_uvec mp_bottom_applies_uvec = mp_bottom_uvec mp_bottom_maxdot = MAX(0.D0, -DSIN(latitude_0_radians - Pi_over_4)) ! Left ("West-longitude") side of the cut: mp_using_left = .TRUE. mp_walls_in_use = mp_walls_in_use + 1 t_vec = mp_pole_2_uvec - 0.1D0 * radians_per_degree * mp_pole_1_uvec CALL DMake_Uvec(t_vec, mp_left_uvec) ! Notes: Because of wrap-around, it points "right". ! It is 0.2D0 degrees of "longitude" away ! from the corresponding "right" limit. t_vec = -mp_pole_2_uvec - 0.05D0 * radians_per_degree * mp_pole_1_uvec CALL DMake_Uvec(t_vec, mp_left_applies_uvec) mp_left_maxdot = 0.D0 ! cut is a great-circle arc ! Right ("East-longitude") side of the cut: mp_using_right = .TRUE. mp_walls_in_use = mp_walls_in_use + 1 t_vec = -mp_pole_2_uvec - 0.1D0 * radians_per_degree * mp_pole_1_uvec CALL DMake_Uvec(t_vec, mp_right_uvec) ! Notes: Because of wrap-around, it points "right". ! It is 0.2D0 degrees of "longitude" away ! from the corresponding "right" limit. t_vec = mp_pole_2_uvec - 0.05D0 * radians_per_degree * mp_pole_1_uvec CALL DMake_Uvec(t_vec, mp_right_applies_uvec) mp_right_maxdot = 0.D0 ! cut is a great-circle arc ! Describe the "cut" edges of the rind (on the sphere), ! to help detect cut-spanning line segments later, in DProcess_L5_Paths: mp_using_cut = .TRUE. mp_cutPole_uvec = cone_pole_uvec mp_intoCut_RC_azimuth_radians = Pi + DRelative_Compass (from_uvec = mp_cutPole_uvec, & & to_uvec = projpoint_uvec) !Note that usage will be via DCOS and DSIN, so +- Two_Pi makes no difference. ! Complete description of rind by finding corners: CALL DCorners (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.D0 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.D0 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.D0 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) = DCOS(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) = DSIN(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) = DCOS(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) = DSIN(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 DSet_Lambert_Conformal_Conic SUBROUTINE DSet_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.D0, 0.D0) at projpoint_uvec, with +y Northward ! and +x Eastward at that point. IMPLICIT NONE REAL*8, INTENT(IN) :: radius_meters, standard_parallel_gap_radians REAL*8, DIMENSION(3), INTENT(IN) :: projpoint_uvec, cone_pole_uvec REAL*8, INTENT(IN), OPTIONAL :: x_projpoint_meters, y_projpoint_meters, & & y_azimuth_radians REAL*8 :: latitude_0_radians, latitude_1_radians, latitude_2_radians, & & Pi_over_4, toconepole_azimuth_radians REAL*8, DIMENSION(3) :: t_vec mp_projection_defined = .TRUE. mp_projection = 'Albers Equal-Area Conic' mp_projection_number = 3 mp_radius_meters = radius_meters CALL DMake_Uvec(projpoint_uvec, mp_projpoint_uvec) ! just for sure! IF (DArc(mp_projpoint_uvec, cone_pole_uvec) > 1.396D0) THEN WRITE (*,"(' ERROR: Distance from projection point to pole of cone exceeds 80 degrees.')") CALL DTraceback END IF ! find the 3 orthogonal uvec's used for all conic projections: CALL DMake_Uvec (cone_pole_uvec, mp_pole_3_uvec) ! just for insurance!!! CALL DCross (mp_pole_3_uvec, mp_projpoint_uvec, t_vec) CALL DMake_Uvec (t_vec, mp_pole_2_uvec) CALL DCross (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 - DArc(mp_projpoint_uvec, cone_pole_uvec) ! Note: "Latitude" is now with respect to conic_pole_uvec as "North pole"! Pi_over_4 = 0.5D0 * Pi_over_2 IF (standard_parallel_gap_radians > 0.0D0) THEN latitude_1_radians = latitude_0_radians - 0.5D0 * standard_parallel_gap_radians latitude_2_radians = latitude_0_radians + 0.5D0 * standard_parallel_gap_radians IF (latitude_2_radians > Pi_over_2) THEN IF (DABS(latitude_2_radians - Pi_over_2) <= 0.01D0) THEN latitude_2_radians = Pi_over_2 ELSE WRITE (*,"(' ERROR: More poleward standard parallel has latitude over 90.')") CALL DTraceback END IF END IF IF (DABS(latitude_1_radians) < 0.01D0) THEN WRITE (*,"(' ERROR: Neither standard parallel may be on equator.')") CALL DTraceback END IF mp_proj3_n = (DSIN(latitude_1_radians) + DSIN(latitude_2_radians)) / 2.D0 mp_proj3_C = DCOS(latitude_1_radians)**2 + 2.D0 * mp_proj3_n * DSIN(latitude_1_radians) ELSE ! only one standard parallel, assumed to be latitude_1 = latitude_0 (projection point). mp_proj3_n = DSIN(latitude_0_radians) mp_proj3_C = DCOS(latitude_0_radians)**2 + 2.D0 * mp_proj3_n * DSIN(latitude_0_radians) END IF mp_proj3_rho0 = mp_radius_meters * DSQRT(mp_proj3_C - 2.D0 * mp_proj3_n * DSIN(latitude_0_radians)) / mp_proj3_n mp_proj3_lambda0 = 0.D0 ! 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 = DRelative_Compass(mp_projpoint_uvec, mp_pole_3_uvec) mp_upright(1,1) = DCOS(toconepole_azimuth_radians) mp_upright(2,2) = mp_upright(1,1) mp_upright(1,2) = DSIN(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) mp_walls_in_use = 0 ! 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_walls_in_use = mp_walls_in_use + 1 mp_top_uvec = mp_pole_3_uvec mp_top_applies_uvec = mp_top_uvec mp_top_maxdot = DSIN(MIN(1.4835D0,latitude_0_radians + Pi_over_4)) ! < 1.D0 END IF ! Low "latitude" (in the possibly-rotated coordinate system) limit: mp_using_bottom = .TRUE. mp_walls_in_use = mp_walls_in_use + 1 mp_bottom_uvec = -mp_pole_3_uvec mp_bottom_applies_uvec = mp_bottom_uvec mp_bottom_maxdot = MAX(0.D0, -DSIN(latitude_0_radians - Pi_over_4)) ! Left ("West-longitude") side of the cut: mp_using_left = .TRUE. mp_walls_in_use = mp_walls_in_use + 1 t_vec = mp_pole_2_uvec - 0.1D0 * radians_per_degree * mp_pole_1_uvec CALL DMake_Uvec(t_vec, mp_left_uvec) ! Notes: Because of wrap-around, it points "right". ! It is 0.2D0 degrees of "longitude" away ! from the corresponding "right" limit. t_vec = -mp_pole_2_uvec - 0.05D0 * radians_per_degree * mp_pole_1_uvec CALL DMake_Uvec(t_vec, mp_left_applies_uvec) mp_left_maxdot = 0.D0 ! cut is a great-circle arc ! Right ("East-longitude") side of the cut: mp_using_right = .TRUE. mp_walls_in_use = mp_walls_in_use + 1 t_vec = -mp_pole_2_uvec - 0.1D0 * radians_per_degree * mp_pole_1_uvec CALL DMake_Uvec(t_vec, mp_right_uvec) ! Notes: Because of wrap-around, it points "right". ! It is 0.2D0 degrees of "longitude" away ! from the corresponding "right" limit. t_vec = mp_pole_2_uvec - 0.05D0 * radians_per_degree * mp_pole_1_uvec CALL DMake_Uvec(t_vec, mp_right_applies_uvec) mp_right_maxdot = 0.D0 ! cut is a great-circle arc ! Describe the "cut" edges of the rind (on the sphere), ! to help detect cut-spanning line segments later, in DProcess_L5_Paths: mp_using_cut = .TRUE. mp_cutPole_uvec = cone_pole_uvec mp_intoCut_RC_azimuth_radians = Pi + DRelative_Compass (from_uvec = mp_cutPole_uvec, & & to_uvec = projpoint_uvec) !Note that usage will be via DCOS and DSIN, so +- Two_Pi makes no difference. ! Complete description of rind by finding corners: CALL DCorners (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.D0 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.D0 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.D0 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) = DCOS(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) = DSIN(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) = DCOS(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) = DSIN(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 DSet_Albers_Equal_Area_Conic SUBROUTINE DSet_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.D0, 0.D0) at projpoint_uvec, with +y Northward ! and +x Eastward at that point. IMPLICIT NONE REAL*8, INTENT(IN) :: radius_meters REAL*8, DIMENSION(3), INTENT(IN) :: projpoint_uvec, cone_pole_uvec REAL*8, INTENT(IN), OPTIONAL :: x_projpoint_meters, y_projpoint_meters, & & y_azimuth_radians REAL*8 :: toconepole_azimuth_radians REAL*8, DIMENSION(3) :: t_vec mp_projection_defined = .TRUE. mp_projection = 'Polyconic' mp_projection_number = 4 mp_radius_meters = radius_meters CALL DMake_Uvec(projpoint_uvec, mp_projpoint_uvec) ! just for sure! IF (DArc(mp_projpoint_uvec, cone_pole_uvec) > 1.396D0) THEN WRITE (*,"(' ERROR: Distance from projection point to pole of cone exceeds 80 degrees.')") CALL DTraceback END IF ! find the 3 orthogonal uvec's used for all conic projections: CALL DMake_Uvec (cone_pole_uvec, mp_pole_3_uvec) ! just for insurance!!! CALL DCross (mp_pole_3_uvec, mp_projpoint_uvec, t_vec) CALL DMake_Uvec (t_vec, mp_pole_2_uvec) CALL DCross (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 - DArc(mp_projpoint_uvec, cone_pole_uvec) ! Note: "Latitude" is now with respect to conic_pole_uvec as "North pole"! mp_proj4_lon0_radians = 0.D0 ! 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 = DRelative_Compass(mp_projpoint_uvec, mp_pole_3_uvec) mp_upright(1,1) = DCOS(toconepole_azimuth_radians) mp_upright(2,2) = mp_upright(1,1) mp_upright(1,2) = DSIN(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) mp_walls_in_use = 3 ! 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.036D0 * mp_pole_1_uvec ! 2 degree tilt CALL DMake_Uvec (t_vec, mp_bottom_uvec) mp_bottom_applies_uvec = mp_bottom_uvec mp_bottom_maxdot = 0.D0 ! 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 = -DSIN(80.01D0*radians_per_degree)*mp_pole_1_uvec & & -DCOS(80.01D0*radians_per_degree)*mp_pole_2_uvec CALL DMake_Uvec(t_vec, mp_left_uvec) mp_left_applies_uvec = mp_left_uvec mp_left_maxdot = 0.D0 ! cut is a great-circle arc ! Right ("East-longitude") side of the cut: mp_using_right = .TRUE. t_vec = -DSIN(80.01D0*radians_per_degree)*mp_pole_1_uvec & & +DCOS(80.01D0*radians_per_degree)*mp_pole_2_uvec CALL DMake_Uvec(t_vec, mp_right_uvec) mp_right_applies_uvec = mp_right_uvec mp_right_maxdot = 0.D0 ! cut is a great-circle arc ! Describe the "cut" edges of the rind (on the sphere), ! to help detect cut-spanning line segments later, in DProcess_L5_Paths: mp_using_cut = .TRUE. mp_cutPole_uvec = cone_pole_uvec mp_intoCut_RC_azimuth_radians = Pi + DRelative_Compass (from_uvec = mp_cutPole_uvec, & & to_uvec = projpoint_uvec) !Note that usage will be via DCOS and DSIN, so +- Two_Pi makes no difference. ! Complete description of rind by finding corners: CALL DCorners (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.D0 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.D0 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.D0 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) = DCOS(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) = DSIN(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) = DCOS(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) = DSIN(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 DSet_Polyconic SUBROUTINE DSet_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.D0, 0.D0) at projpoint_uvec, with +y Northward ! and +x Eastward at that point. IMPLICIT NONE REAL*8, INTENT(IN) :: radius_meters REAL*8, DIMENSION(3), INTENT(IN) :: projpoint_uvec, cone_pole_uvec REAL*8, INTENT(IN), OPTIONAL :: x_projpoint_meters, y_projpoint_meters, & & y_azimuth_radians REAL*8 :: Pi_over_4, toconepole_azimuth_radians REAL*8, DIMENSION(3) :: t_vec mp_projection_defined = .TRUE. mp_projection = 'Geometric Conic' mp_projection_number = 5 mp_radius_meters = radius_meters CALL DMake_Uvec(projpoint_uvec, mp_projpoint_uvec) ! just for sure! IF (DArc(mp_projpoint_uvec, cone_pole_uvec) > 1.396D0) THEN WRITE (*,"(' ERROR: Distance from projection point to pole of cone exceeds 80 degrees.')") CALL DTraceback END IF ! find the 3 orthogonal uvec's used for all conic projections: CALL DMake_Uvec (cone_pole_uvec, mp_pole_3_uvec) ! just for insurance!!! CALL DCross (mp_pole_3_uvec, mp_projpoint_uvec, t_vec) CALL DMake_Uvec (t_vec, mp_pole_2_uvec) CALL DCross (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.D0 mp_proj5_lambda0 = 0.D0 ! by definition of rotated coordinate system mp_proj5_lat0_radians = Pi_over_2 - DArc(mp_projpoint_uvec, mp_pole_3_uvec) mp_proj5_n = DSIN(mp_proj5_lat0_radians) ! ratio of "longitude" angles in unrolled cone to sphere mp_proj5_rtan = mp_radius_meters * DTAN(DArc(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 = DRelative_Compass(mp_projpoint_uvec, mp_pole_3_uvec) mp_upright(1,1) = DCOS(toconepole_azimuth_radians) mp_upright(2,2) = mp_upright(1,1) mp_upright(1,2) = DSIN(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) mp_walls_in_use = 0 ! 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_walls_in_use = mp_walls_in_use + 1 mp_top_uvec = mp_pole_3_uvec mp_top_applies_uvec = mp_top_uvec mp_top_maxdot = DSIN(mp_proj5_lat0_radians + Pi_over_4) ! < 1.D0 END IF ! Low "latitude" (in the possibly-rotated coordinate system) limit: mp_using_bottom = .TRUE. mp_walls_in_use = mp_walls_in_use + 1 mp_bottom_uvec = -mp_pole_3_uvec mp_bottom_applies_uvec = mp_bottom_uvec mp_bottom_maxdot = MAX(0.D0, -DSIN(mp_proj5_lat0_radians - Pi_over_4)) ! Left ("West-longitude") side of the cut: mp_using_left = .TRUE. mp_walls_in_use = mp_walls_in_use + 1 t_vec = mp_pole_2_uvec - 0.1D0 * radians_per_degree * mp_pole_1_uvec CALL DMake_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.05D0 * radians_per_degree * mp_pole_1_uvec CALL DMake_Uvec(t_vec, mp_left_applies_uvec) mp_left_maxdot = 0.D0 ! cut is a great-circle arc ! Right ("East-longitude") side of the cut: mp_using_right = .TRUE. mp_walls_in_use = mp_walls_in_use + 1 t_vec = -mp_pole_2_uvec - 0.1D0 * radians_per_degree * mp_pole_1_uvec CALL DMake_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.05D0 * radians_per_degree * mp_pole_1_uvec CALL DMake_Uvec(t_vec, mp_right_applies_uvec) mp_right_maxdot = 0.D0 ! cut is a great-circle arc ! Complete description of rind by finding corners: CALL DCorners (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.D0 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.D0 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.D0 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) = DCOS(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) = DSIN(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) = DCOS(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) = DSIN(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 DSet_Geometric_Conic SUBROUTINE DSet_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.D0, 0.D0) at projpoint_uvec, with +y Northward ! and +x Eastward at that point. IMPLICIT NONE REAL*8, INTENT(IN) :: radius_meters REAL*8, DIMENSION(3), INTENT(IN) :: projpoint_uvec REAL*8, INTENT(IN), OPTIONAL :: x_projpoint_meters, y_projpoint_meters, & & y_azimuth_radians REAL*8, DIMENSION(3) :: east_uvec, north_uvec mp_projection_defined = .TRUE. mp_projection = 'Stereographic' mp_projection_number = 6 mp_radius_meters = radius_meters CALL DMake_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 DNorthEast_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_walls_in_use = 1 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.001D0 ! Complete description of rind by finding corners: CALL DCorners (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.D0 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.D0 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.D0 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) = DCOS(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) = DSIN(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) = DCOS(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) = DSIN(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 DSet_Stereographic SUBROUTINE DSet_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.D0, 0.D0) at projpoint_uvec, with +y Northward ! and +x Eastward at that point. IMPLICIT NONE REAL*8, INTENT(IN) :: radius_meters REAL*8, DIMENSION(3), INTENT(IN) :: projpoint_uvec REAL*8, INTENT(IN), OPTIONAL :: x_projpoint_meters, y_projpoint_meters, & & y_azimuth_radians REAL*8, 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 DMake_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 DNorthEast_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_walls_in_use = 1 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.001D0 ! Complete description of rind by finding corners: CALL DCorners (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.D0 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.D0 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.D0 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) = DCOS(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) = DSIN(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) = DCOS(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) = DSIN(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 DSet_Lambert_Azimuthal_EqualArea SUBROUTINE DSet_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.D0, 0.D0) at projpoint_uvec, with +y Northward ! and +x Eastward at that point. IMPLICIT NONE REAL*8, INTENT(IN) :: radius_meters REAL*8, DIMENSION(3), INTENT(IN) :: projpoint_uvec REAL*8, INTENT(IN), OPTIONAL :: x_projpoint_meters, y_projpoint_meters, & & y_azimuth_radians REAL*8, DIMENSION(3) :: east_uvec, north_uvec mp_projection_defined = .TRUE. mp_projection = 'Azimuthal Equidistant' mp_projection_number = 8 mp_radius_meters = radius_meters CALL DMake_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 DNorthEast_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_walls_in_use = 1 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.001D0 ! Complete description of rind by finding corners: CALL DCorners (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.D0 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.D0 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.D0 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) = DCOS(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) = DSIN(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) = DCOS(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) = DSIN(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 DSet_Azimuthal_Equidistant SUBROUTINE DSet_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.D0, 0.D0) at projpoint_uvec, with +y Northward ! and +x Eastward at that point. IMPLICIT NONE REAL*8, INTENT(IN) :: radius_meters REAL*8, DIMENSION(3), INTENT(IN) :: projpoint_uvec REAL*8, INTENT(IN), OPTIONAL :: x_projpoint_meters, y_projpoint_meters, & & y_azimuth_radians REAL*8, DIMENSION(3) :: east_uvec, north_uvec mp_projection_defined = .TRUE. mp_projection = 'Orthographic' mp_projection_number = 9 mp_radius_meters = radius_meters CALL DMake_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 DNorthEast_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_walls_in_use = 1 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.01D0 ! (showing a tiny, hazy, over-the-horizon continuation to improve graticules) ! Complete description of rind by finding corners: CALL DCorners (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.D0 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.D0 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.D0 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) = DCOS(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) = DSIN(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) = DCOS(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) = DSIN(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 DSet_Orthographic SUBROUTINE DSet_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.D0, 0.D0) at projpoint_uvec, with +y Northward ! and +x Eastward at that point. IMPLICIT NONE REAL*8, INTENT(IN) :: radius_meters REAL*8, DIMENSION(3), INTENT(IN) :: projpoint_uvec REAL*8, INTENT(IN), OPTIONAL :: x_projpoint_meters, y_projpoint_meters, & & y_azimuth_radians REAL*8, 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 DMake_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 DNorthEast_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 get 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_walls_in_use = 4 mp_using_top = .TRUE. t_vec = -DSIN(70.D0*radians_per_degree) * mp_pole_1_uvec & & +DCOS(70.D0*radians_per_degree) * mp_pole_3_uvec CALL DMake_Uvec (t_vec, mp_top_uvec) mp_top_applies_uvec = mp_top_uvec mp_top_maxdot = 0.0D0 mp_using_left = .TRUE. t_vec = -DSIN(70.D0*radians_per_degree) * mp_pole_1_uvec & & -DCOS(70.D0*radians_per_degree) * mp_pole_2_uvec CALL DMake_Uvec (t_vec, mp_left_uvec) mp_left_applies_uvec = mp_left_uvec mp_left_maxdot = 0.0D0 mp_using_right = .TRUE. t_vec = -DSIN(70.D0*radians_per_degree) * mp_pole_1_uvec & & +DCOS(70.D0*radians_per_degree) * mp_pole_2_uvec CALL DMake_Uvec (t_vec, mp_right_uvec) mp_right_applies_uvec = mp_right_uvec mp_right_maxdot = 0.0D0 mp_using_bottom = .TRUE. t_vec = -DSIN(70.D0*radians_per_degree) * mp_pole_1_uvec & & -DCOS(70.D0*radians_per_degree) * mp_pole_3_uvec CALL DMake_Uvec (t_vec, mp_bottom_uvec) mp_bottom_applies_uvec = mp_bottom_uvec mp_bottom_maxdot = 0.0D0 ! Complete description of rind by finding corners: CALL DCorners (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.D0 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.D0 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.D0 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) = DCOS(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) = DSIN(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) = DCOS(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) = DSIN(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 DSet_Gnomonic SUBROUTINE DNew_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*8, 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 = DNext_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.0D0 ai_pathlib(1:3,0,path) = uvec mp_last_uvec = uvec ! (memory) ELSE ! ai_in_path already WRITE (*,"(' ERROR: Cannot DNew_L45_Path with a path open.')") CALL DTraceback END IF ELSE ! projection not defined yet WRITE (*,"(' ERROR: Must set map projection before DNew_L45_Path.')") CALL DTraceback END IF ELSE ! zoom from L3 <--> L2 not defined WRITE (*,"(' ERROR: Must DSet_Zoom before DNew_L45_Path.')") CALL DTraceback END IF ELSE ! no page open WRITE (*,"(' ERROR: Cannot DNew_L45_Path before DBegin_Page.')") CALL DTraceback END IF ELSE WRITE (*,"(' ERROR: Level ',I2,' not allowed in DNew_L45_Path.')") level CALL DTraceback END IF END SUBROUTINE DNew_L45_Path SUBROUTINE DGreat_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*8, 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 DGreat_To_L45 to antipode.')") CALL DTraceback ELSE ! USUAL CASE: next = ai_segments(path) + 1 ai_segments(path) = next ai_pathlib(4:6,next,path) = 0.0D0 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 DTraceback END IF ELSE WRITE (*,"(' ERROR: Cannot DGreat_To_L45 in path on level ',I2)") & & ai_current_path_level CALL DTraceback END IF ELSE WRITE (*,"(' ERROR: Cannot DGreat_To_L45 before DNew_L45_Path.')") CALL DTraceback END IF END SUBROUTINE DGreat_To_L45 SUBROUTINE DSafe_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 Small_To_L45, this routine will IGNORE any ! arc from the current uvec to itself, because it will be ! treated as a zero-length arc. (Such zero-length arcs ! are sometimes accidentally created during contouring.) ! If you WANT to create a complete small circle by setting ! the to_uvec == from_uvec, then CALL DSmall_To_L45 instead. IMPLICIT NONE REAL*8, DIMENSION(3), INTENT(IN) :: pole_uvec, to_uvec INTEGER :: next, path REAL*8, PARAMETER :: tolerance = 0.04D0 ! 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 (& GCMT) to be !unacceptably inconsistent, due to rounding of trend and plunge to nearest degree. REAL*8 :: arc1, arc2 REAL*8, 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 to debugger arc1 = DArc( pole_uvec, from_uvec ) arc2 = DArc( pole_uvec, to_uvec ) IF (ABS(arc1 - arc2) <= tolerance) THEN ! USUAL CASE: Now, just check for positive length: IF ((to_uvec(1) /= from_uvec(1)).OR.(to_uvec(2) /= from_uvec(2)).OR.(to_uvec(3) /= from_uvec(3))) THEN 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) END IF ! positive length, because to_uvec /= from_uvec ELSE ! inconsistent data; not a small circle! WRITE (*,"(' ERROR: Initial and final radii of small circle in DSmall_To_L45' & &/' are ',F6.4,' and ',F6.4,', which disagree by more than' & &/' the permitted tolerance of ',F6.4,' radians.')") arc1, arc2, tolerance CALL DTraceback END IF ELSE WRITE (*,"(' ERROR: Increase ai_longest from current ',I6)") ai_longest CALL DTraceback END IF ELSE WRITE (*,"(' ERROR: Cannot DSmall_To_L45 in path on level ',I2)") & & ai_current_path_level CALL DTraceback END IF ELSE WRITE (*,"(' ERROR: Cannot DSmall_To_L45 before DNew_L45_Path.')") CALL DTraceback END IF END SUBROUTINE DSafe_Small_To_L45 SUBROUTINE DSmall_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*8, DIMENSION(3), INTENT(IN) :: through_uvec, to_uvec INTEGER :: next, path REAL*8, PARAMETER :: tolerance = 0.001D0 ! limit for small circle arc, in radians (= 6.4 km). REAL*8 :: arc1 REAL*8, 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 DSmall_Through_L45 back to present location; nonunique.')") CALL DTraceback END IF END IF END IF !test for to_uvec very close to from_uvec: vec = from_uvec - to_uvec arc1 = DMagnitude(vec) IF (arc1 < tolerance) THEN CALL DGreat_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 DGreat_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 DGreat_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 DMake_Uvec (vec, uvec1) vec(1:3) = to_uvec(1:3) - through_uvec(1:3) CALL DMake_Uvec (vec, uvec2) CALL DCross (uvec1, uvec2, vec) CALL DMake_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 DTraceback END IF ELSE WRITE (*,"(' ERROR: Cannot DSmall_Through_L45 in path on level ',I2)") & & ai_current_path_level CALL DTraceback END IF ELSE WRITE (*,"(' ERROR: Cannot DSmall_Through_L45 before DNew_L45_Path.')") CALL DTraceback END IF END SUBROUTINE DSmall_Through_L45 SUBROUTINE DSmall_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. ! If you DON'T WANT THIS FEATURE, but prefer that small-circle ! arcs from a uvec to itself should be treated as zero-length ! (and therefore ignored) then CALL DSafe_Small_to_L45. IMPLICIT NONE REAL*8, DIMENSION(3), INTENT(IN) :: pole_uvec, to_uvec INTEGER :: next, path REAL*8, PARAMETER :: tolerance = 0.04D0 ! 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 (& GCMT) to be !unacceptably inconsistent, due to rounding of trend and plunge to nearest degree. REAL*8 :: arc1, arc2 REAL*8, 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 = DArc( pole_uvec, from_uvec ) arc2 = DArc( 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 DSmall_To_L45' & &/' are ',F6.4,' and ',F6.4,', which disagree by more than' & &/' the permitted tolerance of ',F6.4,' radians.')") arc1, arc2, tolerance CALL DTraceback END IF ELSE WRITE (*,"(' ERROR: Increase ai_longest from current ',I6)") ai_longest CALL DTraceback END IF ELSE WRITE (*,"(' ERROR: Cannot DSmall_To_L45 in path on level ',I2)") & & ai_current_path_level CALL DTraceback END IF ELSE WRITE (*,"(' ERROR: Cannot DSmall_To_L45 before DNew_L45_Path.')") CALL DTraceback END IF END SUBROUTINE DSmall_To_L45 SUBROUTINE DEnd_L45_Path (close, stroke, fill, retro) IMPLICIT NONE LOGICAL, INTENT(IN) :: close, stroke, fill LOGICAL, INTENT(IN), OPTIONAL :: retro ! If present and .TRUE., calls DOld_Complex_Process_L5_Paths() instead of DProcess_L5_Paths. INTEGER :: level, path REAL*8, 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: ai_closed(path) = close ai_stroked(path) = stroke ai_filled(path) = fill ai_in_path = .FALSE. mp_last_uvec = (/0.D0,0.D0,0.D0/) ! (undefined) IF (level == 4) THEN CALL DProcess_L4_Paths ELSE ! level 5; decide how to handle this... IF (PRESENT(retro)) THEN IF (retro) THEN CALL DOld_Complex_Process_L5_Paths() ELSE ! OPTIONAL argument "retro" is present but .FALSE. CALL DProcess_L5_Paths() END IF ELSE ! OPTIONAL argument "retro" is not present CALL DProcess_L5_Paths() END IF END IF ELSE !neither stroked nor filled WRITE (*,"(' ERROR: DEnd_L45_path with neither fill nor stroke.')") CALL DTraceback END IF ELSE WRITE (*,"(' ERROR: Cannot DEnd_L45_path when other level open.')") CALL DTraceback END IF ELSE WRITE (*,"(' ERROR: Cannot DEnd_L45_Path when path not open.')") CALL DTraceback END IF END SUBROUTINE DEnd_L45_Path SUBROUTINE DNew_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*8, INTENT(IN) :: r1, r2 INTEGER :: path REAL*8 :: 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 = DNext_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.0D0 IF (level == 6) THEN theta = r1 phi = r2 IF ((theta < 0.D0).OR.(theta > Pi)) THEN WRITE (*,"(' ERROR: Illegal colatitude theta for DNew_L67_Path.')") CALL DTraceback ELSE ai_pathlib(1,0,path) = theta ai_pathlib(2,0,path) = phi mp_last_theta = theta mp_last_phi = phi CALL DThetaPhi_2_Uvec (theta, phi, mp_last_uvec) END IF ELSE ! level == 7 lon = r1 lat = r2 IF ((lat > 90.D0).OR.(lat < -90.D0)) THEN WRITE (*,"(' ERROR: Illegal latitude theta for DNew_L67_Path.')") CALL DTraceback ELSE ai_pathlib(1,0,path) = lon ai_pathlib(2,0,path) = lat mp_last_lon = lon mp_last_lat = lat CALL DLonLat_2_Uvec (lon, lat, mp_last_uvec) END IF END IF ELSE ! ai_in_path already WRITE (*,"(' ERROR: Cannot DNew_L67_Path with a path open.')") CALL DTraceback END IF ELSE ! projection not defined yet WRITE (*,"(' ERROR: Must set map projection before DNew_L67_Path.')") CALL DTraceback END IF ELSE ! zoom from L3 <--> L2 not defined WRITE (*,"(' ERROR: Must DSet_Zoom before DNew_L67_Path.')") CALL DTraceback END IF ELSE ! no page open WRITE (*,"(' ERROR: Cannot DNew_L67_Path before DBegin_Page.')") CALL DTraceback END IF ELSE WRITE (*,"(' ERROR: Level ',I2,' not allowed in DNew_L67_Path.')") level CALL DTraceback END IF END SUBROUTINE DNew_L67_Path SUBROUTINE DGreat_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*8, INTENT(IN) :: r1, r2 INTEGER :: level, next, path REAL*8 :: lat, lon, phi, theta REAL*8, 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.D0).OR.(theta > Pi)) THEN WRITE (*,"(' ERROR: Illegal colatitude to DGreat_To_L67.')") CALL DTraceback 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 DThetaPhi_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 DGreat_To_L67 to antipode.')") CALL DTraceback ELSE ! USUAL CASE: next = ai_segments(path) + 1 ai_segments(path) = next ai_pathlib(3:6,next,path) = 0.0D0 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.D0).OR.(lat < -90.D0)) THEN WRITE (*,"(' ERROR: Illegal latitude to DGreat_To_L67.')") CALL DTraceback 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 DLonLat_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 DGreat_To_L67 to antipode.')") CALL DTraceback ELSE ! USUAL CASE: next = ai_segments(path) + 1 ai_segments(path) = next ai_pathlib(3:6,next,path) = 0.0D0 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 DTraceback END IF ELSE WRITE (*,"(' ERROR: Cannot DGreat_To_L67 in path on level ',I2)") & & ai_current_path_level CALL DTraceback END IF ELSE WRITE (*,"(' ERROR: Cannot DGreat_To_L67 before DNew_L67_Path.')") CALL DTraceback END IF END SUBROUTINE DGreat_To_L67 SUBROUTINE DSmall_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*8, INTENT(IN) :: r1, r2, r3, r4 INTEGER :: level, next, path REAL*8 :: lat, lon, phi, pole_lat, pole_lon, pole_phi, pole_theta, theta REAL*8, 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.D0).OR.(pole_theta > Pi)) THEN WRITE (*,"(' ERROR: Illegal pole colatitude to DSmall_To_L67.')") CALL DTraceback END IF pole_phi = r2 theta = r3 IF ((theta < 0.D0).OR.(theta > Pi)) THEN WRITE (*,"(' ERROR: Illegal colatitude to DSmall_To_L67.')") CALL DTraceback END IF phi = r4 next = ai_segments(path) + 1 ai_segments(path) = next ai_pathlib(5:6,next,path) = 0.0D0 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 DThetaPhi_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.D0).OR.(pole_lat < -90.D0)) THEN WRITE (*,"(' ERROR: Illegal pole latitude to DSmall_To_L67.')") CALL DTraceback END IF lon = r3 lat = r4 IF ((lat > 90.D0).OR.(lat < -90.D0)) THEN WRITE (*,"(' ERROR: Illegal latitude to DSmall_To_L67.')") CALL DTraceback END IF next = ai_segments(path) + 1 ai_segments(path) = next ai_pathlib(5:6,next,path) = 0.0D0 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 DLonLat_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 DTraceback END IF ELSE WRITE (*,"(' ERROR: Cannot DSmall_To_L67 in path on level ',I2)") & & ai_current_path_level CALL DTraceback END IF ELSE WRITE (*,"(' ERROR: Cannot DSmall_To_L67 before DNew_L67_Path.')") CALL DTraceback END IF END SUBROUTINE DSmall_To_L67 SUBROUTINE DEnd_L67_Path (close, stroke, fill, retro) IMPLICIT NONE LOGICAL, INTENT(IN) :: close, stroke, fill LOGICAL, INTENT(IN), OPTIONAL :: retro 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 DGreat_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 DGreat_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 (PRESENT(retro)) THEN IF (level == 6) THEN CALL DProcess_L6_Paths(retro) ELSE ! level 7 CALL DProcess_L7_Paths(retro) END IF ELSE IF (level == 6) THEN CALL DProcess_L6_Paths ELSE ! level 7 CALL DProcess_L7_Paths END IF END IF ELSE !neither stroked nor filled WRITE (*,"(' ERROR: DEnd_L67_path with neither fill nor stroke.')") CALL DTraceback END IF ELSE WRITE (*,"(' ERROR: Cannot DEnd_L67_path when other level open.')") CALL DTraceback END IF ELSE WRITE (*,"(' ERROR: Cannot DEnd_L67_Path when path not open.')") CALL DTraceback END IF END SUBROUTINE DEnd_L67_Path ! =============================================== ! | INTERNAL UTILITY ROUTINES | ! | (in alphabetical order ) | ! =============================================== REAL*8 FUNCTION DArc (from_uvec, to_uvec) ! Returns length of great-circle arc in radians. IMPLICIT NONE REAL*8, DIMENSION(3), INTENT(IN) :: from_uvec, to_uvec REAL*8 :: crossed, dotted REAL*8, DIMENSION(3) :: t_vec CALL DCross (from_uvec, to_uvec, t_vec) crossed = DLength(t_vec) ! >= 0.D0 dotted = DDot(from_uvec, to_uvec) DArc = DATAN2(crossed, dotted) ! 0.D0 to Pi END FUNCTION DArc SUBROUTINE DCircles_Intersect (pole_a_uvec, dot_a, first_a_uvec, last_a_uvec, & & pole_b_uvec, dot_b, first_b_uvec, last_b_uvec, & ! input & overlap, number, point1_uvec, point2_uvec) ! output ! Finds all the points of intersection (0, 1, or 2) between ! arcs of two small circles on a unit sphere. ! Note: The set of small circles includes great circles ! as a special case. The set of arcs includes complete- ! circle arcs as a special case. ! The two arcs are generically named "a" and "b". ! The pole of each is given by a "pole_x_uvec" (unit vector, ! from center of sphere; see module data for definitions). ! The plane containing the small circle is specified by "dot_x", ! its distance from the origin. (If either distance is ! greater than 1.00, there can be no points of intersection.) !"First" and "last" points on each arc are specified by uvecs ! (unit vectors) "first_x_uvec" and "last_x_uvec", ! where one goes counterclockwise around the pole ! first to last (assuming a viewpoint outside the ! unit sphere). If "first" = "last", the arc is a full circle. ! NOTE: Although "first_x_uvec" and "last_x_uvec" will usually ! be points on the corresponding arcs, they need not be. ! Only their azimuths from "pole_x_uvec" are used. In fact, ! if "first_x_uvec == last_x_uvec", then even the azimuths are ! not used, so ANY vector (even a zero vector) can be sent in ! these two positions to signal a complete small circle. ! CAUTION: Results will be strange and unpredictable if the ! the "unit vectors" supplied are not actually 1.00D0 long! ! Values returned are "number" (the count of intersection points) ! and unit vectors for as many points as were found. ! If only one intersection is found, it is always placed in point1_uvec. ! In the special case where the two small circles are the same ! circle, and share some common arc, the logical flag "overlap" ! is set, and the first and last points of the common arc are ! reported, and number = 2 on output. If both of these overlapped ! circles are complete circles, these first and last points are ! the same point. IMPLICIT NONE REAL*8, DIMENSION(3), INTENT(IN) :: pole_a_uvec, first_a_uvec, last_a_uvec, & & pole_b_uvec, first_b_uvec, last_b_uvec !{Following DOUBLE PRECISION was present in the original MODULE Map_Projections.} DOUBLE PRECISION, INTENT(IN) :: dot_a, dot_b LOGICAL, INTENT(OUT) :: overlap INTEGER, INTENT(OUT) :: number REAL*8, DIMENSION(3), INTENT(OUT) :: point1_uvec, point2_uvec !{Following DOUBLE PRECISIONs were present in the original MODULE Map_Projections.} DOUBLE PRECISION :: a_dot_b, alpha, beta, determinant, t DOUBLE PRECISION, DIMENSION(2,2) :: inverse LOGICAL :: ring_a, ring_b ! are the arcs actually complete circles? REAL*8 :: along, antipole_gap, & & first_a_radians, first_b_radians, last_a_radians, & & last_b_radians, min_radius, point_wrt_a_radians, point_wrt_b_radians, & & pole_gap REAL*8, PARAMETER :: tolerance = 0.0000000014D0 REAL*8, DIMENSION(3) :: line_uvec, offset, t_vec ! Check for an easy answer (no intersections): number = 0 overlap = .FALSE. IF (dot_a >= 1.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 = -DRelative_Compass(pole_a_uvec, first_a_uvec) last_a_radians = -DRelative_Compass(pole_a_uvec, last_a_uvec) IF (first_a_radians > last_a_radians) last_a_radians = last_a_radians + Two_Pi END IF IF (.NOT.ring_b) THEN first_b_radians = -DRelative_Compass(pole_b_uvec, first_b_uvec) last_b_radians = -DRelative_Compass(pole_b_uvec, last_b_uvec) IF (first_b_radians > last_b_radians) last_b_radians = last_b_radians + Two_Pi END IF ! Test for special cases of parallel and antipodal poles: t_vec = pole_a_uvec - pole_b_uvec pole_gap = DLength(t_vec) t_vec = pole_a_uvec + pole_b_uvec antipole_gap = DLength(t_vec) IF (pole_gap <= tolerance) THEN ! poles virtually identical IF (dot_a /= dot_b) RETURN ! no intersection ! From here on, assume dot_a == dot_b: IF (ring_a) THEN overlap = .TRUE. number = 2 point1_uvec = first_b_uvec point2_uvec = last_b_uvec ELSE IF (ring_b) THEN overlap = .TRUE. number = 2 point1_uvec = first_a_uvec point2_uvec = last_a_uvec ELSE ! neither circle is complete IF ((first_b_radians >= first_a_radians).AND. & &(first_b_radians <= last_a_radians)) THEN ! 1st of b is in a. point1_uvec = first_b_uvec IF (last_b_radians <= last_a_radians) THEN point2_uvec = last_b_uvec ELSE point2_uvec = last_a_uvec END IF IF ((point1_uvec(1) == point2_uvec(1)).AND. & &(point1_uvec(2) == point2_uvec(2)).AND. & &(point1_uvec(3) == point2_uvec(3))) THEN ! only one common point overlap = .FALSE. number = 1 ELSE ! a common arc overlap = .TRUE. number = 2 END IF ELSE IF ((first_a_radians >= first_b_radians).AND. & &(first_a_radians <= last_b_radians)) THEN ! 1st of a is in b. point1_uvec = first_a_uvec IF (last_a_radians <= last_b_radians) THEN point2_uvec = last_a_uvec ELSE point2_uvec = last_b_uvec END IF IF ((point1_uvec(1) == point2_uvec(1)).AND. & &(point1_uvec(2) == point2_uvec(2)).AND. & &(point1_uvec(3) == point2_uvec(3))) THEN ! only one common point overlap = .FALSE. number = 1 ELSE ! a common arc overlap = .TRUE. number = 2 END IF END IF ! separate arcs of same circle [no ELSE; we RETURN, with number = 0] END IF ! either circle is complete, or neither is ELSE IF (antipole_gap <= tolerance) THEN ! antipodal IF (dot_a /= -dot_b) RETURN ! no intersection ! From here on, assume dot_a == -dot_b: IF (ring_a) THEN overlap = .TRUE. number = 2 point1_uvec = first_b_uvec point2_uvec = last_b_uvec ELSE IF (ring_b) THEN overlap = .TRUE. number = 2 point1_uvec = first_a_uvec point2_uvec = last_a_uvec ELSE ! neither circle is complete ! Because of antipodal relationship, azimuths are confusing. ! Redefine azimuths of arc b in terms of pole a, ! while reversing the direction along arc b to make it ! counterclockwise about a: first_b_radians = -DRelative_Compass(pole_a_uvec, last_b_uvec) last_b_radians = -DRelative_Compass(pole_a_uvec, first_b_uvec) IF (last_b_radians < first_b_radians) last_b_radians = last_b_radians + Two_Pi IF ((first_b_radians >= first_a_radians).AND. & &(first_b_radians <= last_a_radians)) THEN ! 1st of b is in a. point1_uvec = first_b_uvec IF (last_b_radians <= last_a_radians) THEN point2_uvec = last_b_uvec ELSE point2_uvec = last_a_uvec END IF IF ((point1_uvec(1) == point2_uvec(1)).AND. & &(point1_uvec(2) == point2_uvec(2)).AND. & &(point1_uvec(3) == point2_uvec(3))) THEN ! only one common point overlap = .FALSE. number = 1 ELSE ! a common arc overlap = .TRUE. number = 2 END IF ELSE IF ((first_a_radians >= first_b_radians).AND. & &(first_a_radians <= last_b_radians)) THEN ! 1st of a is in b. point1_uvec = first_a_uvec IF (last_a_radians <= last_b_radians) THEN point2_uvec = last_a_uvec ELSE point2_uvec = last_b_uvec END IF IF ((point1_uvec(1) == point2_uvec(1)).AND. & &(point1_uvec(2) == point2_uvec(2)).AND. & &(point1_uvec(3) == point2_uvec(3))) THEN ! only one common point overlap = .FALSE. number = 1 ELSE ! a common arc overlap = .TRUE. number = 2 END IF END IF ! separate arcs on same circle [no ELSE; we RETURN, with number = 0] END IF ! either circle is complete, or neither is ELSE ! **** normal case; unrelated poles ***** ! Each small circle lies in its own plane. ! These two planes intersect in a line in 3-D. ! Find "offset" = (non-unit) vector to point on this line ! which is closest to origin: a_dot_b = DDot(pole_a_uvec, pole_b_uvec) determinant = 1.0D0 - a_dot_b**2 IF (determinant == 0.0D0) THEN WRITE (*,"(' ERROR: Determinant = 0.0D0 in DCircles_Intersect')") CALL DTraceback END IF t = 1.0D0 / determinant inverse(1,1) = t ! t * matrix(2,2) inverse(1,2) = -t * a_dot_b ! -t * matrix(1,2) inverse(2,1) = -t * a_dot_b ! -t * matrix(2,1) inverse(2,2) = t ! t * matrix(1,1) alpha = inverse(1,1)*dot_a + inverse(1,2)*dot_b beta = inverse(2,1)*dot_a + inverse(2,2)*dot_b offset(1) = alpha*pole_a_uvec(1) + beta*pole_b_uvec(1) offset(2) = alpha*pole_a_uvec(2) + beta*pole_b_uvec(2) offset(3) = alpha*pole_a_uvec(3) + beta*pole_b_uvec(3) min_radius = DLength(offset) IF (min_radius == 1.0D0) THEN ! circles osculate at one point ! Check longitudes about poles to see if point is in arcs: IF (.NOT.ring_a) THEN point_wrt_a_radians = -DRelative_Compass(pole_a_uvec, offset) IF (point_wrt_a_radians < first_a_radians) point_wrt_a_radians = point_wrt_a_radians + Two_Pi END IF IF (ring_a .OR. & &((point_wrt_a_radians >= first_a_radians).AND. & & (point_wrt_a_radians <= last_a_radians))) THEN IF (.NOT.ring_b) THEN point_wrt_b_radians = -DRelative_Compass(pole_b_uvec, offset) IF (point_wrt_b_radians < first_b_radians) point_wrt_b_radians = point_wrt_b_radians + Two_Pi END IF IF (ring_b .OR. & &((point_wrt_b_radians >= first_b_radians).AND. & & (point_wrt_b_radians <= last_b_radians))) THEN number = 1 point1_uvec = offset END IF ! in arc b [no ELSE; we RETURN with number = 0] END IF ! in arc a [no ELSE; we RETURN with number = 0] ELSE IF (min_radius < 1.) THEN ! two intersection points between circles ! Find vector parallel to line of intersection of the two ! planes which contain the two small circles: CALL DCross (pole_a_uvec, pole_b_uvec, t_vec) CALL DMake_Uvec (t_vec, line_uvec) ! Now, points on this line are expressed by ! vector "offset" plus any multiple of vector "line_uvec". ! Call the multiple "along". ! Solve the hyperbolic equation that says that the ! radius of a point on this line is 1.00 ! (i.e., it is also a point on the unit sphere): ! min_radius**2 + along**2 = 1.00**2 along = DSQRT ( 1.0D0 - min_radius**2 ) point1_uvec = offset + along * line_uvec point2_uvec = offset - along * line_uvec ! Remember, these are only tentative solutions. ! Check longitudes about poles to see if point1 is in arcs: IF (.NOT.ring_a) THEN point_wrt_a_radians = -DRelative_Compass(pole_a_uvec, point1_uvec) IF (point_wrt_a_radians < first_a_radians) point_wrt_a_radians = point_wrt_a_radians + Two_Pi END IF IF (ring_a .OR. & &((point_wrt_a_radians >= first_a_radians).AND. & & (point_wrt_a_radians <= last_a_radians))) THEN IF (.NOT.ring_b) THEN point_wrt_b_radians = -DRelative_Compass(pole_b_uvec, point1_uvec) IF (point_wrt_b_radians < first_b_radians) point_wrt_b_radians = point_wrt_b_radians + Two_Pi END IF IF (ring_b .OR. & &((point_wrt_b_radians >= first_b_radians).AND. & & (point_wrt_b_radians <= last_b_radians))) THEN number = 1 END IF ! point1 also in arc b: a solution. END IF ! point1 is in arc a ! Check longitudes about poles to see if point2 is in arcs: IF (.NOT.ring_a) THEN point_wrt_a_radians = -DRelative_Compass(pole_a_uvec, point2_uvec) IF (point_wrt_a_radians < first_a_radians) point_wrt_a_radians = point_wrt_a_radians + Two_Pi END IF IF (ring_a .OR. & &((point_wrt_a_radians >= first_a_radians).AND. & & (point_wrt_a_radians <= last_a_radians))) THEN IF (.NOT.ring_b) THEN point_wrt_b_radians = -DRelative_Compass(pole_b_uvec, point2_uvec) IF (point_wrt_b_radians < first_b_radians) point_wrt_b_radians = point_wrt_b_radians + Two_Pi END IF IF (ring_b .OR. & &((point_wrt_b_radians >= first_b_radians).AND. & & (point_wrt_b_radians <= last_b_radians))) THEN IF (number == 1) THEN ! add a second solution number = 2 ELSE ! point2 is the first and only solution number = 1 point1_uvec = point2_uvec END IF ! is this the first or the second solution? END IF ! point2 in arc b: a solution. END IF ! point2 is in arc a END IF ! min_radius = 1.0D0, or < 1.0D0 [no ELSE; we RETURN with number = 0] END IF ! poles parallel, or antipodal, or unrelated END SUBROUTINE DCircles_Intersect LOGICAL FUNCTION DClockways (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*8 :: azim_before, azim_end, azim_start, internal_bend, & & last_lat, last_lon, pazim0, pazim1, size, start_bend, total_bend REAL*8, 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.D0 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 = DArc(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.D0 ! one counterclockwise (negative; left) circuit azim_start = Pi_over_2 + DRelative_Compass (t0_uvec, pole_uvec) azim_end = azim_start ! same trend as start ELSE ! incomplete circle pazim0 = DRelative_Compass (pole_uvec, t0_uvec) pazim1 = DRelative_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 + DRelative_Compass (t0_uvec, pole_uvec) azim_end = Pi_over_2 + DRelative_Compass (t1_uvec, pole_uvec) END IF internal_bend = (pazim1 - pazim0) * DCOS(size) ELSE ! great circle arc; COS(size) = 0.D0 internal_bend = 0.D0 t0_uvec = last_uvec t1_uvec = ai_pathlib(1:3,i,path) CALL DCross (t0_uvec, t1_uvec, t_vec) CALL DMake_Uvec (t_vec, pole_uvec) azim_start = Pi_over_2 + DRelative_Compass (t0_uvec, pole_uvec) azim_end = Pi_over_2 + DRelative_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.D0) 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 DClockways = (total_bend > 0.D0) ELSE WRITE (*,"(' ERROR: Single-point paths cannot have the DClockways property.')") CALL DTraceback END IF ELSE ! wrong path level WRITE (*,"(' ERROR: Wrong path level for argument of DClockways: ',I2)") ai_path_level(path) CALL DTraceback END IF END FUNCTION DClockways REAL*8 FUNCTION DCompass (from_uvec, to_uvec) ! Returns azimuth (in radians, clockwise from North) ! of the great-circle arc from "from_uvec" to "to_uvec", ! measured at location "from_uvec". ! Does NOT work at North or South pole! (See Relative_Compass.) IMPLICIT NONE REAL*8, DIMENSION(3), INTENT(IN) :: from_uvec, to_uvec REAL*8 :: ve, vn REAL*8, DIMENSION(3) :: omega_vec, omega_uvec, & & Phi, Theta, & & v_vec IF ((from_uvec(1) == 0.D0).AND.(from_uvec(2) == 0.D0)) THEN WRITE (*,"(' ERROR: DCompass undefined at N or S pole. Use DRelative_Compass.')") CALL DTraceback ELSE IF ((from_uvec(1) == to_uvec(1)).AND. & &(from_uvec(2) == to_uvec(2)).AND. & &(from_uvec(3) == to_uvec(3))) THEN WRITE (*,"(' ERROR: DCompass bearing from point to itself undefined.')") CALL DTraceback ELSE IF ((from_uvec(1) == -to_uvec(1)).AND. & &(from_uvec(2) == -to_uvec(2)).AND. & &(from_uvec(3) == -to_uvec(3))) THEN WRITE (*,"(' ERROR: DCompass bearing from point to antipode undefined.')") CALL DTraceback ELSE ! Normal case: CALL DCross (from_uvec, to_uvec, omega_vec) CALL DMake_Uvec(omega_vec, omega_uvec) CALL DCross (omega_uvec, from_uvec, v_vec) CALL DLocal_Theta(from_uvec, Theta) CALL DLocal_Phi (from_uvec, Phi) vn = -DDot(v_vec, Theta) ve = DDot(v_vec, Phi) DCompass = DATAN2(ve, vn) END IF END FUNCTION DCompass REAL*8 FUNCTION DConformal_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*8, DIMENSION(3), INTENT(IN) :: uvec REAL*8 :: colatitude, latitude SELECT CASE(mp_projection_number) CASE (1) ! Mercator colatitude = DArc (uvec, mp_pole_3_uvec) DConformal_Deflation = DSIN(colatitude) CASE (2) ! Lambert Conformal Conic ! mp_proj2_F = DCOS(latitude_1_radians) * DTAN(Pi_over_4 + latitude_1_radians/2.D0)**mp_proj2_n / mp_proj2_n latitude = Pi_over_2 - DArc(uvec, mp_pole_3_uvec) DConformal_Deflation = DCOS(latitude)*DTAN(0.5D0*(Pi_over_2 + latitude))**mp_proj2_n / & & (mp_proj2_F * mp_proj2_n) CASE (6) ! Stereographic colatitude = DArc (uvec, mp_projpoint_uvec) DConformal_Deflation = (1.D0 + DCOS(colatitude)) * 0.5D0 CASE DEFAULT DConformal_Deflation = 1.00D0 END SELECT END FUNCTION DConformal_Deflation SUBROUTINE DCorners (projpoint_uvec) ! After DSet_[projection] has defined the (up to) 4 walls of ! the projectable "rind" on the sphere, and before DTrace_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*8, 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*8, DIMENSION(3), INTENT(IN) :: projpoint_uvec INTEGER :: i, number LOGICAL :: overlap, using_wall REAL*8, PARAMETER :: tolerance = 0.00001D0 REAL*8 :: close1, close2, dot1a, dot1b, dot2a, dot2b, error, error1, & & error2, maxdot, sidestep REAL*8, DIMENSION(3), PARAMETER :: north_uvec = (/ 1.D0, 0.D0, 0.D0 /) REAL*8, DIMENSION(3), PARAMETER :: zero_vec = (/ 0.D0, 0.D0, 0.D0 /) ! (This is a useful flag in "CALL DCircles_Intersect" that we ! want to consider the whole circle, not an arc.) REAL*8, 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 DNorthEast_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.D0 <= maxdot < 1.D0 IF ((maxdot < 0.0D0).OR.(maxdot >= 1.0D0)) THEN WRITE (*,"(' ERROR: Violation of limits 0.0D0 <= maxdot < 1.0D0 on rind side ',I1)") i CALL DTraceback END IF ! check that this wall does not exclude projection point! IF (DDot(pole_applies_uvec,projpoint_uvec) > 0.D0) THEN IF (DDot(projpoint_uvec,pole_uvec) > maxdot) THEN WRITE (*,"(' ERROR: Projection point excluded by rind wall ',I1)") i CALL DTraceback END IF END IF ! find "any old" uvec in the wall for a default first/last: sidestep = DSQRT(1.D0 - maxdot**2) CALL DCross (projpoint_uvec, pole_uvec, t_vec) IF (DLength(t_vec) == 0.D0) CALL DCross (north_uvec, pole_uvec, t_vec) IF (DLength(t_vec) == 0.D0) CALL DCross (right_uvec, pole_uvec, t_vec) CALL DMake_Uvec(t_vec, perp1_uvec) CALL DCross (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 DCircles_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 DTraceback END IF IF (number > 0) THEN ! found some possible points ! choose which answer is most reasonable (least error) dot1a = DDot(point1_uvec, mp_top_applies_uvec) dot1b = DDot(point1_uvec, mp_left_applies_uvec) error1 = -MIN(dot1a, dot1b) ! problem if either is <0.D0 error = error1 trial_uvec = point1_uvec IF (number == 2) THEN dot2a = DDot(point2_uvec, mp_top_applies_uvec) dot2b = DDot(point2_uvec, mp_left_applies_uvec) error2 = -MIN(dot2a, dot2b) ! problem if either is <0.D0 IF (DABS(error2 - error1) <= tolerance) THEN ! essentially equal ! choose the one closer to projpoint: close1 = DDot(point1_uvec, projpoint_uvec) close2 = DDot(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 DCircles_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 DTraceback END IF IF (number > 0) THEN ! found some possible points ! choose which answer is most reasonable (least error) dot1a = DDot(point1_uvec, mp_top_applies_uvec) dot1b = DDot(point1_uvec, mp_right_applies_uvec) error1 = -MIN(dot1a, dot1b) ! problem if either is <0.D0 error = error1 trial_uvec = point1_uvec IF (number == 2) THEN dot2a = DDot(point2_uvec, mp_top_applies_uvec) dot2b = DDot(point2_uvec, mp_right_applies_uvec) error2 = -MIN(dot2a, dot2b) ! problem if either is <0.D0 IF (DABS(error2 - error1) <= tolerance) THEN ! essentially equal ! choose the one closer to projpoint: close1 = DDot(point1_uvec, projpoint_uvec) close2 = DDot(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 DCircles_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 DTraceback END IF IF (number == 2) THEN ! found two points ! decide which is right, which is left IF (DDot(point1_uvec,right_uvec) > DDot(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 DCircles_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 DTraceback END IF IF (number > 0) THEN ! found some possible points ! choose which answer is most reasonable (least error) dot1a = DDot(point1_uvec, mp_bottom_applies_uvec) dot1b = DDot(point1_uvec, mp_left_applies_uvec) error1 = -MIN(dot1a, dot1b) ! problem if either is <0.D0 error = error1 trial_uvec = point1_uvec IF (number == 2) THEN dot2a = DDot(point2_uvec, mp_bottom_applies_uvec) dot2b = DDot(point2_uvec, mp_left_applies_uvec) error2 = -MIN(dot2a, dot2b) ! problem if either is <0.D0 IF (DABS(error2 - error1) <= tolerance) THEN ! essentially equal ! choose the one closer to projpoint: close1 = DDot(point1_uvec, projpoint_uvec) close2 = DDot(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 DCircles_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 DTraceback END IF IF (number > 0) THEN ! found some possible points ! choose which answer is most reasonable (least error) dot1a = DDot(point1_uvec, mp_bottom_applies_uvec) dot1b = DDot(point1_uvec, mp_right_applies_uvec) error1 = -MIN(dot1a, dot1b) ! problem if either is <0.D0 error = error1 trial_uvec = point1_uvec IF (number == 2) THEN dot2a = DDot(point2_uvec, mp_bottom_applies_uvec) dot2b = DDot(point2_uvec, mp_right_applies_uvec) error2 = -MIN(dot2a, dot2b) ! problem if either is <0.D0 IF (DABS(error2 - error1) <= tolerance) THEN ! essentially equal ! choose the one closer to projpoint: close1 = DDot(point1_uvec, projpoint_uvec) close2 = DDot(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 DCircles_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 DTraceback END IF IF (number == 2) THEN ! found two points ! decide which is zenith (up), which nadir (down) IF (DDot(point1_uvec,north_uvec) > DDot(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 DCorners SUBROUTINE DCross (a_vec, b_vec, c_vec) ! vector cross product: a x b = c IMPLICIT NONE REAL*8, DIMENSION(3), INTENT(IN) :: a_vec, b_vec REAL*8, DIMENSION(3), INTENT(OUT) :: c_vec c_vec(1) = a_vec(2)*b_vec(3) - a_vec(3)*b_vec(2) c_vec(2) = a_vec(3)*b_vec(1) - a_vec(1)*b_vec(3) c_vec(3) = a_vec(1)*b_vec(2) - a_vec(2)*b_vec(1) END SUBROUTINE DCross REAL*8 FUNCTION DDot (a_vec, b_vec) ! returns scalar (dot) product of two 3-component vectors IMPLICIT NONE REAL*8, DIMENSION(3), INTENT(IN) :: a_vec, b_vec DDot = a_vec(1)*b_vec(1) + a_vec(2)*b_vec(2) + a_vec(3)*b_vec(3) END FUNCTION DDot REAL*8 FUNCTION DEasting(delta_lon_degrees) !returns positive result, 0.D0 to 359.999...D0 IMPLICIT NONE REAL*8, INTENT(IN) :: delta_lon_degrees DEasting = MOD((delta_lon_degrees + 720.D0),360.D0) END FUNCTION DEasting INTEGER FUNCTION DIAbove(x) ! returns first integer >= x, unlike INT() or NINT(): IMPLICIT NONE REAL*8, INTENT(IN) :: x INTEGER answer REAL*8 :: t answer = INT(x) IF (x > 0.D0) THEN t = 1.00D0 * answer IF (x > t) THEN answer = answer + 1 END IF END IF DIAbove = answer END FUNCTION DIAbove CHARACTER*1 FUNCTION DIn_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. ! Input "uvec" is the input unit-vector representing a position ! on the surface of a unit sphere representing the planet. ! Returns 'I' for inside, 'O' for outside, 'B' for on-boundary. IMPLICIT NONE REAL*8, DIMENSION(3), INTENT(IN) :: uvec INTEGER :: i REAL*8 :: upper_limit, value DIn_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.(DDot(uvec, mp_top_applies_uvec) > 0.D0)) THEN value = DDot (uvec, mp_top_uvec) upper_limit = mp_top_maxdot ELSE ! top limit doesn't apply; give a "free pass" value = 0.D0 upper_limit = 1.D0 END IF CASE(2) ! left boundary IF (mp_using_left.AND.(DDot(uvec, mp_left_applies_uvec) > 0.D0)) THEN value = DDot (uvec, mp_left_uvec) upper_limit = mp_left_maxdot ELSE ! left limit doesn't apply; give a "free pass" value = 0.D0 upper_limit = 1.D0 END IF CASE(3) ! right boundary IF (mp_using_right.AND.(DDot(uvec, mp_right_applies_uvec) > 0.D0)) THEN value = DDot (uvec, mp_right_uvec) upper_limit = mp_right_maxdot ELSE ! right limit doesn't apply; give a "free pass" value = 0.D0 upper_limit = 1.D0 END IF CASE(4) ! bottom boundary IF (mp_using_bottom.AND.(DDot(uvec, mp_bottom_applies_uvec) > 0.D0)) THEN value = DDot (uvec, mp_bottom_uvec) upper_limit = mp_bottom_maxdot ELSE ! bottom limit doesn't apply; give a "free pass" value = 0.D0 upper_limit = 1.D0 END IF END SELECT IF (value > upper_limit) THEN DIn_Rind = 'O' ! outside RETURN ELSE IF (value == upper_limit) THEN DIn_Rind = 'B' ! ELSE ! value < upper_limit ! DIn_Rind remains as 'I' or 'B' (whichever it was). END IF END DO END FUNCTION DIn_Rind LOGICAL FUNCTION DL5_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*8, DIMENSION(3), INTENT(IN) :: uvec CHARACTER*1 :: c1 REAL*8 :: x_meters, x_points, y_meters, y_points !------------------------- c1 = DIn_Rind(uvec) IF (c1 == 'O') THEN ! outside rind DL5_In_Window = .FALSE. ELSE CALL DProject (uvec = uvec, x = x_meters, y = y_meters) CALL DMeters_2_Points (x_meters,y_meters, x_points,y_points) c1 = DIn_Window (x_points, y_points) IF (c1 == 'O') THEN ! outside window DL5_In_Window = .FALSE. ELSE DL5_In_Window = .TRUE. END IF END IF END FUNCTION DL5_In_Window REAL*8 FUNCTION DLength(a_vec) IMPLICIT NONE REAL*8, DIMENSION(3), INTENT(IN) :: a_vec DOUBLE PRECISION :: t t = a_vec(1)*a_vec(1) + & & a_vec(2)*a_vec(2) + & & a_vec(3)*a_vec(3) IF (t == 0.0D0) THEN DLength = 0.0D0 ELSE DLength = DSQRT(t) END IF END FUNCTION DLength SUBROUTINE DLocal_Phi (b_, Phi) ! returns local East-pointing unit vector in Cartesian coordinates ! for location b_; not intended to work at the poles! REAL*8, DIMENSION(3), INTENT(IN) :: b_ REAL*8, DIMENSION(3), INTENT(OUT) :: Phi REAL*8, DIMENSION(3) :: temp IF (b_(1) == 0.D0) THEN IF (b_(2) == 0.D0) THEN WRITE (*,"(' ERROR: DLocal_Phi was requested for N or S pole.')") CALL DTraceback END IF END IF temp(1) = - b_(2) temp(2) = b_(1) temp(3) = 0.D0 CALL DMake_Uvec(temp, Phi) END SUBROUTINE DLocal_Phi SUBROUTINE DLocal_Theta (b_, Theta) ! returns local South-pointing unit vector in Cartesian coordinates ! for location b_; not intended to work at the poles! REAL*8, DIMENSION(3), INTENT(IN) :: b_ REAL*8, DIMENSION(3), INTENT(OUT) :: Theta REAL*8, DIMENSION(3) :: temp REAL*8 :: equat, new_equat equat = DSQRT(b_(1)**2 + b_(2)**2) !equatorial component IF (equat == 0.D0) THEN WRITE (*,"(' ERROR: DLocal_Theta was requested for N or S pole.')") CALL DTraceback END IF new_equat = b_(3) ! swap components in a meridional plane temp(3) = - equat ! " temp(1) = new_equat * b_(1) / equat ! partition new equatorial component temp(2) = new_equat * b_(2) / equat ! " CALL DMake_Uvec (temp, Theta) END SUBROUTINE DLocal_Theta SUBROUTINE DLonLat_2_ThetaPhi (lon, lat, theta, phi) ! "lon" is East longitude in degrees. ! "lat" is North latitude in degrees. ! "theta" is co-latitude, from N pole, in radians ! "phi" is East longitude in radians IMPLICIT NONE REAL*8, INTENT(IN) :: lon, lat REAL*8, INTENT(OUT) :: theta, phi IF ((lat > 90.D0).OR.(lat < -90.D0)) THEN WRITE (*,"(' ERROR: Latitude outside legal range.')") CALL DTraceback ELSE theta = radians_per_degree * (90.0D0 - lat) phi = radians_per_degree * lon END IF END SUBROUTINE DLonLat_2_ThetaPhi SUBROUTINE DLonLat_2_Uvec (lon, lat, uvec) IMPLICIT NONE REAL*8, INTENT(IN) :: lon, lat REAL*8, DIMENSION(3), INTENT(OUT) :: uvec REAL*8 :: theta, phi CALL DLonLat_2_ThetaPhi (lon, lat, theta, phi) CALL DThetaPhi_2_Uvec (theta, phi, uvec) END SUBROUTINE DLonLat_2_Uvec REAL*8 FUNCTION DMagnitude (b_) REAL*8, DIMENSION(3), INTENT(IN) :: b_ DMagnitude = DSQRT(b_(1)**2 +b_(2)**2 + b_(3)**2) END FUNCTION DMagnitude ! SUBROUTINE MATHERRQQ( name, length, info, retcode) ! !Provided so that domain errors, underflows, etc. ! ! can be trapped and debugged; otherwise program crashes! ! USE DFLIB ! INTEGER(2) length, retcode ! CHARACTER(length) name ! RECORD /MTH$E_INFO/ info ! PRINT *, "Entered MATHERRQQ" ! PRINT *, "Failing function is: ", name ! PRINT *, "Error type is: ", info.errcode ! IF ((info.ftype == TY$REAL4 ).OR.(info.ftype == TY$REAL8)) THEN ! PRINT *, "Type: REAL" ! PRINT *, "Enter the desired function result: " ! READ(*,*) info.r8res ! retcode = 1 ! ELSE IF ((info.ftype == TY$CMPLX8 ).OR.(info.ftype == TY$CMPLX16)) THEN ! PRINT *, "Type: COMPLEX" ! PRINT *, "Enter the desired function result: " ! READ(*,*) info.c16res ! retcode = 1 ! END IF ! END SUBROUTINE MATHERRQQ SUBROUTINE DMake_Uvec (vector, uvec) ! Shortens or lengthens a three-component vector to make a unit vector. IMPLICIT NONE INTEGER :: i DOUBLE PRECISION, DIMENSION(3), INTENT(IN) :: vector DOUBLE PRECISION, DIMENSION(3), INTENT(OUT):: uvec DOUBLE PRECISION :: factor, size size = DLength(vector) IF (size > 0.0D0) THEN factor = 1.0D0 / size uvec = vector * factor ELSE WRITE (*,"(' ERROR: Cannot DMake_Uvec of (0.0D0, 0.0D0, 0.0D0).')") CALL DTraceback END IF END SUBROUTINE DMake_Uvec SUBROUTINE DMeters_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*8, INTENT(IN) :: x_meters, y_meters REAL*8, 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 DMeters_2_Points SUBROUTINE DNorthEast_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 DRelative_Compass ! and SUBROUTINE DTurn_To call this routine, they can work ! together even when location_uvec is at, or near, a pole. IMPLICIT NONE REAL*8, DIMENSION(3), INTENT(IN) :: location_uvec REAL*8, DIMENSION(3), INTENT(OUT) :: north_uvec, east_uvec REAL*8, PARAMETER :: pole_area = 7.615D-7 ! (0.05 degrees -> radians -> squared) REAL*8 :: equat2 REAL*8, DIMENSION(3) :: t_vec equat2 = location_uvec(1)**2 + location_uvec(2)**2 IF (equat2 > pole_area) THEN ! normal case: t_vec(1) = -location_uvec(2) t_vec(2) = +location_uvec(1) t_vec(3) = 0.0D0 CALL DMake_Uvec (t_vec, east_uvec) ELSE ! very close to N or S pole; act as if on 0E meridian: east_uvec = (/ 0.D0, 1.D0, 0.D0 /) END IF CALL DCross (location_uvec, east_uvec, north_uvec) END SUBROUTINE DNorthEast_Convention SUBROUTINE DOld_Complex_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 DProcess_L4_Paths for ! subdivision into small pieces, conversion to Bezier curves, ! and projection onto the projection plane. IMPLICIT NONE CHARACTER*1 :: abyte, byte, goal_memo CHARACTER*1, DIMENSION(1:ai_longest) :: seg_memo ! O, B, I CHARACTER*1, DIMENSION(0:ai_longest) :: point_memo ! O, B, I !{N.B. Following DOUBLE PRECISION was in the original MODULE Map_Projections.} 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, n_area, new_path, & & new_segments, number, path, segments, 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, trace_clockwise, unchanged LOGICAL, DIMENSION(ai_longest) :: new_segment_used REAL*8 :: azimuth, base_of_t, clockwise_distance, counterclockwise_distance, end_lat, end_lon, far_radians, from_t_wall, phi, & & phi_end, phi_start, phi0, phi1, start_lat, start_lon, t, t_wall, to_t_wall, zeta, zeta0, zeta1 REAL*8, 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*8, DIMENSION(16) :: t_edge, t15_edge, a_edge, b_edge, g_edge REAL*8, DIMENSION(inout_list_limit) :: t_1_to_5, & & a_inout_list, b_inout_list, g_inout_list REAL*8, DIMENSION(1:ai_longest) :: t15_memo DO WHILE (ai_Ln_paths(5) > 0) path = DNext_Path(level = 5) filled = ai_filled(path) segments = ai_segments(path) ! Shortcut complex logic (below) if path has only an initial point, ! but no segments. Pass it through (even though it will have no ! graphical effect) if it falls within the rind. IF (segments == 0) THEN ! path has only an initial point start_uvec(1:3) = ai_pathlib(1:3, 0, path) byte = DIn_Rind(start_uvec) IF ((byte == 'I').OR.(byte == 'B')) THEN !Transform path to level 4, and call for Processing: ai_path_level(path) = 4 ai_Ln_paths(5) = ai_Ln_paths(5) - 1 ai_Ln_paths(4) = ai_Ln_paths(4) + 1 CALL DProcess_L4_Paths() ELSE ! just eliminate this path: ai_path_level(path) = 0 ai_Ln_paths(5) = ai_Ln_paths(5) - 1 ai_total_paths = ai_total_paths - 1 END IF ! inside, or outside CYCLE ! so as to consider the next L5 path in library END IF ! path had no segments, but just an initial point ! 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 = DNext_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) = DIn_Rind (last_uvec) unchanged = (point_memo(0) == 'I').OR.(point_memo(0) == 'B') new_segments = 0 DO i = 1, ai_segments(path) !characterize segment by pole, far_radians (from pole), complete_circle (T/F), and null_arc (T/F): IF (ai_bent(i,path)) THEN goal_uvec = ai_pathlib(4:6,i,path) pole_uvec = ai_pathlib(1:3,i,path) far_radians = DArc (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 DCross (last_uvec, goal_uvec, t_vec) null_arc = (MAX(DABS(t_vec(1)),DABS(t_vec(2)),DABS(t_vec(3))) == 0.0D0) IF (null_arc) THEN ! pole is undefined IF (DABS(goal_uvec(3)) > 0.999D0) THEN t_uvec = (/ 1.0D0, 0.0D0, 0.0D0 /) ELSE t_uvec = (/ 0.0D0, 0.0D0, 1.0D0 /) END IF CALL DCross (t_uvec, goal_uvec, t_vec) !result will be perpendicular to goal_uvec END IF CALL DMake_Uvec (t_vec, pole_uvec) far_radians = Pi_over_2 complete_circle = .FALSE. END IF abyte = DIn_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 = DCOS(DArc(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 = -DRelative_Compass (pole_uvec, last_uvec) phi1 = -DRelative_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.D0 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.D0 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.D0 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.D0 END SELECT CALL DCircles_Intersect & &(pole_a_uvec, dot_a, first_a_uvec, last_a_uvec, & & pole_b_uvec, dot_b, first_b_uvec, last_b_uvec, & ! input & overlap, number, point1_uvec, point2_uvec) ! output IF (number >= 1) THEN unchanged = .FALSE. intersections = intersections + 1 phi = -DRelative_Compass (pole_uvec, point1_uvec) IF (phi < phi0) phi = phi + Two_Pi t = (phi - phi0) / (phi1 - phi0) ! 0.D0 to 1.D0 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 = DRelative_Compass (pole_b_uvec, last_b_uvec) zeta1 = DRelative_Compass (pole_b_uvec, first_b_uvec) IF (zeta1 <= zeta0) zeta1 = zeta1 + Two_Pi zeta = DRelative_Compass (pole_b_uvec, point1_uvec) IF (zeta < zeta0) zeta = zeta + Two_Pi t_wall = (zeta - zeta0) / (zeta1 - zeta0) t_wall = MAX(0.D0, MIN(t_wall, 0.99999D0)) 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 = -DRelative_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 = DRelative_Compass (pole_b_uvec, point2_uvec) IF (zeta < zeta0) zeta = zeta + Two_Pi t_wall = (zeta - zeta0) / (zeta1 - zeta0) t_wall = MAX(0.D0, MIN(t_wall, 0.99999D0)) 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 DSort_Lists(intersections, t_edge, t15_edge, a_edge, b_edge, g_edge) IF (intersections > 0) THEN ! {N.B. See ELSE paragraph far below for case of no intersections.} 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 DTraceback 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.5D0*(phi_start + phi_end) CALL DTurn_To (azimuth, pole_uvec, far_radians, & ! inputs & omega_uvec, t_uvec) seg_memo(new_segments) = DIn_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.0D0) 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 DTraceback 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.5D0*(phi_start + phi1) CALL DTurn_To (azimuth, pole_uvec, far_radians, & ! inputs & omega_uvec, t_uvec) seg_memo(new_segments) = DIn_Rind (t_uvec) point_memo(new_segments) = abyte END IF ! divided segment needs completion ELSE ! no intersections found for this segment !This could be a problem, IFF we went from 'I' to 'O' or from 'O' to 'I': IF (((point_memo(i-1) == 'I').AND.(abyte == 'O')) .OR. & & ((point_memo(i-1) == 'O').AND.(abyte == 'I'))) THEN CALL DUvec_2_LonLat(last_uvec, start_lon, start_lat) CALL DUvec_2_LonLat(goal_uvec, end_lon, end_lat) WRITE (*, "(' ERROR! Segment start point (',F10.5,', ',F9.5,') is ',A1 & & /' and end point (',F10.5,', ',F9.5,') is ',A1)") & & start_lon, start_lat, point_memo(i-1), end_lon, end_lat, abyte WRITE (*, "(' but no intersections with any rind walls were detected.')") CALL DTraceback END IF ! topological impossibility was detected new_segments = new_segments + 1 IF (new_segments >= ai_longest) THEN WRITE (*,"(' ERROR: Path length > ai_longest =',I6)") ai_longest CALL DTraceback 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 seg_memo(new_segments) = abyte ! New in DMap_Projections; old method of computing and classifying midpoint was buggy! 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 DProcess_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 = DClockways(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 DTraceback 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 DSort_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 DNew_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 DSmall_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 DGreat_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 DTraceback END IF !Guard against tracing-rind from/to undefined points: IF ((from_t_wall == 0.0D0).OR.(to_t_wall == 0.0D0)) THEN WRITE (*, "(' ERROR! Undefined start/end point in CALL DTrace_Rind.')") !WRITE (*, "(' mp_rind_tracings = ', I10)") mp_rind_tracings !WRITE (*, "(' mp_ProcL5_entry = ', I10)") mp_ProcL5_entry CALL DTraceback END IF !Be sure to trace around by the shorter way! clockwise_distance = from_t_wall - to_t_wall IF (clockwise_distance < 0.0D0) clockwise_distance = clockwise_distance + mp_walls_in_use counterclockwise_distance = to_t_wall - from_t_wall IF (counterclockwise_distance < 0.0D0) counterclockwise_distance = counterclockwise_distance + mp_walls_in_use trace_clockwise = (clockwise_distance <= counterclockwise_distance) CALL DTrace_Rind(from_t_wall, trace_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 DEnd_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 DTraceback 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 DNew_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 DSmall_To_L45 (pole_uvec, t_uvec) ELSE ! great circle arc t_uvec = ai_pathlib(1:3,i,new_path) CALL DGreat_To_L45 (t_uvec) END IF ELSE ! not in path CALL DNew_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 DSmall_To_L45 (pole_uvec, t_uvec) ELSE t_uvec = ai_pathlib(1:3,i,new_path) CALL DGreat_To_L45 (t_uvec) END IF END IF ! in path or not ELSE ! curve doesn't show IF (ai_in_path) THEN CALL DEnd_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 DEnd_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 DProcess_L4_Paths END DO ! on all L5 paths presented END SUBROUTINE DOld_Complex_Process_L5_Paths SUBROUTINE DPoints_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*8, INTENT(IN) :: x_points, y_points LOGICAL, INTENT(OUT) :: success REAL*8, INTENT(OUT) :: lon, lat REAL*8 :: x_meters, y_meters REAL*8, DIMENSION(3) :: uvec CALL DPoints_2_Meters (x_points,y_points, x_meters,y_meters) CALL DReject (x_meters,y_meters, success, uvec) IF (success) CALL DUvec_2_LonLat (uvec, lon, lat) END SUBROUTINE DPoints_2_LonLat SUBROUTINE DPoints_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*8, INTENT(IN) :: x_points, y_points REAL*8, 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 DPoints_2_Meters SUBROUTINE DProcess_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 DAdobe_Illustrator, ! the test was applied in DCurve_To_L12, and these segments ! do not pass through that user interface.) IMPLICIT NONE INTEGER :: i, path LOGICAL inline1, inline2 REAL*8 :: 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 = DNext_Path(level = 3) CALL DMeters_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 DMeters_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 DMeters_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 DMeters_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 = DSQRT(ux**2 + uy**2) ! > 0.D0 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 = DSQRT(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 = DSQRT(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 DProcess_L2_Paths END DO END SUBROUTINE DProcess_L3_Paths SUBROUTINE DProcess_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*8, INTENT(IN) :: x_meters, y_meters, page_angle_radians, & & lr_fraction, ud_fraction CHARACTER*(*),INTENT(IN) :: text REAL*8 :: x_points, y_points CALL DMeters_2_Points (x_meters,y_meters, x_points,y_points) CALL DProcess_L2_Text (x_points, y_points, page_angle_radians, & & font_points, lr_fraction, ud_fraction, & & text) END SUBROUTINE DProcess_L3_Text SUBROUTINE DProcess_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 DProcess_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 LOGICAL :: null_great_circle REAL*8, PARAMETER :: piece_size = 0.174533D0 ! 10 degrees, in radians REAL*8 :: 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*8, 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 = DNext_Path(level = 4) last_uvec = ai_pathlib(1:3,0,path) ! initialize memory CALL DProject(uvec = last_uvec, x = x0_meters, y = y0_meters) CALL DNew_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 = DArc(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 = DRelative_Compass(ccw_pole_uvec, last_uvec) azim_b = azim_a - Two_Pi ELSE ! less than a complete circle azim_a = DRelative_Compass(ccw_pole_uvec, last_uvec) azim_b = DRelative_Compass(ccw_pole_uvec, goal_uvec) IF (azim_b > azim_a) azim_b = azim_b - Two_Pi ! ccw radians = azim_a - azim_b ! > 0.D0 curviness = DABS( DCOS(small_size_radians) ) pieces = MAX (1, DIAbove(curviness*radians/Pi_over_2) ) END IF ! complete small circle ! Also limit to no more than piece_size arclength. arc_length = radians * DSIN(small_size_radians) pieces = MAX(pieces, DIAbove(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 DCross(last_uvec, goal_uvec, t_vec) crossed = DLength(t_vec) ! > 0.D0? IF (crossed > 0.0D0) THEN ! normal case: null_great_circle = .FALSE. CALL DMake_Uvec(t_vec, ccw_pole_uvec) dotted = DDot(last_uvec, goal_uvec) radians = DATAN2(crossed, dotted) ! > 0.D0 azim_a= DRelative_Compass (ccw_pole_uvec, last_uvec) azim_b = azim_a - radians arc_length = radians pieces = DIAbove(radians/piece_size) ELSE ! rare case. ODDLY, this can happen even though DGreat_to_L45 guards against it! null_great_circle = .TRUE. radians = 0.0D0 arc_length = 0.0D0 azim_a = 0.0D0 azim_b = azim_a pieces = 0 END IF 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 DTurn_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 DProject(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.3516D0 * arc_length / pieces ! Bezier handle on initial end of piece: azimuth = Pi_over_2 + DRelative_Compass (last_piece_uvec, ccw_pole_uvec) CALL DTurn_To (azimuth, last_piece_uvec, & & forward, & ! inputs & omega_uvec, result_uvec) ! outputs CALL DProject (uvec = result_uvec, guide = guide_uvec, & & x = x1_meters, y = y1_meters) ! Bezier handle on final end of piece: azimuth = DRelative_Compass (guide_uvec, ccw_pole_uvec) - Pi_over_2 CALL DTurn_To (azimuth, guide_uvec, & & forward, & ! inputs & omega_uvec, result_uvec) ! outputs CALL DProject (uvec = result_uvec, guide = guide_uvec, & & x = x2_meters, y = y2_meters) CALL DCurve_To_L3 (x1_meters,y1_meters, & & x2_meters,y2_meters, & & x3_meters,y3_meters) ELSE ! piece is a great circle arc IF (.NOT.null_great_circle) THEN ! Bezier handle on initial end of piece: azim_1 = 0.6484D0*azim_0 + 0.3516D0*azim_3 CALL DTurn_To (azim_1, ccw_pole_uvec, & & small_size_radians, & ! inputs & omega_uvec, result_uvec) ! results CALL DProject (uvec = result_uvec, guide = guide_uvec, & & x = x1_meters, y = y1_meters) ! Bezier handle on final end of piece: azim_2 = 0.3516D0*azim_0 + 0.6484D0*azim_3 CALL DTurn_To (azim_2, ccw_pole_uvec, & & small_size_radians, & ! inputs & omega_uvec, result_uvec) ! results CALL DProject (uvec = result_uvec, guide = guide_uvec, & & x = x2_meters, y = y2_meters) CALL DCurve_To_L3 (x1_meters,y1_meters, & & x2_meters,y2_meters, & & x3_meters,y3_meters) END IF END IF 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 DEnd_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 DProcess_L3_Paths !- - - - - - - - - - - - - - - - - - - - - - - - END DO ! on waiting L4 paths END SUBROUTINE DProcess_L4_Paths SUBROUTINE DProcess_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*8, DIMENSION(3), INTENT(IN) :: uvec REAL*8, INTENT(IN) :: angle_radians, & & lr_fraction, ud_fraction LOGICAL, INTENT(IN) :: from_east INTEGER, INTENT(IN) :: font_points CHARACTER*(*), INTENT(IN) :: text REAL*8 :: equat, lon_radians, page_angle_radians, x_meters, x_points, xe_meters, & & xe_points, y_meters, y_points, ye_meters, ye_points REAL*8, DIMENSION(3) :: east_post CALL DProject (uvec = uvec, x = x_meters, y = y_meters) IF (from_east) THEN IF (DABS(uvec(3)) >= 1.D0) 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 = DSQRT(uvec(1)**2 + uvec(2)**2) lon_radians = DATAN2(uvec(2),uvec(1)) east_post(1) = equat * DCOS(lon_radians + radians_per_degree) east_post(2) = equat * DSIN(lon_radians + radians_per_degree) east_post(3) = uvec(3) ! project both points to the (x,y) system: CALL DProject (uvec = east_post, guide = uvec, x = xe_meters, y = ye_meters) ! map both points onto the page: CALL DMeters_2_Points (x_meters,y_meters, x_points,y_points) CALL DMeters_2_Points (xe_meters,ye_meters, xe_points,ye_points) page_angle_radians = DATAN2((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 DProcess_L3_Text (x_meters,y_meters, page_angle_radians, & & font_points, lr_fraction, ud_fraction, & & text) END SUBROUTINE DProcess_L4_Text SUBROUTINE DProcess_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 (up to) 4 predefined limits ! (all of which are great or small circles or half-circles), ! collectively known as the "rind" (projectable area), ! the vertices and segments outside are eliminated. ! At the end, the path(s) are passed to DProcess_L4_Paths for ! subdivision into small pieces, conversion to Bezier curves, ! and projection onto the projection plane, etc. ! ! This subprogram was completely rewritten for MODULE DMap_Projections ! in 2015. The old version was very long and complex, because it ! tried to clip paths and areas exactly at the edge of the rind ! (and add extra points and segments to achieve this). Unfortunately, ! this long code was buggy, and in a complex plot (like contour ! map of a global .grd file) it had a high chance of throwing ! an abend that would stop the whole program. ! ! This new version is much shorter and faster, and (I hope) free ! from abends. The cost is that it leaves the plotted area with ! a slightly ragged boundary, just inside the edge of the rind. ! (However, most users already choose plot parameters so as to ! hide these areas, anyway.) ! ! Note that the old version is still present, under the new name of ! DOld_Complex_Process_L5_Paths. This is now called by ! DEnd_L45_Path, but ONLY if OPTIONAL argument "retro" = .TRUE. ! This method is still used to plot melon-wedges in lower focal ! hemispheres of fault-plane-solutions. IMPLICIT NONE CHARACTER*1 :: byte INTEGER :: i, path, segments LOGICAL :: all_segments_in, all_segments_out, closed, continuation, stroked, filled LOGICAL, DIMENSION(0:ai_longest) :: vertex_is_in LOGICAL, DIMENSION(1:ai_longest) :: segment_is_in REAL*8, DIMENSION(0:ai_longest) :: fromCutPole_RC_azimuth_radians REAL*8, DIMENSION(3) :: goal_uvec, pole_uvec, start_uvec, vertex_uvec DO WHILE (ai_Ln_paths(5) > 0) path = DNext_Path(level = 5) closed = ai_closed(path) stroked = ai_stroked(path) filled = ai_filled(path) segments = ai_segments(path) ! Shortcut complex logic (below) if path has only an initial point, ! but no segments. Pass it through (even though it will have no ! graphical effect) if it falls within the rind. IF (segments == 0) THEN ! path has only an initial point start_uvec(1:3) = ai_pathlib(1:3, 0, path) byte = DIn_Rind(start_uvec) IF ((byte == 'I').OR.(byte == 'B')) THEN !Transform path to level 4, and call for Processing: ai_path_level(path) = 4 ai_Ln_paths(5) = ai_Ln_paths(5) - 1 ai_Ln_paths(4) = ai_Ln_paths(4) + 1 CALL DProcess_L4_Paths() ELSE ! just eliminate this path: ai_path_level(path) = 0 ai_Ln_paths(5) = ai_Ln_paths(5) - 1 ai_total_paths = ai_total_paths - 1 END IF ! inside, or outside CYCLE ! so as to consider the next L5 path in library END IF ! path had no segments, but just an initial point !--------------------------------------------------------- !Begin normal processing (path with at least 1 segment): !Decide which vertices are in rind, vs. outside: DO i = 0, segments IF (i == 0) THEN vertex_uvec(1:3) = ai_pathlib(1:3, 0, path) ELSE IF (ai_bent(i, path)) THEN ! a "Small_to_" vertex_uvec(1:3) = ai_pathlib(4:6, i, path) ELSE ! a "Great_to_" vertex_uvec(1:3) = ai_pathlib(1:3, i, path) END IF byte = DIn_Rind(vertex_uvec) vertex_is_in(i) = ((byte == 'I').OR.(byte == 'B')) IF (mp_using_cut) THEN ! pre-measure azimuth from the cut-pole: fromCutPole_RC_azimuth_radians(i) = DRelative_Compass (from_uvec = mp_cutPole_uvec, & & to_uvec = vertex_uvec) !and note that DRelative_Compass never abends, even for 2 == uvecs. END IF ! mp_using_cut, in this map projection? END DO ! i = 0, segments (are vertices inside or outside?) !Decide which segments are (completely) within the rind: DO i = 1, segments segment_is_in(i) = (vertex_is_in(i-1).AND.vertex_is_in(i)) END DO ! i = 1, segments (are segments within the rind?) !---------------------------------------------------------------- !Detect 2 simple cases: all segments in? Or, all segments out? all_segments_in = .TRUE. ! just initializing; may change below... all_segments_out = .TRUE. ! just initializing; may change below... DO i = 1, segments IF (segment_is_in(i)) THEN all_segments_out = .FALSE. ELSE ! this segment is out all_segments_in = .FALSE. END IF END DO ! i = 1, segments !---------------------------------------------------------------- ! If mp_using_cut, check for segments that would span across it: IF (mp_using_cut) THEN IF (.NOT.all_segments_out) THEN ! some are "in" and need checking: DO i = 1, segments !Is starting point on the "cut" side of the planet? IF (DCOS(fromCutPole_RC_azimuth_radians(i-1) - mp_intoCut_RC_azimuth_radians) > 0.0D0) THEN !Is ending point on the "cut" side of the planet? IF (DCOS(fromCutPole_RC_azimuth_radians(i) - mp_intoCut_RC_azimuth_radians) > 0.0D0) THEN !Test for case where start and end are on opposite sides of cut: IF ((DSIN(fromCutPole_RC_azimuth_radians(i-1) - mp_intoCut_RC_azimuth_radians) * & & DSIN(fromCutPole_RC_azimuth_radians(i) - mp_intoCut_RC_azimuth_radians)) <= 0.0D0) THEN !Eliminate this segment from map: segment_is_in(i) = .FALSE. all_segments_in = .FALSE. END IF ! bad, cut-spanning segment! END IF ! ending point is on "cut" side of planet END IF ! starting point is on "cut" side of planet END DO ! i = 1, segments END IF ! (.NOT.all_segments_out), so some segments are "in"(?) END IF ! mp_using_cut !---------------------------------------------------------------- ! Finally, take appropriate action on this L5 path: IF (all_segments_in) THEN ! just reclassify path as Level-4: ai_path_level(path) = 4 ai_Ln_paths(5) = ai_Ln_paths(5) - 1 ai_Ln_paths(4) = ai_Ln_paths(4) + 1 ELSE IF (all_segments_out) THEN ! just eliminate this path: ai_path_level(path) = 0 ai_Ln_paths(5) = ai_Ln_paths(5) - 1 ai_total_paths = ai_total_paths - 1 ELSE ! path is partly-in, but partly-out: !Copy any in-rind segments of L5 path to new L4 path(s): !(Note that in-rind segments are copied, and perhaps ! regrouped into a greater number of L4 paths, ! but the segments themselves are not modified.) closed = .FALSE. ! regardless of original intent; segments have been removed! DO i = 1, segments IF (segment_is_in(i)) THEN continuation = ai_in_path ! copying to local variable (NOT changed by DNew_L45_Path) IF (continuation) THEN ! continue path with one more segment IF (ai_bent(i, path)) THEN pole_uvec(1:3) = ai_pathlib(1:3, i, path) goal_uvec(1:3) = ai_pathlib(4:6, i, path) CALL DSmall_to_L45(pole_uvec, goal_uvec) ELSE goal_uvec(1:3) = ai_pathlib(1:3, i, path) CALL DGreat_To_L45(goal_uvec) END IF ELSE ! not a continuation; begin new path on L4: IF (i == 1) THEN start_uvec(1:3) = ai_pathlib(1:3, 0, path) ELSE IF (ai_bent(i-1, path)) THEN start_uvec(1:3) = ai_pathlib(4:6, i-1, path) ELSE ! .NOT.ai_bent(i-1, path), so i>1 and previous segment was a Great_to start_uvec(1:3) = ai_pathlib(1:3, i-1, path) END IF CALL DNew_L45_Path(4, start_uvec) END IF ELSE ! .NOT.segment_is_in(i); segment is out. IF (ai_in_path) THEN ! terminate previous L4 path CALL DEnd_L45_Path(closed, stroked, filled) END IF ! (no ELSE is needed here) END IF ! segment is in, or out END DO ! i = 1, segments IF (ai_in_path) CALL DEnd_L45_Path(closed, stroked, filled) !Now, delete original L5 path to clean-up: ai_path_level(path) = 0 ai_Ln_paths(5) = ai_Ln_paths(5) - 1 ai_total_paths = ai_total_paths - 1 END IF ! all_segments_in, all_segments_out, or a mixture CALL DProcess_L4_Paths() END DO ! WHILE (ai_Ln_paths(5) > 0) END SUBROUTINE DProcess_L5_Paths SUBROUTINE DProcess_L5_Text (uvec, angle_radians, from_east, & & font_points, lr_fraction, ud_fraction, & & text) ! Check whether uvec is DIn_Rind, and if so, pass on to L4. IMPLICIT NONE REAL*8, DIMENSION(3), INTENT(IN) :: uvec REAL*8, 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 = DIn_Rind (uvec) IF ((result == 'I').OR.(result == 'B')) THEN CALL DProcess_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 DProcess_L5_Text SUBROUTINE DProcess_L6_Paths (retro) ! Apply simple in-place (theta,phi) -> uvec transformation, ! and pass on to next level. IMPLICIT NONE LOGICAL, INTENT(IN), OPTIONAL :: retro INTEGER :: i, path REAL*8 :: phi, phi2, theta, theta2 REAL*8, DIMENSION(3) :: uvec DO WHILE (ai_Ln_paths(6) > 0) path = DNext_Path(level = 6) theta = ai_pathlib(1,0,path) phi = ai_pathlib(2,0,path) CALL DThetaPhi_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 DThetaPhi_2_Uvec (theta2,phi2, uvec) ai_pathlib(4:6,i,path) = uvec END IF CALL DThetaPhi_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 IF (PRESENT(retro)) THEN IF (retro) THEN CALL DOld_Complex_Process_L5_Paths() ELSE CALL DProcess_L5_Paths() END IF ELSE CALL DProcess_L5_Paths() END IF END DO ! while L6 paths remain END SUBROUTINE DProcess_L6_Paths SUBROUTINE DProcess_L6_Text (theta, phi, angle_radians, from_east, & & font_points, lr_fraction, ud_fraction, & & text) ! Transform position coordinates, and pass on. IMPLICIT NONE REAL*8, 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*8, DIMENSION(3) :: uvec CALL DThetaPhi_2_Uvec (theta, phi, uvec) CALL DProcess_L5_Text (uvec, angle_radians, from_east, & & font_points, lr_fraction, ud_fraction, & & text) END SUBROUTINE DProcess_L6_Text SUBROUTINE DProcess_L7_Paths (retro) ! Apply simple in-place (lon,lat) -> (theta,phi) transformation, ! and pass on to next level. IMPLICIT NONE LOGICAL, INTENT(IN), OPTIONAL :: retro INTEGER :: i, path REAL*8 :: lat, lon, phi, theta DO WHILE (ai_Ln_paths(7) > 0) path = DNext_Path(level = 7) lon = ai_pathlib(1,0,path) lat = ai_pathlib(2,0,path) CALL DLonLat_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 DLonLat_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 DLonLat_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 IF (PRESENT(retro)) THEN CALL DProcess_L6_Paths(retro) ELSE CALL DProcess_L6_Paths END IF END DO ! while L7 paths remain END SUBROUTINE DProcess_L7_Paths SUBROUTINE DProcess_L7_Text (lon, lat, angle_radians, from_east, & & font_points, lr_fraction, ud_fraction, & & text) ! Transform position coordinates and pass on. IMPLICIT NONE REAL*8, 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*8 :: theta, phi CALL DLonLat_2_ThetaPhi (lon, lat, theta, phi) CALL DProcess_L6_Text (theta, phi, angle_radians, from_east, & & font_points, lr_fraction, ud_fraction, & & text) END SUBROUTINE DProcess_L7_Text SUBROUTINE DProject (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 DMap_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 DProject to test inputs, ! give warning messages, create abends for projections ! to infinity or the complex plane, et cetera! ! Other programs (especially DProcess_L5_Paths) will ! guard against ever calling DProject for troublesome ! points outside the walls of the "rind" of projectable ! surface points on the sphere. Therefore, DProject ! can concentrate on simplicity and speed. !================================================================ IMPLICIT NONE REAL*8, DIMENSION(3), INTENT(IN) :: uvec REAL*8, DIMENSION(3), INTENT(IN), OPTIONAL :: guide REAL*8, INTENT(OUT) :: x, y REAL*8 :: alpha, arc_radians, argument_of_DLOG, argument_of_DTAN, 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 = DDot (uvec, mp_pole_1_uvec) dot2 = DDot (uvec, mp_pole_2_uvec) dot3 = DDot (uvec, mp_pole_3_uvec) SELECT CASE (mp_projection_number) CASE (1) ! Mercator east_radians = DATAN2(dot2, dot1) equat = DSQRT(dot1**2 + dot2**2) IF (PRESENT(guide)) THEN dot1 = DDot (guide, mp_pole_1_uvec) dot2 = DDot (guide, mp_pole_2_uvec) guide_east = DATAN2(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 = DATAN2(dot3, equat) argument_of_DTAN = (up_radians * 0.5D0) + Pi_over_4 IF (ABS(DCOS(argument_of_DTAN)) > 1.0D-10) THEN ! angle is far enough away from Pi_over_2 (or -Pi_over_2) so that DTAN() can be safely evaluated. argument_of_DLOG = DTAN(argument_of_DTAN) ELSE IF (argument_of_DTAN > 0.0D0) THEN ! close to Pi_over_2; DTAN would be +inf. argument_of_DLOG = HUGE(Pi) ELSE ! close to -Pi_over_2; DTAN would be -inf. argument_of_DLOG = -HUGE(Pi) END IF IF (argument_of_DLOG > 0.0D0) THEN beta = mp_radius_meters * DLOG(argument_of_DLOG) ELSE beta = -HUGE(Pi) END IF ! 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 = DSQRT(dot1**2 + dot2**2) latitude_radians = DATAN2(dot3, equat) IF (equat > 0.D0) THEN longitude_radians = DATAN2(dot2, dot1) ELSE longitude_radians = mp_proj2_lambda0 END IF IF (PRESENT(guide)) THEN dot1 = DDot (guide, mp_pole_1_uvec) dot2 = DDot (guide, mp_pole_2_uvec) equat2 = dot1**2 + dot2**2 IF (equat2 > 0.D0) THEN guide_east = DATAN2(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 (DABS(latitude_radians - Pi_over_2) > 0.001D0) THEN rho = mp_radius_meters * mp_proj2_F / & & DTAN((Pi_over_2 + latitude_radians)/2.D0)**mp_proj2_n ELSE ! avoid DTAN(Pi_over_2) ---> Domain Error (crash) rho = 0.D0 END IF alpha = rho * DSIN(theta) beta = mp_proj2_rho0 - rho * DCOS(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 = DSQRT(dot1**2 + dot2**2) latitude_radians = DATAN2(dot3, equat) IF (equat > 0.D0) THEN longitude_radians = DATAN2(dot2, dot1) ELSE longitude_radians = mp_proj3_lambda0 END IF IF (PRESENT(guide)) THEN dot1 = DDot (guide, mp_pole_1_uvec) dot2 = DDot (guide, mp_pole_2_uvec) equat2 = dot1**2 + dot2**2 IF (equat2 > 0.D0) THEN guide_east = DATAN2(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 * DSQRT(MAX(0.D0, mp_proj3_C - 2.D0 * mp_proj3_n * DSIN(latitude_radians))) / mp_proj3_n alpha = rho * DSIN(theta) beta = mp_proj3_rho0 - rho * DCOS(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 = DSQRT(dot1**2 + dot2**2) latitude_radians = DATAN2(dot3, equat) IF (equat > 0.D0) THEN longitude_radians = DATAN2(dot2, dot1) ELSE longitude_radians = mp_proj4_lon0_radians END IF IF (PRESENT(guide)) THEN dot1 = DDot (guide, mp_pole_1_uvec) dot2 = DDot (guide, mp_pole_2_uvec) equat2 = dot1**2 + dot2**2 IF (equat2 > 0.D0) THEN guide_east = DATAN2(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) * DSIN(latitude_radians) IF (DABS(latitude_radians - Pi_over_2) > 0.0001D0) THEN alpha = mp_radius_meters * DSIN(E) / DTAN(latitude_radians) ELSE ! avoid DTAN(Pi_over_2) ---> Domain error (crash)! alpha = 0.D0 END IF beta = mp_radius_meters * (latitude_radians - mp_proj4_lat0_radians & & + (1.D0 - DCOS(E)) / DTAN(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 = DSQRT(dot1**2 + dot2**2) latitude_radians = DATAN2(dot3, equat) IF (equat > 0.D0) THEN longitude_radians = DATAN2(dot2, dot1) ELSE longitude_radians = 0.D0 END IF IF (PRESENT(guide)) THEN dot1 = DDot (guide, mp_pole_1_uvec) dot2 = DDot (guide, mp_pole_2_uvec) equat2 = dot1**2 + dot2**2 IF (equat2 > 0.D0) THEN guide_east = DATAN2(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 (DABS((latitude_radians - mp_proj5_lat0_radians) - Pi_over_2) > 0.0001D0) THEN rho = mp_proj5_rtan - mp_radius_meters * DTAN(latitude_radians - mp_proj5_lat0_radians) ELSE ! avoid DTAN(Pi_over_2) ---> Domain error (crash)! rho = 0.D0 END IF alpha = rho * DSIN(theta) beta = mp_proj5_ypole - rho * DCOS(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.D0) THEN azimuth = DATAN2 (dot2, dot3) arc_radians = DArc (mp_pole_1_uvec, uvec) radius_meters = 2.D0 * mp_radius_meters * DTAN(0.5D0 * arc_radians) xpp = radius_meters * DSIN(azimuth) ypp = radius_meters * DCOS(azimuth) ELSE xpp = 0.0D0 ypp = 0.0D0 END IF CASE (7) ! Lambert Azimuthal Equal-Area equat2 = dot2**2 + dot3**2 IF (equat2 > 0.D0) THEN azimuth = DATAN2 (dot2, dot3) arc_radians = DArc (mp_pole_1_uvec, uvec) radius_meters = 2.D0 * mp_radius_meters * DSIN(0.5D0 * arc_radians) xpp = radius_meters * DSIN(azimuth) ypp = radius_meters * DCOS(azimuth) ELSE xpp = 0.0D0 ypp = 0.0D0 END IF CASE (8) ! Azimuthal Equidistant equat2 = dot2**2 + dot3**2 IF (equat2 > 0.D0) THEN azimuth = DATAN2 (dot2, dot3) arc_radians = DArc (mp_pole_1_uvec, uvec) radius_meters = mp_radius_meters * arc_radians xpp = radius_meters * DSIN(azimuth) ypp = radius_meters * DCOS(azimuth) ELSE xpp = 0.0D0 ypp = 0.0D0 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.D0) THEN azimuth = DATAN2 (dot2, dot3) arc_radians = DArc (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 * DTAN(arc_radians) xpp = radius_meters * DSIN(azimuth) ypp = radius_meters * DCOS(azimuth) ELSE xpp = 0.0D0 ypp = 0.0D0 END IF CASE DEFAULT WRITE (*,"(' ERROR: Map projection [',A,'] is not supported.')") mp_projection CALL DTraceback 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 DProject SUBROUTINE DReject (x_meters, y_meters, success, uvec) ! Opposite of "SUBROUTINE DProject": ! 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*8, INTENT(IN) :: x_meters, y_meters LOGICAL, INTENT(OUT) :: success REAL*8, DIMENSION(3), INTENT(OUT) :: uvec INTEGER :: i LOGICAL :: converged, converging REAL*8 :: 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*8, 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 (DABS(east_radians) > Pi) THEN success = .FALSE. RETURN END IF up_radians = Pi_over_2 - 2.D0 * DATAN(DEXP(-beta/mp_radius_meters)) IF (DABS(up_radians) > 1.221381D0) THEN ! 69.98 degrees success = .FALSE. ELSE success = .TRUE. dot1 = DCOS(up_radians) * DCOS(east_radians) dot2 = DCOS(up_radians) * DSIN(east_radians) dot3 = DSIN(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.D0, mp_proj2_n) * DSQRT(alpha**2 + (mp_proj2_rho0 - beta)**2) IF (mp_proj2_n > 0.D0) THEN theta = DATAN2(alpha, (mp_proj2_rho0 - beta)) ELSE theta = DATAN2(-alpha, (-mp_proj2_rho0 + beta)) END IF longitude_radians = (theta / mp_proj2_n) + mp_proj2_lambda0 IF (rho == 0.D0) THEN latitude_radians = SIGN(Pi_over_2, mp_proj2_n) ELSE latitude_radians = 2.D0 * DATAN((mp_radius_meters * mp_proj2_f / rho)**(1.D0/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*DCOS(latitude_radians)*DCOS(longitude_radians) + & & mp_pole_2_uvec*DCOS(latitude_radians)*DSIN(longitude_radians) + & & mp_pole_3_uvec*DSIN(latitude_radians) success = .TRUE. ! provisionally, unless changed below: IF (mp_using_top) THEN IF (DDot(uvec,mp_top_uvec) > mp_top_maxdot) success = .FALSE. END IF IF (mp_using_bottom) THEN IF (DDot(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 = DSQRT(alpha**2 + (mp_proj3_rho0 - beta)**2) IF (mp_proj3_n > 0.D0) THEN theta = DATAN2(alpha, (mp_proj3_rho0 - beta)) ELSE theta = DATAN2(-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.D0 * mp_proj3_n) IF (DABS(argument) > 1.D0) THEN success = .FALSE. ELSE ! usual case latitude_radians = DASIN(argument) ! Note: both of the the above may be in a rotated system: uvec = mp_pole_1_uvec*DCOS(latitude_radians)*DCOS(longitude_radians) + & & mp_pole_2_uvec*DCOS(latitude_radians)*DSIN(longitude_radians) + & & mp_pole_3_uvec*DSIN(latitude_radians) success = .TRUE. ! unless changed below: IF (mp_using_top) THEN IF (DDot(uvec,mp_top_uvec) > mp_top_maxdot) success = .FALSE. END IF IF (mp_using_bottom) THEN IF (DDot(uvec,mp_bottom_uvec) > mp_bottom_maxdot) success = .FALSE. END IF END IF ! DABS(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.D0 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.0D0 iterate: DO i = 1, 10 tanlat = DTAN(latitude_radians) new_lat_radians = latitude_radians - (A*(latitude_radians*tanlat + 1.D0) & & - latitude_radians - 0.5D0*(latitude_radians**2 + B)*tanlat) & & / (((latitude_radians - A)/tanlat) - 1.D0) change = DABS(new_lat_radians - latitude_radians) IF (change > 1.D0) EXIT iterate ! probably diverging converging = (change < last_change) IF (converging.AND.(change < 1.D-5)) EXIT iterate latitude_radians = new_lat_radians ! memory for next iteration last_change = change END DO iterate converged = converging .AND. (change < 1.D-4) IF (converged) THEN IF (latitude_radians < 0.D0) THEN ! wrong hemisphere success = .FALSE. ELSE IF (latitude_radians == 0.D0) 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 = ((DASIN(alpha*tanlat/mp_radius_meters)) & & / DSIN(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*DCOS(latitude_radians)*DCOS(longitude_radians) + & & mp_pole_2_uvec*DCOS(latitude_radians)*DSIN(longitude_radians) + & & mp_pole_3_uvec*DSIN(latitude_radians) IF (mp_using_left) THEN IF (DDot(uvec,mp_left_uvec) > mp_left_maxdot) success = .FALSE. END IF IF (mp_using_right) THEN IF (DDot(uvec,mp_right_uvec) > mp_right_maxdot) success = .FALSE. END IF IF (mp_using_bottom) THEN IF (DDot(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 = DSQRT(alpha**2 + (mp_proj5_ypole - beta)**2) theta = DATAN2(alpha, (mp_proj5_ypole - beta)) longitude_radians = (theta / mp_proj5_n) + mp_proj5_lambda0 IF (rho == 0.D0) THEN latitude_radians = SIGN(Pi_over_2, mp_proj5_n) ELSE latitude_radians = mp_proj5_lat0_radians + DATAN2(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*DCOS(latitude_radians)*DCOS(longitude_radians) + & & mp_pole_2_uvec*DCOS(latitude_radians)*DSIN(longitude_radians) + & & mp_pole_3_uvec*DSIN(latitude_radians) success = .TRUE. IF (mp_using_top) THEN IF (DDot(uvec,mp_top_uvec) > mp_top_maxdot) success = .FALSE. END IF IF (mp_using_bottom) THEN IF (DDot(uvec,mp_bottom_uvec) > mp_bottom_maxdot) success = .FALSE. END IF CASE (6) ! Stereographic radius_meters = DSQRT (xpp**2 + ypp**2) IF (radius_meters > 0.D0) THEN azimuth = DATAN2(xpp, ypp) arc_radians = 2.0D0 * DATAN(0.5D0 * radius_meters / mp_radius_meters) CALL DTurn_To (azimuth, mp_pole_1_uvec, arc_radians, & ! inputs & omega_uvec, uvec) ELSE azimuth = 0.D0 arc_radians = 0.D0 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 = DSQRT (xpp**2 + ypp**2) IF (radius_meters > 0.D0) THEN azimuth = DATAN2(xpp, ypp) argument = 0.5D0 * radius_meters / mp_radius_meters argument = MIN(argument, 0.99D0) ! prevent out-of-range error arc_radians = 2.0D0 * DASIN(argument) CALL DTurn_To (azimuth, mp_pole_1_uvec, arc_radians, & ! inputs & omega_uvec, uvec) ELSE azimuth = 0.D0 arc_radians = 0.D0 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 = DSQRT (xpp**2 + ypp**2) IF (radius_meters > 0.D0) THEN azimuth = DATAN2(xpp, ypp) arc_radians = radius_meters / mp_radius_meters CALL DTurn_To (azimuth, mp_pole_1_uvec, arc_radians, & ! inputs & omega_uvec, uvec) ELSE azimuth = 0.D0 arc_radians = 0.D0 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 = DSQRT (xpp**2 + ypp**2) IF (radius_meters <= mp_radius_meters) THEN success = .TRUE. dot2 = xpp / mp_radius_meters dot3 = ypp / mp_radius_meters dot1 = DSQRT (1.00D0 - 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 = DSQRT (xpp**2 + ypp**2) IF (radius_meters > 0.D0) THEN azimuth = DATAN2(xpp, ypp) arc_radians = DATAN(radius_meters / mp_radius_meters) CALL DTurn_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 DTraceback END SELECT END SUBROUTINE DReject REAL*8 FUNCTION DRelative_Compass (from_uvec, to_uvec) ! At most points, works exactly like DCompass: ! returns azimuth (in radians, clockwise from North) ! of the great-circle arc from "from_uvec" to "to_uvec", ! measured at location from_uvec. ! However, unlike DCompass, it does not crash at the N and S poles, ! where Theta and Phi are undefined. Instead, it uses ! SUBROUTINE DNorthEast_Convention to make an ! arbitrary choice of axes, so that RELATIVE azimuths can ! be measured from pivot points at the poles, by multiple calls. ! Another safety feature is that this version does NOT abend ! if to_uvec == from_uvec (or is antipodal to it); ! in this case an arbitrary azimuth of 0.0 is returned. IMPLICIT NONE REAL*8, DIMENSION(3), INTENT(IN) :: from_uvec, to_uvec REAL*8 :: ve, vn REAL*8, DIMENSION(3) :: omega_vec, omega_uvec, & & north_uvec, east_uvec, & & v_vec IF ((from_uvec(1) == to_uvec(1)).AND. & &(from_uvec(2) == to_uvec(2)).AND. & &(from_uvec(3) == to_uvec(3))) THEN !WRITE (*,"(' ERROR: Compass bearing from point to itself undefined.')") !CALL DTraceback DRelative_Compass = 0.0D0 ELSE IF ((from_uvec(1) == -to_uvec(1)).AND. & &(from_uvec(2) == -to_uvec(2)).AND. & &(from_uvec(3) == -to_uvec(3))) THEN !WRITE (*,"(' ERROR: Compass bearing from point to antipode undefined.')") !CALL DTraceback DRelative_Compass = 0.0D0 ELSE ! Normal case: CALL DCross (from_uvec, to_uvec, omega_vec) CALL DMake_Uvec(omega_vec, omega_uvec) CALL DCross (omega_uvec, from_uvec, v_vec) CALL DNorthEast_Convention (from_uvec, north_uvec, east_uvec) vn = DDot(v_vec, north_uvec) ve = DDot(v_vec, east_uvec) DRelative_Compass = DATAN2(ve, vn) END IF END FUNCTION DRelative_Compass SUBROUTINE DRestore_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 DSave_mp_State() before CALL DRestore_mp_State()')") CALL DTraceback END IF ! state was previously saved, or not END SUBROUTINE DRestore_mp_State SUBROUTINE DSave_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 DSet_Zoom before DSave_mp_State')") CALL DTraceback END IF ! a state of mp is defined, or not END SUBROUTINE DSave_mp_State SUBROUTINE DThetaPhi_2_LonLat (theta, phi, lon, lat) ! Converts from theta (co-latitude, from N pole, in radians) ! and phi (longitude, E from Greenwich, in radians) ! to "lon" (East longitude, in degrees; West is negative) ! and "lat" (North latitude, in degrees; South is negative). IMPLICIT NONE REAL*8, INTENT(IN) :: theta, phi REAL*8, INTENT(OUT) :: lon, lat lat = 90.0D0 - degrees_per_radian * DABS(theta) lat = MAX (lat, -90.0D0) lon = degrees_per_radian * phi IF (lon > 180.0D0) lon = lon - 360.0D0 IF (lon <= -180.0D0) lon = lon + 360.0D0 END SUBROUTINE DThetaPhi_2_LonLat SUBROUTINE DThetaPhi_2_Uvec (theta, phi, uvec) ! Converts from theta (co-latitude, from N pole) and ! phi (longitude, E from Greenwich) [both in radians] ! to a 3-component Cartesian unit vector, which points ! from the center of the unit sphere to a surface point. ! Its first axis outcrops at (0E, 0N). ! Its second axis outcrops at (90E, 0N). ! Its third axis outcrops at 90N. IMPLICIT NONE REAL*8, INTENT(IN) :: theta, phi REAL*8, DIMENSION(3), INTENT(OUT) :: uvec REAL*8 :: equat uvec(3) = DCOS(theta) equat = DSIN(theta) uvec(1) = equat * DCOS(phi) uvec(2) = equat * DSIN(phi) END SUBROUTINE DThetaPhi_2_Uvec SUBROUTINE DTrace_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 initial 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*8, INTENT(IN) :: from_t_wall, to_t_wall REAL*8, DIMENSION(3), INTENT(IN) :: to_point_uvec INTEGER :: current_wall, end_wall, i, start_wall LOGICAL first_pass, one_step REAL*8 :: end_f, start_f REAL*8, 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 DSmall_To_L45 (pole_uvec, to_point_uvec) ELSE ! ccw around the boundary is ccw around the anti-pole pole_uvec = -pole_uvec CALL DSmall_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 rind-edge-following small-circles 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 choose other target corner 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 DSmall_To_L45 (pole_uvec, to_point_uvec) RETURN ELSE ! go to corner point CALL DSmall_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 ! counterclockwise current_wall = mp_next_wall(1,current_wall) END IF END DO END IF END SUBROUTINE DTrace_Rind SUBROUTINE DTurn_To (azimuth_radians, base_uvec, far_radians, & ! input & omega_uvec, result_uvec) ! output ! Computes uvec "result_uvec" (a 3-component Cartesian unit ! vector from the center of the planet) which results from ! rotating along a great circle beginning at "base_uvec" ! for an angle of "far_radians", in the initial direction ! given by azimuth_radians" (clockwise, from North). ! Also returned is "omega_uvec", the pole of rotation. ! NOTE: At the poles, azimuth is undefined. Near the ! poles, it is defined but numerically unstable. Therefore, ! Turn_To uses the same SUBROUTINE DNorthEast_Convention as ! FUNCTION DRelative_Compass does, so they can work together ! to find internal points on a small circle (as in DProcess_L4_Paths). IMPLICIT NONE REAL*8, INTENT(IN) :: azimuth_radians, far_radians REAL*8, DIMENSION(3), INTENT(IN) :: base_uvec REAL*8, DIMENSION(3), INTENT(OUT) :: omega_uvec, result_uvec REAL*8, PARAMETER :: pole_width = 1.745D-4 ! 0.01 degrees REAL*8 :: complement, cos_size, e_part, n_part, sin_size REAL*8, DIMENSION(3) :: east_uvec, north_uvec, t_uvec REAL*8, DIMENSION(3,3) :: rotation_matrix CALL DNorthEast_Convention (base_uvec, north_uvec, east_uvec) e_part = -DCOS(azimuth_radians) n_part = DSIN(azimuth_radians) omega_uvec(1) = e_part*east_uvec(1) + n_part*north_uvec(1) omega_uvec(2) = e_part*east_uvec(2) + n_part*north_uvec(2) omega_uvec(3) = e_part*east_uvec(3) + n_part*north_uvec(3) cos_size = DCOS(far_radians) sin_size = DSIN(far_radians) complement = 1.00D0 - cos_size rotation_matrix(1,1) = cos_size + complement*omega_uvec(1)*omega_uvec(1) rotation_matrix(1,2) = complement*omega_uvec(1)*omega_uvec(2) - sin_size*omega_uvec(3) rotation_matrix(1,3) = complement*omega_uvec(1)*omega_uvec(3) + sin_size*omega_uvec(2) rotation_matrix(2,1) = complement*omega_uvec(2)*omega_uvec(1) + sin_size*omega_uvec(3) rotation_matrix(2,2) = cos_size + complement*omega_uvec(2)*omega_uvec(2) rotation_matrix(2,3) = complement*omega_uvec(2)*omega_uvec(3) - sin_size*omega_uvec(1) rotation_matrix(3,1) = complement*omega_uvec(3)*omega_uvec(1) - sin_size*omega_uvec(2) rotation_matrix(3,2) = complement*omega_uvec(3)*omega_uvec(2) + sin_size*omega_uvec(1) rotation_matrix(3,3) = cos_size + complement*omega_uvec(3)*omega_uvec(3) !Copy base_uvec in case user of this routine plans to change it: t_uvec = base_uvec result_uvec(1) = rotation_matrix(1,1)*t_uvec(1) + & & rotation_matrix(1,2)*t_uvec(2) + & & rotation_matrix(1,3)*t_uvec(3) result_uvec(2) = rotation_matrix(2,1)*t_uvec(1) + & & rotation_matrix(2,2)*t_uvec(2) + & & rotation_matrix(2,3)*t_uvec(3) result_uvec(3) = rotation_matrix(3,1)*t_uvec(1) + & & rotation_matrix(3,2)*t_uvec(2) + & & rotation_matrix(3,3)*t_uvec(3) END SUBROUTINE DTurn_To SUBROUTINE DUvec_2_LonLat (uvec, lon, lat) IMPLICIT NONE REAL*8, DIMENSION(3), INTENT(IN) :: uvec REAL*8, INTENT(OUT) :: lon, lat REAL*8 :: theta, phi CALL DUvec_2_ThetaPhi (uvec, theta, phi) CALL DThetaPhi_2_LonLat (theta, phi, lon, lat) END SUBROUTINE DUvec_2_LonLat SUBROUTINE DUvec_2_ThetaPhi (uvec, theta, phi) ! converts from Cartesian unit vector to theta (colatitude) ! and phi (longitude), both in radians IMPLICIT NONE REAL*8, DIMENSION(3), INTENT(IN) :: uvec REAL*8, INTENT(OUT) :: theta, phi REAL*8 :: equat, equat2 equat2 = uvec(1)*uvec(1) + uvec(2)*uvec(2) IF (equat2 == 0.0D0) THEN phi = 0.0D0 ! actually undefined; provide default 0.0D0 IF (uvec(3) > 0.0D0) THEN theta = 0.0D0 ! N pole ELSE theta = Pi ! S pole END IF ELSE equat = DSQRT(equat2) theta = DATAN2(equat, uvec(3)) phi = DATAN2(uvec(2), uvec(1)) END IF END SUBROUTINE DUvec_2_ThetaPhi REAL*8 FUNCTION DVector_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*8, DIMENSION(3), INTENT(IN) :: site_uvec, vector REAL*8, DIMENSION(3) :: theta_uvec, phi_uvec REAL*8 :: N_part, E_part CALL DLocal_Theta(site_uvec, theta_uvec) CALL DLocal_Phi (site_uvec, phi_uvec) N_part = -DDot(theta_uvec, vector) E_part = DDot(phi_uvec, vector) IF ((N_part /= 0.0D0).OR.(E_part /= 0.0D0)) THEN DVector_Azimuth = DATAN2(E_part, N_part) ELSE DVector_Azimuth = 0.0D0 END IF END FUNCTION DVector_Azimuth END MODULE DMap_Projections !===============================================