#!/bin/bash

# edit clean.dat and remove the bad point in San Joaquin Valley by hand
#awk '$1>-119.4 && $1<-119.2 && $2>35.61 && $2<35.63 {print $0}' junk2.dat
#-119.349       35.623  0.8 0.34 -0.83 0.30     5.41244009103   -5.4890783

rm junk*
rm GPS_resid_east.dat GPS_resid_north.dat
awk '{print $1,$2,$3,$4}' < gps.dat > junk0.dat
#awk '{print $1,$2,$3+15.7,$5,$4-17.3,$6 }' < clean.dat > junk0.dat
#
#   change the model to pixel registration
#
grdsample ve0.grd -Gve.grd -r
grdsample vn0.grd -Gvn.grd -r
#
#  extract the model from the grid
#
gmt grdtrack junk0.dat -Gve.grd > junk1.dat
gmt grdtrack junk1.dat -Gvn.grd > junk2.dat

#
# gridding the clean.dat
#
awk '{print ($1,$2,$3-$5)}' < junk2.dat | gmt blockmedian -R-125/-113.333/32/40.336 -I1m -r | awk '{print ($1, $2, $3)}'  >  GPS_resid_east.dat
awk '{print ($1,$2,$4-$6)}' < junk2.dat | gmt blockmedian -R-125/-113.333/32/40.336 -I1m -r | awk '{print ($1, $2, $3)}'  >  GPS_resid_north.dat
#
# combine the east and north residuals
#
gmt convert -A GPS_resid_east.dat GPS_resid_north.dat -o0-2,5 > blk.lluv
#
#  do the gridding
#
gmt gpsgridder `gmt grdinfo -I ve.grd` `gmt grdinfo -I- ve.grd` -Gtmp_%s.nc blk.lluv -r -fg -Vl -Fd8 -Cn250+eigen.txt -S0.35
mv tmp_u.nc ve_resid_surface.grd
mv tmp_v.nc vn_resid_surface.grd
#
#  now grid the residuals and add back to the velocity
#
gmt grdmath ve.grd ve_resid_surface.grd ADD = ve_tot.grd
gmt grdmath vn.grd vn_resid_surface.grd ADD = vn_tot.grd

./plot_Lat_Long.com -125 -114 31.8 40.5 0.65 ve_resid_surface 1
./plot_Lat_Long.com -125 -114 31.8 40.5 0.65 vn_resid_surface 1

./plot_Lat_Long.com -125 -114 31.8 40.5 0.65 vn_tot 5
./plot_Lat_Long.com -125 -114 31.8 40.5 0.65 ve_tot 5

#
#rm GPS_resid_east.dat GPS_resid_north.dat
rm junk*
gmt ps2raster ve_tot_ll.ps -A -Tj
gmt ps2raster vn_tot_ll.ps -A -Tj
gmt ps2raster ve_resid_surface_ll.ps -A -Tj
gmt ps2raster vn_resid_surface_ll.ps -A -Tj