#
# sumatran_quake_magnetic_thickness.gmt
#
#IN1='/home/purucker/MF3/working_model_MF3.xyz'
IN1='working_model_MF3.xyz'
OUT1='sumatran_quake_mag_thickness.ps'
#
gmtset MEASURE_UNIT cm
gmtset BASEMAP_TYPE fancy

#
cat <<END > mag.cpt
#	cpt file created by: grd2cpt -Csealand -Z thick.grd
#COLOR_MODEL = RGB
#
0	140	102	255	4	102	132	255
4	102	132	255	8	102	202	255
8	102	202	255	12	102	255	240
12	102	255	240	14	102	255	171
14	102	255	171	17	102	255	102
17	102	255	102	19	171	255	102
19	171	255	102	22	255	238	166
22	255	238	166	25	255	184	166
25	255	184	166	30	255	181	203
30	255	181	203	40	255	204	243
B	140	102	255
F	255	204	243
N	128	128	128
END
#
/bin/gawk '{print $2,$1,$7}' $IN1 >thick.xyz
surface thick.xyz -Gthick.grd -R75/115/-12/13 -I1/1
grdsample thick.grd -Gthick1.grd -I0.05 -R
grdimage thick1.grd -JM22.3 -Bf5a10WeSn -R -K \
   -Cmag.cpt -Y4 > $OUT1
pscoast -R -JM -Di -O -K -W4 -A1000 >> $OUT1
cat <<END >fault_hypocenter.xy
95.85 3.3
END
psxy fault_hypocenter.xy -R -JM -Sa0.6 -G255/0/0 -O -K >> $OUT1
cat <<END >fault_break1.xy
95.1 2
93 4.5
END
cat <<END >fault_break2.xy
93 4.5
92.5 9
END
psxy fault_break1.xy -R -JM -W7/255/0/0 -O -K  >> $OUT1
psxy fault_break2.xy -R -JM -W7/255/0/0 -O -K  >> $OUT1
#
psscale -D12/-1/15/1h -Cmag.cpt -O -L -B::/:"km": >> $OUT1
#
rm mag.cpt
rm thick.xyz
rm thick.grd
rm thick1.grd
rm fault_break1.xy
rm fault_break2.xy
rm fault_hypocenter.xy
#rm tgrad.grd