#!/bin/bash

#
# prepare the smith_konter strain map
#
awk '{print $3}' <  eyy.dat > tmp1
awk '{print $3}' <  exy.dat > tmp2
paste exx.dat tmp1 tmp2 > all.dat
awk '{print $1,$2,sqrt($3*$3+$4*$4+2*$5*$5)*1.e9 }' < all.dat > secinv.gmt
gmt xyz2grd secinv.gmt -Gjunk.grd -I0.02 `gmt grdinfo grad.grd -I0.1/0.1` -NNaN -V
gmt grdsample junk.grd -nb -I.01 `gmt grdinfo grad.grd -I0.1/0.1` -Gmask.grd -V
gmt surface secinv.gmt -I.01 `gmt grdinfo grad.grd -I0.1/0.1` -T.5 -Gsurf.grd -V -Ll0.
gmt grdmath surf.grd mask.grd OR = secinv.grd -V 
gmt grdsample -nb secinv.grd `gmt grdinfo grad.grd -I0.1/0.1` -V -I0.00833333/0.00833333 -F -Gstrainrate.grd
#
#  clean up
#
rm eyy exy all.dat secinv.gmt secinv.grd  surf.grd mask.grd junk.grd
rm tmp1 tmp2 
PROJ="-R-125/-114/31.0/42 -Jm0.65"
gmt gmtset COLOR_MODEL rgb
gmt makecpt -Cno_green -Qo -T1/1000/1 -D -Z > strain.cpt
gmt makecpt -Cgray -T-1.4/.7/.1 -Z > grad.cpt
#
#  mask the deep ocean areas of the strain grid.
#
gmt grdmath topo.grd -2000 GT 0. NAN = mask.grd
cp strainrate.grd tmp.grd
gmt grdmath tmp.grd mask.grd OR = strainratemask.grd
#
#   image the topo and strain grids
#
gmt grdimage grad.grd $PROJ -Cgrad.cpt -Xc -Yc -P -K -V > strain.ps
gmt grdimage strainratemask.grd -Igrad.grd $PROJ -Cstrain.cpt -Q -O -K -V >> strain.ps
gmt pscoast $PROJ -W1 -Di -A200 -Ba2f1g0WSen:."$1": -O -V -N1 -N2 -K >> strain.ps
gmt psxy fault_data.gmt $PROJ -W1 -O -V -K >> strain.ps
gmt psscale -Cstrain.cpt -D1/1.5/1/0.1h -L -O  -K >> strain.ps
gmt pstext $PROJ -O <<!>> strain.ps
-124.2 32 12 0 5 1 2nd invariant (nanostrain/yr)
!
gmt ps2raster strain.ps -A -Tj