#
#   script to add horizontal strain rate tensor from one of the models
#   the units of stain rate are nanostrain/yr
#
grdtrack hmax_G05N30_rect_gmt.dat -G../08_average/exx_$1.grd > tmp1.dat
grdtrack tmp1.dat -G../08_average/eyy_$1.grd > tmp2.dat
grdtrack tmp2.dat -G../08_average/exy_$1.grd > tmp3.dat
#
#   compute the azimuth
#
grdmath 2 ../08_average/exy_$1.grd MUL ../08_average/exx_$1.grd ../08_average/eyy_$1.grd SUB ATAN2 2 DIV = a.grd
grdtrack tmp3.dat -Ga.grd > tmp4.dat
#
#   compute the second invariant
#
awk '{print $1,$2,$3,$4,$5,$6,$7,$8,$9,-57.295*$10,sqrt($7*$7+$8*$8+2*$9*$9)}' < tmp4.dat > hmax_G05N30_$1.dat
#
# columns of hmax_G05N30_$1.dat are:
#  lon lat x-scale y-scale azi azi_err exx eyy exy eazi emag
#
#  dump the angles when the strain rate is significant and compute statistics
#
awk '{ if($11 > 10) print $10-$5}' < hmax_G05N30_$1.dat > azi.dat
echo -n $1 "mean, std, " >> stats.dat
awk '{sum+=$1; array[NR]=$1} END {for(x=1;x<=NR;x++){sumsq+=((array[x]-(sum/NR))**2);}print sum/NR, sqrt(sumsq/NR)}' < azi.dat >> stats.dat
#
#  clean up
#
rm tmp* a.grd azi.dat