Cross section Sumbagsel

#!/bin/csh
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 <> $psfile
$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 <> $psfile
$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 <> $psfile
$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 <> $psfile
98.0    -2.5
101.5   1.0
EOF
psxy -J$size -R$region -W0.5 -P -O -K <> $psfile
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 <> $psfile
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 <> $psfile
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