#!/bin/sh

#####
#       BOURNE SHELL SCRIPT FOR PLOTTING CROSS-SECTIONS
#####

#####
#       INPUT/OUTPUT FILES FROM O92UTIL
#####
# Input source fault file
FLT_FILE="fault.dat" # EVLO EVLA EVDP STR DIP RAK SLIP WID LEN
# Output normal/shear/Coulomb stress file
STS_FILE="coul.out" # STLO STLA STDP STRS
# Cross-section file
XSEC_FILE="xsec.dat"

#####
#       GMT PLOTTING VARIABLES
#####
# Map projection (use 'man psbasemap' to see options)
PROJ="-JMx0.02i -P"
# Map limits (-RXMIN/XMAX/YMIN/YMAX)
LIMS="-R-100/100/-100/0"
# Output PostScript file name
PSFILE="xsection.ps"

#####
#       GMT PLOTTING COMMANDS
#####
# Generate color palette for plotting Coulomb stresses
makecpt -Cno_green -T-1e5/1e5/1e4 -D > stress.cpt

# Convert stress output to NetCDF grid file
# -Ixincr/yincr specifies the grid increments, and should be the same
# as the increment used in the grid command
paste $STS_FILE $XSEC_FILE |\
    awk '{print $5,-$6,$4}' |\
    xyz2grd -Gstress.grd $LIMS -I0.02/0.02
# Plot stress grid, using color palette generated above
grdimage stress.grd $PROJ $LIMS -Cstress.cpt -K > $PSFILE

# Draw map outline and label axes
psbasemap $PROJ $LIMS -Ba20:"Distance (km)":/:"Depth (km)":WeSn -O >> $PSFILE

