#!/bin/sh

#####
#       BOURNE SHELL SCRIPT FOR PLOTTING RESOLVED STRESS COMPONENTS
#####

#####
#       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="normal.out" # STLO STLA STDP STRS

#####
#       GMT PLOTTING VARIABLES
#####
# Map projection (use 'man psbasemap' to see options)
PROJ="-JM4i -P"
# Map limits (-RXMIN/XMAX/YMIN/YMAX)
LIMS="-R31/33/40/42"
# Output PostScript file name
PSFILE="stress.ps"

#####
#       GMT PLOTTING COMMANDS
#####
# Initialize plot
echo 0 0 | gmt psxy $PROJ $LIMS -K > $PSFILE

# Generate color palette for plotting Coulomb stresses
gmt 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
awk '{print $1,$2,$4}' $STS_FILE |\
    gmt xyz2grd -Gstress.grd $LIMS -I0.02/0.02
# Plot stress grid, using color palette generated above
gmt grdimage stress.grd $PROJ $LIMS -Cstress.cpt -K -O >> $PSFILE

# Draw coastline (-W) and national boundaries (-N1)
gmt pscoast $PROJ $LIMS -Dh -W0.75p -N1/0.5p -K -O >> $PSFILE

# Plot focal mechanisms of input faults
awk '{print $1,$2,$3,$4,$5,$6,5}' $FLT_FILE |\
    gmt psmeca $PROJ $LIMS -Sa0.5i -W1p -L1p -Ggrey -K -O >> $PSFILE
# Plot horizontal projection of rectangular input faults
# To convert degrees to radians, multiply by pi/180 = 0.01745
awk '{print $1,$2,$4,$9,$8*cos($5*0.017)}' $FLT_FILE |\
    gmt psxy $PROJ $LIMS -SJ -W3p,darkgreen -K -O >> $PSFILE

# Draw map outline and label axes
gmt psbasemap $PROJ $LIMS -Ba1 -BWeSn -K -O >> $PSFILE

# Finalize plot
echo 0 0 | gmt psxy $PROJ $LIMS -O >> $PSFILE
