MODULE DFlat_Section
!
!Initial letter "D" indicated DOUBLE PRECISION (or REAL*8) version; all CALL arguments
!are converted to REAL*8, and "D" is added to all FUNCTION/SUBROUTINE names. 2015.07.27.
!
!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 DAdobe_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 DMap_Projections
! MODULE DMap_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.
!Upgrade to REAL*8 ("D" version) 2015.07.27.
USE DAdobe_Illustrator ! provided as file DAdobe_Illustrator.f90
USE DMap_Projections ! provided as file DMap_Projections.f90
USE DMap_Tools ! provided as file DMap_Tools.f90
IMPLICIT NONE
LOGICAL :: fs_zoom_defined = .FALSE.
REAL*8 :: 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*8 :: fs_last_FSL3_r_meters = -999999.D0, & ! last "pen" position,
& fs_last_FSL3_z_meters = -999999.D0 ! used to decide which lineto's
! are degenerate; -999999.D0 is "undefined"
REAL*8, DIMENSION(3) :: fs_rzv_2_xyz_constants, fs_xyz_2_rzv_constants
REAL*8, DIMENSION(3, 3) :: fs_rzv_2_xyz_partials, fs_xyz_2_rzv_partials
CONTAINS
SUBROUTINE DCircle_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*8, INTENT(IN) :: r_meters, z_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_FSL3 before DBegin_Page.')")
CALL DTraceback
END IF
IF (.NOT.fs_zoom_defined) THEN
WRITE (*,"(' ERROR: Cannot DCircle_on_FSL3 before DSet_Section_Zoom.')")
CALL DTraceback
END IF
IF (ai_in_path) THEN
WRITE (*,"(' ERROR: Cannot DCircle_on_FSL3 with another path already open.')")
CALL DTraceback
END IF
IF (radius_meters < 0.0D0) THEN
WRITE (*,"(' ERROR: Negative radius to DCircle_on_FSL3: ',1P,E9.2)") radius_meters
CALL DTraceback
END IF
CALL DFSMeters_2_Points (r_meters,z_meters, x_points,y_points)
radius_points = 2834.6D0 * (radius_meters / fs_scale_denominator)
CALL DCircle_on_L12 (2, x_points,y_points, radius_points, stroke, fill)
END SUBROUTINE DCircle_on_FSL3
SUBROUTINE DEnd_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 DProcess_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 DLine_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 DProcess_FSL3_Paths
ELSE !neither stroked nor filled
WRITE (*,"(' ERROR: DEnd_FSL3_path with neither fill nor stroke.')")
CALL DTraceback
END IF
ELSE
WRITE (*,"(' ERROR: Cannot DEnd_FSL3_path when other level open:',I3)") level
CALL DTraceback
END IF
ELSE
WRITE (*,"(' ERROR: Cannot DEnd_FSL3_Path when path not open.')")
CALL DTraceback
END IF
END SUBROUTINE DEnd_FSL3_Path
SUBROUTINE DFSMeters_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*8, INTENT(IN) :: r_meters, z_meters
REAL*8, 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 DFSMeters_2_Points
SUBROUTINE DLine_To_FSL3 (r_meters, z_meters)
IMPLICIT NONE
REAL*8, 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 DTraceback
END IF
ELSE
WRITE (*,"(' ERROR: Cannot DLine_To_FSL3 in path on level ',I2)") &
& ai_current_path_level
CALL DTraceback
END IF
ELSE
WRITE (*,"(' ERROR: Cannot DLine_To_FSL3 before DNew_FSL3_Path.')")
CALL DTraceback
END IF
END SUBROUTINE DLine_To_FSL3
SUBROUTINE DNew_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*8, 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 = 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.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 DNew_FSL3_Path with a path open.')")
CALL DTraceback
END IF
ELSE ! zoom not defined
WRITE (*,"(' ERROR: Must DSet_Zoom before DNew_FSL3_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_FSL3_Path
SUBROUTINE DProcess_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*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 DFSMeters_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 DFSMeters_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 DFSMeters_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 DFSMeters_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.0D0
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_FSL3_Paths
SUBROUTINE DSection_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*8, PARAMETER :: tick_length_points = 4.0D0
REAL*8 :: 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 DSet_Line_Style (width_points = 0.5D0, dashed = .FALSE.)
CALL DSet_Fill_or_Pattern (use_pattern = .FALSE., color_name = 'foreground')
section_length_m = fs_scale_denominator * (ai_window_x2_points - ai_window_x1_points) / 2834.6D0
dr = 1.0D5; dr_km = 100
steps = DInt_Below(section_length_m / dr)
IF (steps == 0) THEN
dr = 1.0D4
dr_km = 10
steps = DInt_Below(section_length_m / dr)
END IF
IF (steps == 0) THEN
dr = 1.0D3
dr_km = 1
steps = DInt_Below(section_length_m / dr)
END IF
y2p = ai_window_y1_points; y1p = y2p - tick_length_points
CALL DBegin_Group
x1p = ai_window_x1_points
CALL DNew_L12_Path(1, x1p, y2p)
CALL DLine_to_L12(x1p, y1p)
CALL DEnd_L12_Path(close = .FALSE., stroke = .TRUE., fill = .FALSE.)
CALL DWrite_L1_Text (x_points = x1p, y_points = y1p, angle_radians = 0.0D0, &
& font_points = font_size_points, lr_fraction = 0.5D0, ud_fraction = 1.0D0, &
& text = '0')
now_at_km = 0
DO i = 1, steps
x1p = x1p + 2834.6D0 * dr / fs_scale_denominator
now_at_km = now_at_km + dr_km
WRITE (c4, "(I4)") now_at_km
CALL DNew_L12_Path(1, x1p, y2p)
CALL DLine_to_L12(x1p, y1p)
CALL DEnd_L12_Path(close = .FALSE., stroke = .TRUE., fill = .FALSE.)
CALL DWrite_L1_Text (x_points = x1p, y_points = y1p, angle_radians = 0.0D0, &
& font_points = font_size_points, lr_fraction = 0.5D0, ud_fraction = 1.0D0, &
& text = TRIM(ADJUSTL(c4)))
END DO
CALL DWrite_L1_Text (x_points = x1p, y_points = y1p, angle_radians = 0.0D0, &
& font_points = font_size_points, lr_fraction = -1.2D0, ud_fraction = 1.0D0, &
& text = 'km')
CALL DEnd_Group
! Tick marks and numbers on left (depth scale)
CALL DSet_Line_Style (width_points = 0.5D0, dashed = .FALSE.)
CALL DSet_Fill_or_Pattern (use_pattern = .FALSE., color_name = 'foreground')
section_depth_m = fs_scale_denominator * (ai_window_y2_points - ai_window_y1_points) / 2834.6D0
dz = 1.0D4; dz_km = 10
steps = DInt_Below(1.0001D0 * section_depth_m / dz)
IF (steps == 0) THEN
dz = 1.0D3
dz_km = 1
steps = DInt_Below(1.0001D0 * section_depth_m / dz)
END IF
x2p = ai_window_x1_points; x1p = x2p - tick_length_points
CALL DBegin_Group
y1p = ai_window_y2_points
CALL DNew_L12_Path(1, x2p, y1p)
CALL DLine_to_L12(x1p, y1p)
CALL DEnd_L12_Path(close = .FALSE., stroke = .TRUE., fill = .FALSE.)
CALL DWrite_L1_Text (x_points = x1p, y_points = y1p, angle_radians = 0.0D0, &
& font_points = font_size_points, lr_fraction = 1.0D0, ud_fraction = 0.4D0, &
& text = '0')
now_at_km = 0
DO i = 1, steps
y1p = y1p - 2834.6D0 * dz / fs_scale_denominator
now_at_km = now_at_km + dz_km
WRITE (c4, "(I4)") now_at_km
CALL DNew_L12_Path(1, x2p, y1p)
CALL DLine_to_L12(x1p, y1p)
CALL DEnd_L12_Path(close = .FALSE., stroke = .TRUE., fill = .FALSE.)
CALL DWrite_L1_Text (x_points = x1p, y_points = y1p, angle_radians = 0.0D0, &
& font_points = font_size_points, lr_fraction = 1.0D0, ud_fraction = 0.4D0, &
& text = TRIM(ADJUSTL(c4)))
END DO
CALL DWrite_L1_Text (x_points = x1p, y_points = y1p, angle_radians = 0.0D0, &
& font_points = font_size_points, lr_fraction = 1.0D0, ud_fraction = 2.0D0, &
& text = 'km')
CALL DEnd_Group
! the bounding frame line:
IF (PRESENT(draw_box).AND.draw_box) THEN
CALL DSet_Line_Style (width_points = 0.5D0, dashed = .FALSE.)
CALL DSet_Stroke_Color ('foreground')
CALL DNew_L12_Path (1, ai_window_x1_points, ai_window_y1_points)
CALL DLine_To_L12 (ai_window_x2_points, ai_window_y1_points)
CALL DLine_To_L12 (ai_window_x2_points, ai_window_y2_points)
CALL DLine_To_L12 (ai_window_x1_points, ai_window_y2_points)
CALL DLine_To_L12 (ai_window_x1_points, ai_window_y1_points)
CALL DEnd_L12_Path (close = .TRUE., stroke = .TRUE., fill = .FALSE.)
END IF ! draw_box
END SUBROUTINE DSection_Km_Frame
SUBROUTINE DSet_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*8, INTENT(IN) :: scale_denominator, &
& section_x1_m, section_y1_m, &
& section_x2_m, section_y2_m, &
& x1_points, y2_points
REAL*8 :: 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.0D0 * 39.370079D0)
!------------------------------------------------------------------
!Matrix for converting (x, y, z) to (r, z, v):
section_length_m = DSQRT((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.0D0
!coefficents of (x, y, & z) in building of z:
fs_xyz_2_rzv_partials(2, 1) = 0.0D0
fs_xyz_2_rzv_partials(2, 2) = 0.0D0
fs_xyz_2_rzv_partials(2, 3) = 1.0D0
!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.0D0
fs_xyz_2_rzv_constants(1) = 0.0D0 - 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.0D0 ! = z when (x = 0, y = 0, z = 0)
fs_xyz_2_rzv_constants(3) = 0.0D0 - 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.0D0
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.0D0
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.0D0
fs_rzv_2_xyz_partials(3, 2) = 1.0D0
fs_rzv_2_xyz_partials(3, 3) = 0.0D0
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.0D0 ! = z when (r = 0, z = 0, v = 0)
!------------------------------------------------------------------------------
!Conversion from r to points-right:
fs_r_2_pointsRight_partial = 2834.6D0 / scale_denominator
fs_r_2_pointsRight_constant = x1_points
!Convertsion from z to points-up:
fs_z_2_pointsUp_partial = 2834.6D0 / scale_denominator
fs_z_2_pointsUp_constant = y2_points
!Conversion from pointsRight to r:
fs_pointsRight_2_r_partial = scale_denominator / 2834.6D0
fs_pointsRight_2_r_constant = 0.0D0 - fs_pointsRight_2_r_partial * x1_points
!Conversion from pointsUp to z:
fs_pointsUp_2_z_partial = scale_denominator / 2834.6D0
fs_pointsUp_2_z_constant = 0.0D0 - fs_pointsUp_2_z_partial * y2_points
END SUBROUTINE DSet_Section_Zoom
SUBROUTINE DStress_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 DFlat_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 DVector_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*8, INTENT(IN) :: s_rr, s_rz, s_zz, s_vv, ref_pressure_SI, ref_diameter_points, r, z
REAL*8 :: from_r, from_z, meters_per_point, radius, s1, s2, to_r, to_z, u1r, u1z, u2r, u2z
meters_per_point = (3.527777D-4) * fs_scale_denominator
CALL DBegin_Group
! Section-perpendicular normal stress (s_vv):
IF (s_vv < 0.0D0) THEN ! circle for compression
radius = (0.5D0*ref_diameter_points)*(-s_vv / ABS(ref_pressure_Si))
radius = radius * meters_per_point
CALL DCircle_on_FSL3 (r, z, radius, .TRUE., .FALSE.)
ELSE IF (s_vv > 0.0D0) THEN ! equilateral triangle for tension
radius = (0.5D0*ref_diameter_points)*(s_vv / ABS(ref_pressure_Si))
radius = radius * meters_per_point
CALL DNew_FSL3_Path (r, z + radius)
CALL DLine_to_FSL3 (r - 0.8660D0 * radius, z - 0.5D0 * radius)
CALL DLine_to_FSL3 (r + 0.8660D0 * radius, z - 0.5D0 * radius)
CALL DLine_to_FSL3 (r, z + radius)
CALL DEnd_FSL3_Path (close = .TRUE., stroke = .TRUE., fill = .FALSE.)
END IF ! circle or triangle for section-perpendicular stress
! Find principal stresses in section plane:
CALL DPrincipal_Axes_22 (s_rr, s_rz, s_zz, &
& s1, s2, u1r,u1z, u2r,u2z)
! Plot sigma1 (in 2-D section plane) axis
radius = (0.5D0*ref_diameter_points)*(ABS(s1) / ABS(ref_pressure_Si))
radius = radius * meters_per_point
IF (s1 < 0.0D0) THEN ! compression: converging arrows
from_r = r + radius * u1r
from_z = z + radius * u1z
CALL DVector_in_Section (from_r, from_z, r, z)
from_r = r - (from_r - r)
from_z = z - (from_z - z)
CALL DVector_in_Section (from_r, from_z, r, z)
ELSE IF (s1 > 0.0D0) THEN ! extension; diverging arrows
to_r = r + radius * u1r
to_z = z + radius * u1z
CALL DVector_in_Section (r, z, to_r, to_z)
to_r = r - (to_r - r)
to_z = z - (to_z - z)
CALL DVector_in_Section (r, z, to_r, to_z)
END IF
! Plot sigma2 (in 2-D section plane) axis
radius = (0.5D0*ref_diameter_points)*(ABS(s2) / ABS(ref_pressure_Si))
radius = radius * meters_per_point
IF (s2 < 0.0D0) THEN ! compression: converging arrows
from_r = r + radius * u2r
from_z = z + radius * u2z
CALL DVector_in_Section (from_r, from_z, r, z)
from_r = r - (from_r - r)
from_z = z - (from_z - z)
CALL DVector_in_Section (from_r, from_z, r, z)
ELSE IF (s2 > 0.0D0) THEN ! extension; diverging arrows
to_r = r + radius * u2r
to_z = z + radius * u2z
CALL DVector_in_Section (r, z, to_r, to_z)
to_r = r - (to_r - r)
to_z = z - (to_z - z)
CALL DVector_in_Section (r, z, to_r, to_z)
END IF
CALL DEnd_Group
END SUBROUTINE DStress_in_Section
SUBROUTINE DTensor_xyz_2_rzv(xyz_tensor, rzv_tensor)
!Makes use of coordinate transformation previously defined by
!DSet_Section_Zoom to convert a symmetric tensor (e.g., stress)
!from (x, y, z) coordinates to (r, z, v) coordinates.
IMPLICIT NONE
REAL*8, DIMENSION(3, 3), INTENT(IN) :: xyz_tensor
REAL*8, 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.0D0 ! 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 DTensor_xyz_2_rzv
SUBROUTINE DVector_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*8, INTENT(IN) :: from_r, from_z, to_r, to_z
REAL*8, PARAMETER :: headlength = 0.1666667D0, headhalfwidth = 0.0666667D0
REAL*8 :: dr, dz, leftr, leftz, rightr, rightz
CALL DBegin_Group
! shaft of arrow
CALL DNew_FSL3_Path (from_r,from_z)
CALL DLine_to_FSL3 (to_r,to_z)
CALL DEnd_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 DNew_FSL3_Path (leftr,leftz)
CALL DLine_to_FSL3 (to_r,to_z)
CALL DLine_to_FSL3 (rightr,rightz)
CALL DEnd_FSL3_Path(close = .FALSE., stroke = .TRUE., fill = .FALSE.)
CALL DEnd_Group
END SUBROUTINE DVector_in_Section
END MODULE DFlat_Section