#!/bin/sh
#
gmt gmtset PS_MEDIA 620x500
#PROJ="-R-123/-114/32.0/37.5 -Jm2.2c"
PROJ="-R-123/-114/32.0/37.5 -JM7i"
SCALE="-S12c"
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.
#
grdsample topo.grd -Rexx_mean.grd -Gtopo.nc
gmt grdgradient topo.nc -A0 -Nt.4 -Ggrad.nc -V
grdmath topo.nc -100 GT 0. NAN = mask.nc
grdmath  exx_mean.grd mask.nc OR = txx.nc
grdmath  exy_mean.grd mask.nc OR = txy.nc
grdmath  eyy_mean.grd mask.nc OR = tyy.nc
#
#   image the topo and strain grids
#
grdimage grad.nc $PROJ -Cgrad.cpt -X1.5 -Y1.5 -K -V >tensor_mean.ps
#
#  resample the grids
#
grdsample txx.nc -Gexx.nc -I50+/60+ 
grdsample txy.nc -Gexy.nc -I50+/60+ 
grdsample tyy.nc -Geyy.nc -I50+/60+
#
#  create the prigrdipal stresses and azimuth
#
grdmath 2 exy.nc MUL exx.nc eyy.nc SUB ATAN2 2 DIV = a.nc
grdmath exx.nc eyy.nc SUB 2 POW 0.25 MUL exy.nc 2 POW ADD SQRT exx.nc eyy.nc ADD 2 DIV ADD = e1.nc
grdmath exx.nc eyy.nc SUB 2 POW 0.25 MUL exy.nc 2 POW ADD SQRT NEG exx.nc eyy.nc ADD 2 DIV ADD = e2.nc
grdmath e1.nc a.nc COS MUL = e1_x.nc
grdmath e1.nc a.nc SIN MUL = e1_y.nc
grdmath e2.nc a.nc SIN NEG MUL = e2_x.nc
grdmath e2.nc a.nc COS MUL = e2_y.nc
grdclip e1.nc  -Gmask.nc -Sb0/0 -Sa0/1
grdmath mask.nc e1_x.nc MUL = tx.nc
grdmath mask.nc e1_y.nc MUL = ty.nc
grdvector tx.nc ty.nc $PROJ $SCALE -Q+jc -W1,red -K -O >> tensor_mean.ps
grdclip e1.nc  -Gmask.nc -Sa0/0 -Sb0/1
grdmath mask.nc e1_x.nc MUL = tx.nc
grdmath mask.nc e1_y.nc MUL = ty.nc
grdvector tx.nc ty.nc -R -JM $SCALE -Q+jc -W1,blue  -K -O >> tensor_mean.ps
grdclip e2.nc  -Gmask.nc -Sb0/0 -Sa0/1
grdmath mask.nc e2_x.nc MUL = tx.nc
grdmath mask.nc e2_y.nc MUL = ty.nc
grdvector tx.nc ty.nc -R -JM $SCALE -Q+jc -W1,red -K -O >> tensor_mean.ps
grdclip e2.nc  -Gmask.nc -Sa0/0 -Sb0/1
grdmath mask.nc e2_x.nc MUL = tx.nc
grdmath mask.nc e2_y.nc MUL = ty.nc
grdvector tx.nc ty.nc -R -JM $SCALE -Q+jc -W1,blue  -K -O >> tensor_mean.ps
pscoast $PROJ -W2 -Di -A200 -B2WSen:."mean": -O -V -N1 -N2 -K >> tensor_mean.ps
psxy fault_data.gmt $PROJ -W.5 -O -V -K >> tensor_mean.ps
#psxy ../profiles/profiles.gmt -m $PROJ -W2,white  -O -V -K >> tensor_mean.ps
#psscale -Cstrain.cpt -D5/1.5/8/0.3h -L -O  -K >> tensor_mean.ps
#pstext $PROJ -O <<!>> tensor_mean.ps
#-124.2 32 12 0 5 1 2nd invariant (nanostrain/yr)
#!
rm *.nc