#!/bin/sh

#####
#	BOURNE SHELL SCRIPT FOR PLOTTING EFFECTIVE STRESS
#####

#####
#	INPUT/OUTPUT FILES FROM O92UTIL
#####
# Input source fault file
FLT_FILE="fault.dat" # EVLO EVLA EVDP STR DIP RAK SLIP WID LEN
# Output effective stress file
STS_FILE="estress.out" # STLO STLA STDP ESTRS

#####
#       GMT PLOTTING VARIABLES
#####
# Map projection (use 'man psbasemap' to see options)
PROJ="-JM4i -P"
# Map limits (-RXMIN/XMAX/YMIN/YMAX)
LIMS="-R95/97/20/22"
# Output PostScript file name
PSFILE="estress.ps"

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

# Generate color palette for plotting effective stresses
gmt makecpt -Crainbow -T0/1e6/1e5 -D > estress.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,sqrt($4)}' $STS_FILE |\
    gmt xyz2grd -Gestress.grd $LIMS -I0.2/0.2
# Plot effective stress grid, using color palette generated above
gmt grdimage estress.grd $PROJ $LIMS -Cestress.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 -Ba1WeSn -K -O >> $PSFILE

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

