gmtset ANNOT_FONT_SIZE_PRIMARY 10p HEADER_FONT_SIZE 18p PLOT_DEGREE_FORMAT ddd:mm:ssF
set cptfile = color.cpt
set psfile = south_sumatra.ps
set size = M11c
set region = 97/103/-4/3
makecpt -Cglobe -Z > $cptfile
set AZ = 45
#Padang_30/9/2010
set Alon = 98.0
set Alat = -2.5
set Adep = 0
echo $Alon $Alat $Adep >! A.d
set lon_X = 99.67
set lat_X = -0.79
set dep_X = 77.8
set lon_Xplus4 = 101.67
set lat_Xplus0 = -0.79
set lat_Xmin1 = -1.29
echo $lon_X $lat_X $dep_X >! mainshock_X.d
#Draw bathymetry and basemap
grdgradient indo.nc -A30/270 -Gintens.grd -Nt0.30 -V
grdimage indo.nc -R$region -J$size -C$cptfile -Iintens.grd -K -P -Ba1g0.5NWse -X2 -Y10 > $psfile
pscoast -R$region -J$size -P -C$cptfile -Dl -W1 -O -K -Lf108/-6/25/200 >> $psfile
#Draw cross centre
awk -F, '{ print $1, $2}' A.d | psxy -J$size -R$region -Sx.5/.5 -W2::0 -P -O -K >> $psfile
# Draw tectonical structures
awk -F, '{ print $1, $2}' trench.gmt | psxy -R$region -J$size -Gdarkblue -St0.2 -O -K >> $psfile
awk -F, '{ print $1, $2}' transform.gmt | psxy -R$region -J$size -Gdarkblue -St0.2 -O -K >> $psfile
awk -F, '{ print $1, $2}' ridge.gmt | psxy -R$region -J$size -Gdarkblue -St0.2 -O -K >> $psfile
#Draw usgs data
awk -F, '{ print $1, $2}' gempa_kecil.lst | psxy -J$size -R$region -Sc0.1c -W0.1 -H -P -O -K >> $psfile
awk -F, '{ print $1, $2}' gempa_sedang.lst | psxy -J$size -R$region -Sc0.15c -W0.15 -H -P -O -K >> $psfile
awk -F, '{ print $1, $2}' gempa_besar.lst | psxy -J$size -R$region -Sc0.2c -W0.2 -H -P -O -K >> $psfile
#Projection with the center
awk -F, '{ print $1, $2, $3}' trench.gmt | \
project -Q -C$Alon/$Alat -A$AZ -Fxyzpqrs -L-10/550 -W-50/50 >! xyzpqrs_trench
awk -F, '{ print $1, $2, $3}' transform.gmt | \
project -Q -C$Alon/$Alat -A$AZ -Fxyzpqrs -L-10/550 -W-50/50 >! xyzpqrs_transform
awk -F, '{ print $1, $2, $3}' ridge.gmt | \
project -Q -C$Alon/$Alat -A$AZ -Fxyzpqrs -L-10/550 -W-50/50 >! xyzpqrs_ridge
awk -F, '{ print $1, $2, $3}' gempa_kecil.lst | \
project -Q -C$Alon/$Alat -A$AZ -Fxyzpqrs -L-10/550 -W-50/50 >! xyzpqrs_kecil
awk -F, '{ print $1, $2, $3}' gempa_sedang.lst | \
project -Q -C$Alon/$Alat -A$AZ -Fxyzpqrs -L-10/550 -W-50/50 >! xyzpqrs_sedang
awk -F, '{ print $1, $2, $3}' gempa_besar.lst | \
project -Q -C$Alon/$Alat -A$AZ -Fxyzpqrs -L-10/550 -W-50/50 >! xyzpqrs_besar
project -Q -C$Alon/$Alat -A$AZ mainshock_X.d -Fxyzpqrs -L-10/550 -W-50/50 >! xyzpqrs_X
#Plot earthquakes projected (red dots) and mainshock (yellow star) on map
awk '{print $1, $2}' xyzpqrs_kecil | psxy -J$size -R$region -Sc0.1c -W0.1/red -P -O -K >> $psfile
awk '{print $1, $2}' xyzpqrs_sedang | psxy -J$size -R$region -Sc0.15c -W0.15/red -P -O -K >> $psfile
awk '{print $1, $2}' xyzpqrs_besar | psxy -J$size -R$region -Sc0.2c -W0.2/red -P -O -K >> $psfile
awk '{print $1, $2}' mainshock_X.d | psxy -J$size -R$region -Sa0.4c -Gyellow -Wthin -P -O -K >> $psfile
#Focalmec lines for true epicentres
psxy -R$region -J$size -P -O -K -W2 <
$lon_X $lat_X
$lon_Xplus4 $lat_Xplus0
EOF
# Modified focalmec sphere from Global CMT
psmeca -R$region -J$size -Sm0.5c/-1 -Gred -Wthin -O -K <
$lon_Xplus4 $lat_Xplus0 $dep_X 1.76 -0.76 -0.99 0.66 -0.99 -1.94 27 X Y 200909301016A
EOF
pstext -R$region -J$size -Gblack -Wwhite -O -K <
$lon_Xplus4 $lat_Xmin1 9 0 1 LC 2009/09/30
EOF
#Create points of cross section line on map
psxy -J$size -R$region -W0.5 -P -O -K <
98.0 -2.5
101.5 1.0
EOF
psxy -J$size -R$region -W0.5 -P -O -K <
98.318 -2.818
101.813 0.677
101.177 1.313
97.682 -2.182
98.318 -2.818
EOF
pstext -J$size -R$region -Gdarkblue -P -O -K <
98.0 -2.5 14 0 1 LT A
EOF
#DRAW CROSSEC DIAGRAM
psbasemap -JX9.5/-5 -R0/550/0/300 -B50/30NWes -P -O -K -Y-7 >> $psfile
pstext -JX -R -Gdarkblue -P -O -K <
0 0 14 0 1 LT A
EOF
#Plot stars of structures, earthquakes, and mainshocks on crossec diagram
awk '{print $4, $3}' xyzpqrs_trench | \
psxy -JX -R -Gdarkblue -Sc0.2 -P -O -K >> $psfile
awk '{print $4, $3}' xyzpqrs_transform | \
psxy -JX -R -Gdarkblue -Sc0.2 -P -O -K >> $psfile
awk '{print $4, $3}' xyzpqrs_ridge | \
psxy -JX -R -Gdarkblue -Sc0.2 -P -O -K >> $psfile
awk '{print $4, $3}' xyzpqrs_kecil | \
psxy -JX -R -Sc0.1c -W0.1/red -P -O -K >> $psfile
awk '{print $4, $3}' xyzpqrs_sedang | \
psxy -JX -R -Sc0.15c -W0.15/red -P -O -K >> $psfile
awk '{print $4, $3}' xyzpqrs_besar | \
psxy -JX -R -Sc0.2c -W0.2/red -P -O -K >> $psfile
awk '{print $4, $3}' xyzpqrs_X | \
psxy -JX -R -Sa0.4c -Gyellow -Wthin -P -O -K >> $psfile
psscale -C$cptfile -Iintens.grd -D0.4c/1/10/0.15 -B1000 -P -O -K -X11.5c -Y12c >> $psfile
# Create legend of magnitude size
cat << EOF >! magnitude.legend
N 3
S 0.025i c 0.025i red 0.01p 0.050i mag 3-5
S 0.050i c 0.050i red 0.01p 0.075i mag 5-6
S 0.075i c 0.075i red 0.01p 0.100i mag 6-9
>
EOF
pslegend magnitude.legend -Dx0.1i/-0.2i/6.8i/1.5i/TC -Jx0.4i -R0/8/0/8 -O -X-4c -Y-5 >> $psfile
\rm -f *.d *.xy xyzpqrs_*
gs -sDEVICE=x11 $psfile
