#!/bin/bash 

gmt set FONT_ANNOT_PRIMARY 12p,Helvetica COLOR_NAN 255 MAP_FRAME_TYPE fancy  MAP_TICK_LENGTH -0.2 \
MAP_ANNOT_OFFSET_PRIMARY 0.25c MAP_LABEL_OFFSET 0.25c FORMAT_GEO_MAP ddd:mm:ssF

for sta in BEL GHU
do
gmt xyz2grd $sta.DAT -Gmigrs.grd -I0.5/0.5 -R0/450/0/80
gmt grdmath migrs.grd 2.5 MUL = tmp.grd; /bin/mv tmp.grd migrs.grd
gmt grdfilter migrs.grd -D0 -Fg10/10 -Gtmp.grd ; mv tmp.grd migrs.grd
gmt grdsample migrs.grd -Gtmp.grd -I0.1/0.1 ;  mv tmp.grd migrs.grd
mv migrs.grd $sta.grd
done

gmt psbasemap -JX16/-4 -R -B0 -K -X2.5 -Y20 -P > gridplot.ps
gmt psxy -R0/450/0/80 -JX16/-4 -Ba50f25:"Distance (km)":/a20f5:"Depth (km)":WSEn -O -K < /dev/null >> gridplot.ps
gmt makecpt -C../polar.cpt -T-1/1/0.1 -Z -I > test.cpt

gmt grdimage BEL.grd -Ctest.cpt -JX16/-4 -K -O >> gridplot.ps
gmt grdimage GHU.grd -Ctest.cpt -JX16/-4 -K -O >> gridplot.ps


echo 70 -0.04 BEL | gmt pstext -R0/450/0/0.5 -JX16/-3.8 -F+f12p,Helvetica-Bold+jBLM+a90 -O -N -K >> gridplot.ps
echo 387 -0.04 GHU | gmt pstext -R0/450/0/0.5 -JX16/-3.8 -F+f12p,Helvetica-Bold+jBLM+a90 -O -N >> gridplot.ps

rm tmp.input topo.xy gmt.conf gmt.history

gmt psconvert gridplot.ps -A -Tg
gs gridplot.ps