#!/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/100/33/41'
set filename = 'atfclose'
set projection = '-JM8'

###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 -R -JM -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  -O -K -V \
>> $filename.ps

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


#our search parameters are:<ul>
#catalog=CNSS
#<li>start_time=1910/01/01,00:00:00
#<li>end_time=1999/12/01,00:00:00
#<li>minimum_latitude=19
#<li>maximum_latitude=45
#<li>minimum_longitude=74
#<li>maximum_longitude=130
#<li>minimum_magnitude=3.0
#<li>maximum_magnitude=9
#<li>minimum_depth=0
#<li>maximum_depth=200
#<li>event_type=E
#
#http://quake.geo.berkeley.edu/cnss/catalog-search.html

#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 -R -JM -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 -R -JM -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 -R -JM -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 -R -JM -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 -R -JM -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 -R -JM -Sa -G0/0/0 -W0.2/0/0/0 -L -O -K >>$filename.ps
80.25 39.5 0.15
END

pstext << END -R -JM -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

gs $filename.ps