#!/bin/csh
###THis is based on my Pamirs plots, but is designed to do the 
###Altyn Tagh Fault
###This one is designed to overlay the seismicity on the topo
###Hopefully this will do the focalmechs as well

#Here are some variables to save some trouble
set region = '-R80/98/33/41'
set filename = 'atfseissmall'
set projection = '-JM5.8 -P'

###These need to be done only the first time (or whenever you change the 
#range of the grd.

##First make the grd from the gtopo30 data set:
#here is the location:
#http://edcwww.cr.usgs.gov/landdaac/gtopo30/gtopo30.html
#note that you need to be sure that the grdraster.info file is 
#set up correctly.  On alai it is in /disk2/opt/GMT3.1/lib/grdraster.info
#these are the western grids
#echo 'making the southern grd from the raster'
#grdraster 1 -R75/100/30/40 -G/tmp/atfS.grd -I0.5m -V
#echo 'making the northern grd from the raster'
#grdraster 3 -R75/100/40/43 -G/tmp/atfN.grd -I0.5m -V

##Then paste them together:
#echo 'pasting the grds'
#grdpaste /tmp/atfS.grd /tmp/atfN.grd  \
#-G/tmp/atfW.grd -V

##Then make a -S10c -ent (intensity or illumination) file:
#echo 'making the intensity map'
#grdgradient /tmp/atfW.grd -G/tmp/atf.int -A315 -V -Ne0.6


/bin/rm $filename.ps

###Basemap
echo 'plotting basemap'
psbasemap $region  $projection -Ba5  -K > $filename.ps

###Topography
echo 'plotting topo'
grdimage \
/tmp/atfW.grd \
-C/disk2/opt/igmt_1.0/colormaps/topo.cpt \
-I/tmp/atf.int \
-V $region $projection -O -K -U$filename>>$filename.ps

##Coastlines and whatnot
echo 'plotting coastlines and whatnot'
pscoast -R -JM  -Di -N1 -I1/1/0/0/255 \
-S0/0/255 $region $projection -O -K -V \
>> $filename.ps

##historic seismicity
echo 'plotting historic seismicity with stars'
nawk '{print $4, $3, 0.05*($5-4)}' \
histseis.txt | \
psxy $region $projection -Sa -G0/0/0 -W0.2/0/0/0 -L -O -K >>$filename.ps

##paleoseismicity
#these data come from 
#Database of historic Earthquakes reported in the Altun Active Fault Zone
#monograph, page 222.	
echo 'plotting paleoseismic events older than 2500 yrs black'
nawk '{if($1>2500) print $2, $3, 0.05*($4-4)}' \
ssb | \
psxy $region $projection -Sd -G0/0/0 -W0.2/0/0/0 -L -O -K >>$filename.ps
echo 'plotting paleoseismic events younger than 2500 yrs white'
nawk '{if($1<=2500) print $2, $3, 0.05*($4-4)}' \
ssb | \
psxy $region $projection -Sd -G255/255/255 -W0.2/255/255/255 -L -O -K >>$filename.ps


#note that much of the plotting code originally comes
#from Steve Gao (http://stop.at/gao)

echo 'plotting shallow events black'
## Shallow earthquakes
nawk '{if($5 < 50) print $4, $3, 0.05*($6-4)}' \
CNSSevents | \
psxy $region $projection -Sc -G0/0/0 -W0.2/0/0/0 -L -O -K >>$filename.ps

echo 'plotting intermediate events red'
## earthquakes with intermediate depths
nawk '{if($5 >= 50 && $5 <200) print $4, $3, 0.05*($6-4)}' \
CNSSevents | \
psxy $region $projection -Sc -G255/0/0 -W0.2/255/0/0 -L -O -K >>$filename.ps

echo 'plotting deep events red-green'
## Deep earthquakes
nawk '{if($5 >= 200) print $4, $3, 0.05*($6-4)}' \
CNSSevents | \
psxy $region $projection -Sc -G255/255/0 -W0.2/255/255/0 -L -O -K >>$filename.ps

#Magnitude scale
echo 'plotting magnitude scale'

#this is just a mask
psxy << END $region $projection -O -G255 -K -V -A -L >> $filename.ps
80.01 41
83.25 41
83.25 39.3
80.01 39.3
END

psxy << END $region $projection -Sc -G255/0/0 -W0.2/255/0/0 -L -O -K >>$filename.ps
80.25 40.7 0.05
80.25 40.4 0.10
80.25 40.1 0.15
80.25 39.8 0.20
END

psxy << END $region $projection -Sa -G0/0/0 -W0.2/0/0/0 -L -O -K >>$filename.ps
80.25 39.5 0.15
END

pstext << END $region $projection -O -K >>$filename.ps
80.8 40.9 10 0 1 MC M scale
80.6 40.7 10 0 1 MC 5
80.6 40.4 10 0 1 MC 6
80.6 40.1 10 0 1 MC 7
80.6 39.8 10 0 1 MC 8
80.6 39.5 8 0 1 ML Historic earthquake
END

#this is just a mask
psxy << END $region $projection -O -G255 -K -V -A -L >> $filename.ps
98 41
95 41
95 40
98 40
END

psxy << END $region $projection -Sd -G0/0/0 -W0.2/0/0/0 -L -O -K >>$filename.ps
95.35 40.9 0.10
95.35 40.6 0.15
95.35 40.3 0.20
END

pstext << END $region $projection -O -K >>$filename.ps
95.7 40.9 10 0 1 MC 6
95.7 40.6 10 0 1 MC 7
95.7 40.3 10 0 1 MC 8
95.9 40.5 8 0 1 ML White <=2500
95.9 40.3 8 0 1 ML  years BP
END

gs $filename.ps