MODULE Flat_Section !Support routines for any (external) main program that creates a "section" diagram, !displaying data on one particular vertical planar cut through the lithosphere. !A "flat section" is one where the sea-level line at the top is straight, not curved. !That is, the section can be thought of as being cut down into a map-projection of !the lithosphere onto a virtual flat-Earth. ! !This module is to be used in conjunction with Peter Bird's other mapping & graphical ! MODULE Adobe_Illustrator, !because many of its routines will function without modification in this new context. !However, it contains replacements for some routines in other modules: ! MODULE Map_Projections ! MODULE Map_Tools !which were designed specifically for creating maps. Each replacement has a !modified name, including "_Section" or "FS", to avoid confusion. !State variables associated with this module all have names beginning with "fs_". !Note that module Map_Projections will typically still be needed to create a section, ! because the (lon, lat) coordinates of points along the section depend on the ! current map projection, and must be computed to look up data values. ! !Here are a few conceptual distinctions between my code for maps and my code for sections: !(1) The available coordinate systems for maps include round-Earth (lon, lat; or uvec) ! and flat-Earth (x, y, z) and graphical (points-to-right, points-up from lower left). ! The available coordinate systems for sections include one orthogonal Cartesian ! right-handed set of axes, in units of real-Earth-sized meters: ! r = "right" = horizontal distance along a section, measured from the left end. ! z = "up" = vertically upward, both in Nature and on the flat section page; ! defined as 0. at mean sea level, and positive upwards. ! v = "viewer" = orthogonal to the section, and toward the viewer. ! A matrix-based translation of vectors and tensors from map-based (x, y, z) ! (as in application FlatMaxwell) to/from flat-section (r, z, v) coordinates is ! part of the new code in this module; fortunately, the transformation law ! does not vary across the section. ! The same graphical system of (points-to-right, points-up) is still available, ! based on the definition that there are 72 points in an inch, ! and taking the graphical origin as the lower-left corner of the page. ! As in my map-making modules, there is some choice of "levels" at which to plot: ! Level 1 provides access to the whole graphical page, but only in point units. ! Level 2 also defines locations using points, but objects plotted to level 2 ! will be automatically truncated if they extend outside the rectangular ! "window" of the plane section. (This window is also the rectangular area ! that may be occupied by optional colored bitmaps in some sections.) ! Level 3 allows plotting in (r, z) coordinates in units of meters, with ! plotted objects automatically truncated to stay in the window (as on level 2). ! This module has no higher levels (4 ~ 7) corresponding to the round-Earth ! coordinates of my map-making modules. ! !By Peter Bird, UCLA, 2013.10+. !Copyright 2013 by Peter Bird and the Regents of the University of California. USE Adobe_Illustrator ! provided as Adobe_Illustrator.f90 USE Map_Projections ! provided as Map_Projections.f90 USE Map_Tools ! provided as Map_Tools.f90 IMPLICIT NONE LOGICAL :: fs_zoom_defined = .FALSE. REAL :: fs_meters_per_point, & & fs_pointsRight_2_r_constant, fs_pointsRight_2_r_partial, & & fs_pointsUp_2_z_constant, fs_pointsUp_2_z_partial, & & fs_r_2_pointsRight_constant, fs_r_2_pointsRight_partial, & & fs_scale_denominator, fs_z_2_pointsUp_constant, fs_z_2_pointsUp_partial REAL :: fs_last_FSL3_r_meters = -999999., & ! last "pen" position, & fs_last_FSL3_z_meters = -999999. ! used to decide which lineto's ! are degenerate; -999999. is "undefined" REAL, DIMENSION(3) :: fs_rzv_2_xyz_constants, fs_xyz_2_rzv_constants REAL, DIMENSION(3, 3) :: fs_rzv_2_xyz_partials, fs_xyz_2_rzv_partials CONTAINS SUBROUTINE Circle_on_FSL3 (r_meters,z_meters, radius_meters, stroke, fill) ! Draws a complete circle centered at (r_meters,z_meters) on section, ! with radius "radius_meters". ! Note that stroke width and color (if used) and ! fill color or pattern (if used) must be predefined. IMPLICIT NONE REAL, INTENT(IN) :: r_meters, z_meters, radius_meters LOGICAL, INTENT(IN) :: stroke, fill REAL :: radius_points, x_points, y_points IF (.NOT.ai_page_open) THEN WRITE (*,"(' ERROR: Cannot Circle_on_FSL3 before Begin_Page.')") CALL Traceback END IF IF (.NOT.fs_zoom_defined) THEN WRITE (*,"(' ERROR: Cannot Circle_on_FSL3 before Set_Section_Zoom.')") CALL Traceback END IF IF (ai_in_path) THEN WRITE (*,"(' ERROR: Cannot Circle_on_FSL3 with another path already open.')") CALL Traceback END IF IF (radius_meters < 0.0) THEN WRITE (*,"(' ERROR: Negative radius to Circle_on_FSL3: ',1P,E9.2)") radius_meters CALL Traceback END IF CALL FSMeters_2_Points (r_meters,z_meters, x_points,y_points) radius_points = 2834.6 * (radius_meters / fs_scale_denominator) CALL Circle_on_L12 (2, x_points,y_points, radius_points, stroke, fill) END SUBROUTINE Circle_on_FSL3 SUBROUTINE End_FSL3_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 ((fs_last_FSL3_r_meters /= ai_pathlib(1,0,path)).OR. & &(fs_last_FSL3_z_meters /= ai_pathlib(2,0,path))) THEN CALL Line_To_FSL3(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. fs_last_FSL3_r_meters = -9.E25 ! (undefined) fs_last_FSL3_z_meters = -9.E25 CALL Process_FSL3_Paths ELSE !neither stroked nor filled WRITE (*,"(' ERROR: End_FSL3_path with neither fill nor stroke.')") CALL Traceback END IF ELSE WRITE (*,"(' ERROR: Cannot End_FSL3_path when other level open:',I3)") level CALL Traceback END IF ELSE WRITE (*,"(' ERROR: Cannot End_FSL3_Path when path not open.')") CALL Traceback END IF END SUBROUTINE End_FSL3_Path SUBROUTINE FSMeters_2_Points (r_meters,z_meters, x_points,y_points) ! Takes one point from level FSL3 (section, in meters) ! to level 1 or 2 (page, in points) IMPLICIT NONE REAL, INTENT(IN) :: r_meters, z_meters REAL, INTENT(OUT) :: x_points, y_points x_points = fs_r_2_pointsRight_constant + & & fs_r_2_pointsRight_partial * r_meters y_points = fs_z_2_pointsUp_constant + & & fs_z_2_pointsUp_partial * z_meters END SUBROUTINE FSMeters_2_Points SUBROUTINE Line_To_FSL3 (r_meters, z_meters) IMPLICIT NONE REAL, INTENT(IN) :: r_meters, z_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 (r_meters == fs_last_FSL3_r_meters) THEN IF (z_meters == fs_last_FSL3_z_meters) RETURN END IF ! USUAL CASE: next = ai_segments(path) + 1 ai_segments(path) = next ai_pathlib(1:6,next,path) = 0.0 ai_pathlib(1,next,path) = r_meters ai_pathlib(2,next,path) = z_meters ai_bent(next,path) = .FALSE. fs_last_FSL3_r_meters = r_meters ! (memory) fs_last_FSL3_z_meters = z_meters ELSE WRITE (*,"(' ERROR: Increase ai_longest from current ',I6)") ai_longest CALL Traceback END IF ELSE WRITE (*,"(' ERROR: Cannot Line_To_FSL3 in path on level ',I2)") & & ai_current_path_level CALL Traceback END IF ELSE WRITE (*,"(' ERROR: Cannot Line_To_FSL3 before New_FSL3_Path.')") CALL Traceback END IF END SUBROUTINE Line_To_FSL3 SUBROUTINE New_FSL3_Path (r_meters, z_meters) ! Level FSL3 is data on the cross-section plane, which is ! an imaginary planet-sized plane with dimensions in meters. ! The coordinates are (r, z, v), where r is to Right !(from 0 at the left end of the section), z is up !(from sea level; therefore z <= 0. within section), ! and v is perpendicular to the section, toward the ! Viewer (used only for expressing stress tensors). IMPLICIT NONE REAL, INTENT(IN) :: r_meters, z_meters INTEGER :: path IF (ai_page_open) THEN IF (fs_zoom_defined) THEN IF (.NOT. ai_in_path) THEN ! USUAL CASE: ai_in_path = .TRUE. path = Next_Free_Path() ! Note: this function will check for full library. ai_current_path_index = path ai_current_path_level = 3 ai_path_level(path) = 3 ai_Ln_paths(3) = ai_Ln_paths(3) + 1 ai_total_paths = ai_total_paths + 1 ai_segments(path) = 0 ai_pathlib(1:6,0,path) = 0.0 ai_pathlib(1,0,path) = r_meters ai_pathlib(2,0,path) = z_meters fs_last_FSL3_r_meters = r_meters ! (memory) fs_last_FSL3_z_meters = z_meters ELSE ! ai_in_path already WRITE (*,"(' ERROR: Cannot New_FSL3_Path with a path open.')") CALL Traceback END IF ELSE ! zoom not defined WRITE (*,"(' ERROR: Must Set_Zoom before New_FSL3_Path.')") CALL Traceback END IF ELSE ! no page open WRITE (*,"(' ERROR: Cannot New_L3_Path before Begin_Page.')") CALL Traceback END IF END SUBROUTINE New_FSL3_Path SUBROUTINE Process_FSL3_Paths () ! Perform an in-place translation from meters (r, z) to points (right, up), ! and pass on to the next level. ! Another action performed at this point is to replace all ! curved segments which are exactly or nearly straight ! with straight-line segments. ! This action will catch all segments from higher levels, too. ! (This is necessary because in Module Adobe_Illustrator, ! the test was applied in Curve_To_L12, and these segments ! do not pass throught that user interface.) IMPLICIT NONE INTEGER :: i, path LOGICAL inline1, inline2 REAL :: crossx, crossy, dot, & & last_x_points, last_y_points, length, offline, & & ux, uy, vx, vy, & & x0_points, x1_points, x2_points, x3_points, & & y0_points, y1_points, y2_points, y3_points DO WHILE (ai_Ln_paths(3) > 0) path = Next_Path(level = 3) CALL FSMeters_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 FSMeters_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 FSMeters_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 FSMeters_2_Points(ai_pathlib(5,i,path),ai_pathlib(6,i,path),x3_points,y3_points) ai_pathlib(5,i,path) = x3_points ai_pathlib(6,i,path) = y3_points ! Simplify degenerate curveto's !(those which are straight within a tolerance ! of 1 point) by converting them to lineto's. !(This will also cover the case of handle points ! that erroneously coincide with end points.) IF ((x3_points /= last_x_points).OR. & &(y3_points /= last_y_points)) THEN ! straight line is defined; test can proceed: ux = x3_points - last_x_points uy = y3_points - last_y_points length = SQRT(ux**2 + uy**2) ! > 0. ux = ux / length ! unit vector of baseline uy = uy / length vx = x1_points - last_x_points ! lever1 vector vy = y1_points - last_y_points dot = vx * ux + vy * uy ! length of // component crossx = vx - dot * ux ! perpendicular component crossy = vy - dot * uy offline = SQRT(crossx**2 + crossy**2) inline1 = (offline <= ai_min_amplitude_points) vx = x3_points - x2_points ! -lever2 vector vy = y3_points - y2_points dot = vx * ux + vy * uy ! length of // component crossx = vx - dot * ux ! perpendicular component crossy = vy - dot * uy offline = SQRT(crossx**2 + crossy**2) inline2 = (offline <= ai_min_amplitude_points) IF (inline1.AND.inline2) THEN ! replace curved segment with straight ai_bent(i,path) = .FALSE. ai_pathlib(1,i,path) = x3_points ai_pathlib(2,i,path) = y3_points END IF END IF ! in-line test was possible last_x_points = x3_points last_y_points = y3_points ELSE ! segment presented as straight last_x_points = x1_points last_y_points = y1_points END IF ! segment presented as bent / straight END DO ! segments in path ai_Ln_paths(3) = ai_Ln_paths(3) - 1 ai_path_level(path) = 2 ai_Ln_paths(2) = ai_Ln_paths(2) + 1 CALL Process_L2_Paths END DO END SUBROUTINE Process_FSL3_Paths SUBROUTINE Section_Km_Frame (draw_box) ! Adds marginal information to a flat section, including ! km depth and length (r) scales on the left and bottom of the ! section, respectively. ! If optional argument draw_box is present and .TRUE., ! a rectangular box frame is drawn around the section window. IMPLICIT NONE LOGICAL, OPTIONAL, INTENT(IN) :: draw_box CHARACTER*1 :: c1 CHARACTER*2 :: c2 CHARACTER*3 :: c3 CHARACTER*4 :: c4 INTEGER :: dr_km, dz_km, font_size_points, i, now_at_km, steps REAL, PARAMETER :: tick_length_points = 4.0 REAL :: dr, dz, section_depth_m, section_length_m, x1p, x2p, y1p, y2p !text labels on depth and length scales and degree marks font_size_points = NINT(ai_lonlatlabel_points) ! Tick marks and numbers on bottom (length scale) CALL Set_Line_Style (width_points = 0.5, dashed = .FALSE.) CALL Set_Fill_or_Pattern (use_pattern = .FALSE., color_name = 'foreground') section_length_m = fs_scale_denominator * (ai_window_x2_points - ai_window_x1_points)/2834.6 dr = 1.0E5; dr_km = 100 steps = Int_Below(section_length_m / dr) IF (steps == 0) THEN dr = 1.0E4 dr_km = 10 steps = Int_Below(section_length_m / dr) END IF IF (steps == 0) THEN dr = 1.0E3 dr_km = 1 steps = Int_Below(section_length_m / dr) END IF y2p = ai_window_y1_points; y1p = y2p - tick_length_points CALL Begin_Group x1p = ai_window_x1_points CALL New_L12_Path(1, x1p, y2p) CALL Line_to_L12(x1p, y1p) CALL End_L12_Path(close = .FALSE., stroke = .TRUE., fill = .FALSE.) CALL Write_L1_Text (x_points = x1p, y_points = y1p, angle_radians = 0., & & font_points = font_size_points, lr_fraction = 0.5, ud_fraction = 1.0, & & text = '0') now_at_km = 0 DO i = 1, steps x1p = x1p + 2834.6 * dr / fs_scale_denominator now_at_km = now_at_km + dr_km WRITE (c4, "(I4)") now_at_km CALL New_L12_Path(1, x1p, y2p) CALL Line_to_L12(x1p, y1p) CALL End_L12_Path(close = .FALSE., stroke = .TRUE., fill = .FALSE.) CALL Write_L1_Text (x_points = x1p, y_points = y1p, angle_radians = 0., & & font_points = font_size_points, lr_fraction = 0.5, ud_fraction = 1.0, & & text = TRIM(ADJUSTL(c4))) END DO CALL Write_L1_Text (x_points = x1p, y_points = y1p, angle_radians = 0., & & font_points = font_size_points, lr_fraction = -1.2, ud_fraction = 1.0, & & text = 'km') CALL End_Group ! Tick marks and numbers on left (depth scale) CALL Set_Line_Style (width_points = 0.5, dashed = .FALSE.) CALL Set_Fill_or_Pattern (use_pattern = .FALSE., color_name = 'foreground') section_depth_m = fs_scale_denominator * (ai_window_y2_points - ai_window_y1_points)/2834.6 dz = 1.0E4; dz_km = 10 steps = Int_Below(1.0001 * section_depth_m / dz) IF (steps == 0) THEN dz = 1.0E3 dz_km = 1 steps = Int_Below(1.0001 * section_depth_m / dz) END IF x2p = ai_window_x1_points; x1p = x2p - tick_length_points CALL Begin_Group y1p = ai_window_y2_points CALL New_L12_Path(1, x2p, y1p) CALL Line_to_L12(x1p, y1p) CALL End_L12_Path(close = .FALSE., stroke = .TRUE., fill = .FALSE.) CALL Write_L1_Text (x_points = x1p, y_points = y1p, angle_radians = 0., & & font_points = font_size_points, lr_fraction = 1.0, ud_fraction = 0.4, & & text = '0') now_at_km = 0 DO i = 1, steps y1p = y1p - 2834.6 * dz / fs_scale_denominator now_at_km = now_at_km + dz_km WRITE (c4, "(I4)") now_at_km CALL New_L12_Path(1, x2p, y1p) CALL Line_to_L12(x1p, y1p) CALL End_L12_Path(close = .FALSE., stroke = .TRUE., fill = .FALSE.) CALL Write_L1_Text (x_points = x1p, y_points = y1p, angle_radians = 0., & & font_points = font_size_points, lr_fraction = 1.0, ud_fraction = 0.4, & & text = TRIM(ADJUSTL(c4))) END DO CALL Write_L1_Text (x_points = x1p, y_points = y1p, angle_radians = 0., & & font_points = font_size_points, lr_fraction = 1.0, ud_fraction = 2.0, & & text = 'km') CALL End_Group ! the bounding frame line: IF (PRESENT(draw_box).AND.draw_box) THEN CALL Set_Line_Style (width_points = 0.5, dashed = .FALSE.) CALL Set_Stroke_Color ('foreground') CALL New_L12_Path (1, ai_window_x1_points, ai_window_y1_points) CALL Line_To_L12 (ai_window_x2_points, ai_window_y1_points) CALL Line_To_L12 (ai_window_x2_points, ai_window_y2_points) CALL Line_To_L12 (ai_window_x1_points, ai_window_y2_points) CALL Line_To_L12 (ai_window_x1_points, ai_window_y1_points) CALL End_L12_Path (close = .TRUE., stroke = .TRUE., fill = .FALSE.) END IF ! draw_box END SUBROUTINE Section_Km_Frame SUBROUTINE Set_Section_Zoom (scale_denominator, & & section_x1_m, section_y1_m, & & section_x2_m, section_y2_m, & & x1_points, y2_points) !Defines relationships between 3 different coordinate systems that !are relevant to a flat section: ! *(x, y[, z]) system, in meters, of the map projection plane; ! *(r, z[, v]) system, in meters, of the flat section plane; ! *(points-right, points-up) system of the page language. !(Note that the relationship between (x, y[, z]) and (lon, lat) ! is set elsewhere, by routines in MODULE Map_Projections.) !The input coordinate pairs (section_x1_m, section_y1_m) ! and (section_x2_m, section_y2_m) are in the map-projection- ! plane (x, y[, z]) system, and identify the left and right ! ends of the flat section on the map. !In the section-specific (r, z[, v]) system (also in meters), ! r begins at 0. at the left end and increases right; ! z is 0. at sea level and increases upward; ! v is 0. in the section plane and increases toward the viewer. !Coordinates in page language are in points (72 points = 1 inch) ! and are measured from the lower-left corner of the page, ! not considering any unprintable margins. !Input arguments (x1_points, y2_points) are in this page-language ! system, and identify the printed location of (r = 0, z = 0). IMPLICIT NONE REAL, INTENT(IN) :: scale_denominator, & & section_x1_m, section_y1_m, & & section_x2_m, section_y2_m, & & x1_points, y2_points REAL :: section_length_m !Store values in module-level variables: fs_zoom_defined = .TRUE. fs_scale_denominator = scale_denominator fs_meters_per_point = scale_denominator / (72. * 39.370079) !-------------------------------------------------------------- !Matrix for converting (x, y, z) to (r, z, v): section_length_m = SQRT((section_x1_m - section_x2_m)**2 + (section_y1_m - section_y2_m)**2) !coefficients of (x, y, & z) in building of r: fs_xyz_2_rzv_partials(1, 1) = (section_x2_m - section_x1_m)/section_length_m fs_xyz_2_rzv_partials(1, 2) = (section_y2_m - section_y1_m)/section_length_m fs_xyz_2_rzv_partials(1, 3) = 0. !coefficents of (x, y, & z) in building of z: fs_xyz_2_rzv_partials(2, 1) = 0. fs_xyz_2_rzv_partials(2, 2) = 0. fs_xyz_2_rzv_partials(2, 3) = 1. !coefficients of (x, y, & z) in building of v: fs_xyz_2_rzv_partials(3, 1) = (section_y2_m - section_y1_m)/section_length_m fs_xyz_2_rzv_partials(3, 2) = (section_x1_m - section_x2_m)/section_length_m fs_xyz_2_rzv_partials(3, 3) = 0. fs_xyz_2_rzv_constants(1) = 0. - fs_xyz_2_rzv_partials(1, 1) * section_x1_m - & & fs_xyz_2_rzv_partials(1, 2) * section_y1_m! = r when (x = 0, y = 0, z = 0) fs_xyz_2_rzv_constants(2) = 0. ! = z when (x = 0, y = 0, z = 0) fs_xyz_2_rzv_constants(3) = 0. - fs_xyz_2_rzv_partials(3, 1) * section_x1_m - & & fs_xyz_2_rzv_partials(3, 2) * section_y1_m ! = v when (x = 0, y = 0, z = 0) !----------------------------------------------------------------------------------------------------------- !Matrix for converting (r, z, v) to (x, y, z): !coefficients of (r, z, & v) in building of x: fs_rzv_2_xyz_partials(1, 1) = (section_x2_m - section_x1_m)/section_length_m fs_rzv_2_xyz_partials(1, 2) = 0. fs_rzv_2_xyz_partials(1, 3) = (section_y2_m - section_y1_m)/section_length_m !coefficents of (r, z, & v) in building of y: fs_rzv_2_xyz_partials(2, 1) = (section_y2_m - section_y1_m)/section_length_m fs_rzv_2_xyz_partials(2, 2) = 0. fs_rzv_2_xyz_partials(2, 3) = (section_x1_m - section_x2_m)/section_length_m !coefficients of (r, z, & v) in building of z: fs_rzv_2_xyz_partials(3, 1) = 0. fs_rzv_2_xyz_partials(3, 2) = 1. fs_rzv_2_xyz_partials(3, 3) = 0. fs_rzv_2_xyz_constants(1) = section_x1_m ! = x when (r = 0, z = 0, v = 0) fs_rzv_2_xyz_constants(2) = section_y1_m ! = y when (r = 0, z = 0, v = 0) fs_rzv_2_xyz_constants(3) = 0. ! = z when (r = 0, z = 0, v = 0) !------------------------------------------------------------------------------ !Conversion from r to points-right: fs_r_2_pointsRight_partial = 2834.6 / scale_denominator fs_r_2_pointsRight_constant = x1_points !Convertsion from z to points-up: fs_z_2_pointsUp_partial = 2834.6 / scale_denominator fs_z_2_pointsUp_constant = y2_points !Conversion from pointsRight to r: fs_pointsRight_2_r_partial = scale_denominator / 2834.6 fs_pointsRight_2_r_constant = 0. - fs_pointsRight_2_r_partial * x1_points !Conversion from pointsUp to z: fs_pointsUp_2_z_partial = scale_denominator / 2834.6 fs_pointsUp_2_z_constant = 0. - fs_pointsUp_2_z_partial * y2_points END SUBROUTINE Set_Section_Zoom SUBROUTINE Stress_in_Section (r, z, s_rr, s_rz, s_zz, s_vv, & & ref_pressure_SI, ref_diameter_points) ! Plots a single 3-axis symbol at (r,z) for the state of: ! * stress [in Pa], OR ! * stress anomaly [(stress + reference pressure), in Pa]. ! Location coordinates (r, z) are in Earth-scale meters, with: ! r increasing to right from the left end of the section; and ! z increasing upward from sea level (so, z <= 0. in most CALLS). ! Draws on Flat Section Level 3 (FSL3), so symbols will be truncated ! or omitted if they extend/fall outside the graphical window. ! The stress components provided must be in the (r, z, v) coordinate ! system, where v is orthogonal to the section and points toward ! the viewer. (This is a right-handed Cartesian system.) ! The convention is that positive normal stress components indicate ! (relative) tension, and are shown by diverging arrows. ! Negative normal stress components indicate (relative) compression ! and are shown with converging arrows. ! Negative or (relatively) compressive s_vv stress is shown ! by circles. ! Positive or (relatively) tensile s_vv stress is shown by ! by equilateral triangles; the radius from the center of the triangle ! to any vertex is comparable to the radius of the circle used for ! negative s_vv stresses. ! Positive shear stress s_rz tends to elongate bodies in the +r/+z ! (upper right) quadrant and shorten them in the +r/-z quadrant. ! Note that components of shear traction on the section plane are ! not considered or plotted by this routine; you could add them ! with another CALL to Vector_in_Section. ! The scaling of the symbols is given by the parameters ! ref_pressure_SI and ref_diameter_points. ! An (integrated?) isotropic pressure (anomaly?) of ref_pressure_SI ! [ units are generically "SI" for either Pa or N/m ] ! is always plotted with a symbol of diameter ref_diameter_points. ! [ One point is 1/72 inch. ] IMPLICIT NONE REAL, INTENT(IN) :: s_rr, s_rz, s_zz, s_vv, ref_pressure_SI, ref_diameter_points, r, z REAL :: from_r, from_z, meters_per_point, radius, s1, s2, to_r, to_z, u1r, u1z, u2r, u2z meters_per_point = (3.527777E-4) * fs_scale_denominator CALL Begin_Group ! Section-perpendicular normal stress (s_vv): IF (s_vv < 0.) THEN ! circle for compression radius = (0.5*ref_diameter_points)*(-s_vv/ABS(ref_pressure_Si)) radius = radius * meters_per_point CALL Circle_on_FSL3 (r, z, radius, .TRUE., .FALSE.) ELSE IF (s_vv > 0.) THEN ! equilateral triangle for tension radius = (0.5*ref_diameter_points)*(s_vv/ABS(ref_pressure_Si)) radius = radius * meters_per_point CALL New_FSL3_Path (r, z + radius) CALL Line_to_FSL3 (r - 0.8660 * radius, z - 0.5 * radius) CALL Line_to_FSL3 (r + 0.8660 * radius, z - 0.5 * radius) CALL Line_to_FSL3 (r, z + radius) CALL End_FSL3_Path (close = .TRUE., stroke = .TRUE., fill = .FALSE.) END IF ! circle or triangle for section-perpendicular stress ! Find principal stresses in section plane: CALL Principal_Axes_22 (s_rr, s_rz, s_zz, & & s1, s2, u1r,u1z, u2r,u2z) ! Plot sigma1 (in 2-D section plane) axis radius = (0.5*ref_diameter_points)*(ABS(s1)/ABS(ref_pressure_Si)) radius = radius * meters_per_point IF (s1 < 0.) THEN ! compression: converging arrows from_r = r + radius * u1r from_z = z + radius * u1z CALL Vector_in_Section (from_r, from_z, r, z) from_r = r - (from_r - r) from_z = z - (from_z - z) CALL Vector_in_Section (from_r, from_z, r, z) ELSE IF (s1 > 0.) THEN ! extension; diverging arrows to_r = r + radius * u1r to_z = z + radius * u1z CALL Vector_in_Section (r, z, to_r, to_z) to_r = r - (to_r - r) to_z = z - (to_z - z) CALL Vector_in_Section (r, z, to_r, to_z) END IF ! Plot sigma2 (in 2-D section plane) axis radius = (0.5*ref_diameter_points)*(ABS(s2)/ABS(ref_pressure_Si)) radius = radius * meters_per_point IF (s2 < 0.) THEN ! compression: converging arrows from_r = r + radius * u2r from_z = z + radius * u2z CALL Vector_in_Section (from_r, from_z, r, z) from_r = r - (from_r - r) from_z = z - (from_z - z) CALL Vector_in_Section (from_r, from_z, r, z) ELSE IF (s2 > 0.) THEN ! extension; diverging arrows to_r = r + radius * u2r to_z = z + radius * u2z CALL Vector_in_Section (r, z, to_r, to_z) to_r = r - (to_r - r) to_z = z - (to_z - z) CALL Vector_in_Section (r, z, to_r, to_z) END IF CALL End_Group END SUBROUTINE Stress_in_Section SUBROUTINE Tensor_xyz_2_rzv(xyz_tensor, rzv_tensor) !Makes use of coordinate transformation previously defined by !Set_Section_Zoom to convert a symmetric tensor (e.g., stress) !from (x, y, z) coordinates to (r, z, v) coordinates. IMPLICIT NONE REAL, DIMENSION(3, 3), INTENT(IN) :: xyz_tensor REAL, DIMENSION(3, 3), INTENT(OUT) :: rzv_tensor INTEGER :: i, j, m, n !Refer to Peter Bird's notes "Equilibrium" for ESS 246 "Stress in the Lithosphere" ! and realize that matrix R of those notes is replaced by fs_xyz_2_rzv_partials. rzv_tensor = 0. ! whole matrix; initializing before sums below: DO i = 1, 3 DO j = 1, 3 DO m = 1, 3 DO n = 1, 3 rzv_tensor(i, j) = rzv_tensor(i, j) + fs_xyz_2_rzv_partials(i, m) * xyz_tensor(m, n) * fs_xyz_2_rzv_partials(j, n) END DO END DO END DO END DO END SUBROUTINE Tensor_xyz_2_rzv SUBROUTINE Vector_in_Section (from_r,from_z, to_r,to_z) ! Draws simple vector in a section (from_r,from_z)----->(to_r,to_z). ! Coordinates (r, z) are in Earth-scale meters, with r increasing ! from left end of section, and z increasing up from sea level, respectively. ! NOTE: Only writes to Flat Section Level 3 (FSL3), because Map_Tools ! contains earlier Vector_in_Plane which I can ! CALLed for non-section uses. IMPLICIT NONE REAL, INTENT(IN) :: from_r, from_z, to_r, to_z REAL, PARAMETER :: headlength = 0.1666667, headhalfwidth = 0.0666667 REAL :: dr, dz, leftr, leftz, rightr, rightz CALL Begin_Group ! shaft of arrow CALL New_FSL3_Path (from_r,from_z) CALL Line_to_FSL3 (to_r,to_z) CALL End_FSL3_Path(close = .FALSE., stroke = .TRUE., fill = .FALSE.) ! arrowhead: dr = to_r - from_r dz = to_z - from_z leftr = to_r - headlength * dr - headhalfwidth * dz leftz = to_z - headlength * dz + headhalfwidth * dr rightr = to_r - headlength * dr + headhalfwidth * dz rightz = to_z - headlength * dz - headhalfwidth * dr CALL New_FSL3_Path (leftr,leftz) CALL Line_to_FSL3 (to_r,to_z) CALL Line_to_FSL3 (rightr,rightz) CALL End_FSL3_Path(close = .FALSE., stroke = .TRUE., fill = .FALSE.) CALL End_Group END SUBROUTINE Vector_in_Section END MODULE Flat_Section