#!/bin/csh
# & rubah menjadi <
gmtset ANNOT_FONT_SIZE_PRIMARY 10p HEADER_FONT_SIZE 18p PLOT_DEGREE_FORMAT ddd:mm:ssF
set region = 90/107/-10/7
set size = M11c
set psfile = 4cross.ps
set grdfile=./sumatra.nc
set cptfile=color.cpt
makecpt -Cglobe -Z > $cptfile
set AZ = 45
#Cross box near subduction
set boxlon = 99.0
set boxlat = -4.0
set boxdep = 0
echo $boxlon $boxlat $boxdep >! box.d
#Draw bathymetry and basemap
grdgradient $grdfile -A30/270 -Gintens.grd -Nt0.30 -V
grdimage $grdfile -R$region -J$size -C$cptfile -Iintens.grd -K -P -Ba4f1NWse -X2 -Y17 >! $psfile
gmtset ANOT_FONT_SIZE 9
pscoast -R$region -J$size -C$cptfile -Dl -W1 -P -Lf92.0/-9/40/60 -Tf92.0/-8.65/10i/3 -O -K >> $psfile
#Draw cross centre
awk -F, '{ print $1, $2}' box.d | psxy -J$size -R$region -Sx.5/.5 -W2::0 -P -O -K >> $psfile
# Draw tectonical structures
awk -F, '{ print $1, $2}' ../focalproposal/trench.gmt | psxy -R$region -J$size -Gdarkblue -Sc0.05 -O -K >> $psfile
awk -F, '{ print $1, $2}' ../focalproposal/transform.gmt | psxy -R$region -J$size -Gdarkblue -Sc0.05 -O -K >> $psfile
awk -F, '{ print $1, $2}' ../focalproposal/ridge.gmt | psxy -R$region -J$size -Gdarkblue -Sc0.05 -O -K >> $psfile
#Draw usgs data
awk -F, '{ print $1, $2}' ../focalproposal/gempa_kecil.lst | psxy -J$size -R$region -Sc0.1c -W0.1 -H -P -O -K >> $psfile
awk -F, '{ print $1, $2}' ../focalproposal/gempa_sedang.lst | psxy -J$size -R$region -Sc0.15c -W0.15 -H -P -O -K >> $psfile
awk -F, '{ print $1, $2}' ../focalproposal/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}' ../focalproposal/trench.gmt | project -Q -C$boxlon/$boxlat -A$AZ -Fxyzpqrs -L-10/550 -W-50/50 >! xyzpqrs_trench
awk -F, '{ print $1, $2, $3}' ../focalproposal/transform.gmt | project -Q -C$boxlon/$boxlat -A$AZ -Fxyzpqrs -L-10/550 -W-50/50 >! xyzpqrs_transform
awk -F, '{ print $1, $2, $3}' ../focalproposal/ridge.gmt | project -Q -C$boxlon/$boxlat -A$AZ -Fxyzpqrs -L-10/550 -W-50/50 >! xyzpqrs_ridge
awk -F, '{ print $1, $2, $3}' ../focalproposal/gempa_kecil.lst | project -Q -C$boxlon/$boxlat -A$AZ -Fxyzpqrs -L-10/550 -W-50/50 >! xyzpqrs_kecil
awk -F, '{ print $1, $2, $3}' ../focalproposal/gempa_sedang.lst | project -Q -C$boxlon/$boxlat -A$AZ -Fxyzpqrs -L-10/550 -W-50/50 >! xyzpqrs_sedang
awk -F, '{ print $1, $2, $3}' ../focalproposal/gempa_besar.lst | project -Q -C$boxlon/$boxlat -A$AZ -Fxyzpqrs -L-10/550 -W-50/50 >! xyzpqrs_besar
#Plot projected earthquakes (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
#Mainshock will appear in front of all earthquakes
#BE CAREFULL !! MAINSHOCK FILE MUST CONSIST OF LON LAT DEPTH
awk '{print $1, $2}' ../focalproposal/cobi.d | psxy -J$size -R$region -Sa0.4c -Gyellow -Wthin -P -O -K >> $psfile
project -Q -C$boxlon/$boxlat -A$AZ ../focalproposal/cobi.d -Fxyzpqrs -L-10/550 -W-50/50 >! xyzpqrs_cobi
#Focalmec lines taken from true epicenter
sed -n '2,3p' ../focalproposal/test.txt | psxy -R$region -J$size -Ggray -W1 -O -K >> $psfile
sed -n '5,6p' ../focalproposal/test.txt | psxy -R$region -J$size -Ggray -W1 -O -K >> $psfile
sed -n '8,9p' ../focalproposal/test.txt | psxy -R$region -J$size -Ggray -W1 -O -K >> $psfile
....................................
#Modified location of Global CMT's focalmec
psmeca ../focalproposal/focaltop_cmt.gmt -R$region -J$size -Sm0.25c/-1 -Gred -W1 -O -K >> $psfile
pstext ../focalproposal/cibi.txt -R$region -J$size -Gblack -O -K >> $psfile
#Drawing box and line of cross section on map
psxy box.txt -J$size -R$region -W0.5 -P -O -K >> $psfile
awk '{print $1, $2}' line.txt | psxy -J$size -R$region -W0.5 -P -O -K >> $psfile
awk '{print $1, $2, 14, 0, 1, "LT", $3 }' line.txt | pstext -J$size -R$region -Gdarkblue -P -O -K >> $psfile
#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 C
520 0 14 0 1 LT C'
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_cobi | 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.04i red 0.010p 0.1i mag 3.5-4.5
S 0.020i c 0.06i red 0.015p 0.1i mag 4.5-5.5
S 0.015i c 0.08i red 0.020p 0.1i mag 5.5-9.5
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
