#!/bin/sh
#
#PROJ="-R-125/-114/31.0/42.0 -Jm1.7"
PROJ="-R-122.5/-114/31.5/37.5 -Jm2.2"
#PROJ="-R-121.0/-114.0/34.0/39.5 -Jm2.7"
SCALE="-S12i"
gmtset COLOR_MODEL rgb
#makecpt -Cno_green -Qo -T1/1000/1 -D -Z > strain.cpt
makecpt -Cgray -T-1.8/.5/.1 -Z > grad.cpt
#
#  mask the deep ocean areas of the strain grid.
#
grdmath topo.grd -2000 GT 0. NAN = mask.grd
grdmath  $1_exx.grd mask.grd OR = txx.grd
grdmath  $1_exy.grd mask.grd OR = txy.grd
grdmath  $1_eyy.grd mask.grd OR = tyy.grd
#
#   image the topo and strain grids
#
grdimage grad.grd $PROJ -Cgrad.cpt -X1.5 -Y1.5 -K -V >$1_tensor.ps
#
#  resample the grids
#
grdsample txx.grd -Gexx.grd -I80+/100+ 
grdsample txy.grd -Gexy.grd -I80+/100+ 
grdsample tyy.grd -Geyy.grd -I80+/100+
#
#  create the principal stresses and azimuth
#
grdmath 2 exy.grd MUL exx.grd eyy.grd SUB ATAN2 2 DIV = a.grd
grdmath exx.grd eyy.grd SUB 2 POW 0.25 MUL exy.grd 2 POW ADD SQRT exx.grd eyy.grd ADD 2 DIV ADD = e1.grd
grdmath exx.grd eyy.grd SUB 2 POW 0.25 MUL exy.grd 2 POW ADD SQRT NEG exx.grd eyy.grd ADD 2 DIV ADD = e2.grd
grdmath e1.grd a.grd COS MUL = e1_x.grd
grdmath e1.grd a.grd SIN MUL = e1_y.grd
grdmath e2.grd a.grd SIN NEG MUL = e2_x.grd
grdmath e2.grd a.grd COS MUL = e2_y.grd
grdclip e1.grd  -Gmask.grd -Sb0/0 -Sa0/1
grdmath mask.grd e1_x.grd MUL = tx.grd
grdmath mask.grd e1_y.grd MUL = ty.grd
grdvector tx.grd ty.grd $PROJ $SCALE -Q+jc -W1,red -K -O >> $1_tensor.ps
grdclip e1.grd  -Gmask.grd -Sa0/0 -Sb0/1
grdmath mask.grd e1_x.grd MUL = tx.grd
grdmath mask.grd e1_y.grd MUL = ty.grd
grdvector tx.grd ty.grd -R -Jm $SCALE -Q+jc -W1,blue  -K -O >> $1_tensor.ps
grdclip e2.grd  -Gmask.grd -Sb0/0 -Sa0/1
grdmath mask.grd e2_x.grd MUL = tx.grd
grdmath mask.grd e2_y.grd MUL = ty.grd
grdvector tx.grd ty.grd -R -Jm $SCALE -Q+jc -W1,red -K -O >> $1_tensor.ps
grdclip e2.grd  -Gmask.grd -Sa0/0 -Sb0/1
grdmath mask.grd e2_x.grd MUL = tx.grd
grdmath mask.grd e2_y.grd MUL = ty.grd
grdvector tx.grd ty.grd -R -Jm $SCALE -Q+jc -W1,blue  -K -O >> $1_tensor.ps
pscoast $PROJ -W2 -Di -A200 -B2WSen:."$1": -O -V -N1 -N2 -K >> $1_tensor.ps
psxy fault_data.gmt $PROJ -W.5 -O -V -K >> $1_tensor.ps
#psxy ../profiles/profiles.gmt -m $PROJ -W2,white  -O -V -K >> $1_tensor.ps
#psscale -Cstrain.cpt -D5/1.5/8/0.3h -L -O  -K >> $1_tensor.ps
#pstext $PROJ -O <<!>> $1_tensor.ps
#-124.2 32 12 0 5 1 2nd invariant (nanostrain/yr)
#!