#!/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 is the big map:
#Here are some variables to save some trouble
set region = '-R74/130/19/45'
set filename = 'atfall'
set projection = '-JM6.8'
set offsets = '-X1 -Y6'

###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 -R70/100/19/40 -G/tmp/atfS.grd -I2.5m -V
echo 'making the northern grd from the raster'
grdraster 3 -R70/100/40/45 -G/tmp/atfN.grd -I2.5m -V
echo 'making the northeast raster'
grdraster 6 -R100/130/40/45 -G/tmp/atfNE.grd -I2.5m -V
echo 'making the southeast raster'
grdraster 5 -R100/130/19/40 -G/tmp/atfSE.grd -I2.5m -V

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

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


/bin/rm $filename.ps

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

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

##Coastlines and whatnot
echo 'plotting coastlines and whatnot'
pscoast $region $projection -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.04*($5-4)}' \
histseis.txt | \
psxy $region $projection -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.04*($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.04*($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.04*($6-4)}' \
CNSSevents | \
psxy $region $projection -Sc -G255/255/0 -W0.2/255/255/0 -L -O -K \
>>$filename.ps

#Magnitude scale
echo 'plotting magnistude scale'

#this is just a mask
psxy << END $region $projection -O -G255 -K -V -A -L >> $filename.ps
110 44.95
129.95 44.95
129.95 40.5
110 40.5
END

psxy << END $region $projection -Sc -G255/0/0 -W0.2/255/0/0 -L -O -K \
>>$filename.ps
127 44 0.04
127 43 0.08
127 42 0.12
127 41 0.16
END

psxy << END $region $projection -Sc -G0/0/0 -W0.2/0/0/0 -L -O -K \
>>$filename.ps
110.5 44.5 0.12
END
psxy << END $region $projection -Sc -G255/0/0 -W0.2/255/0/0 -L -O -K \
>>$filename.ps
110.5 43.5 0.12
END
psxy << END $region $projection -Sc -G255/0/0 -W0.2/255/0/0 -L -O -K \
>>$filename.ps
110.5 42.5 0.12
END
psxy << END $region $projection -Sa -G0/0/0 -W0.2/0/0/0 -L -O -K \
>>$filename.ps
110.5 41.25 0.12
END

pstext << END $region $projection -O -K >>$filename.ps
127.75 44.5 10 0 1 MC M scale
128 44 10 0 1 MC 5
128 43 10 0 1 MC 6
128 42 10 0 1 MC 7
128 41 10 0 1 MC 8
111.5 44.5 10 0 1 ML Depth<50km
111.5 43.5 10 0 1 ML Depth>=50km and <200km
111.5 42.5 10 0 1 ML Depth>=200km
111.5 41.25 10 0 1 ML Historic earthquake
END

echo 'plotting rectangle of the blow up area'
psxy << END $region $projection -O -K -V -A -L >> $filename.ps
80 33
98 33
98 41
80 41
END



#Here is the smaller map
#Here are some variables to save some trouble
set region = '-R80/98/33/41'
set projection = '-JM6.8'
set offsets = '-X0 -Y-5'

###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


###Basemap
echo 'plotting basemap'
psbasemap $region  $projection -Ba5  -K -O -P $offsets>> $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 >>$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


#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 $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

gs $filename.ps