#!/bin/bash
#
# plot my example data
#
input_file=ex.data;output_file=ex.ps
scale=50 # scale for earthquake
plot_topo=2 # 2: shade 1: topo 0: coastline
# allow for temporary files
tmpn=/tmp/$USER.$HOST.$$
trap "rm -f $tmpn.* ; exit" 0 1 2 15
# make a colorscale for eq depth
#makecpt -T0/500/50 -Chaxby > $tmpn.cpt
#makecpt -T0/500/50 -D -Cpolar -Z > $tmpn.cpt
makecpt -Qo -T1/300/3 -D -Cpolar -Z > $tmpn.cpt
# topo
#makecpt -Ctopo > mytopo.cpt 
makecpt -Cgray -T-6000/3000/100 > mytopo.cpt 
#colorscale=$tmpn.cpt
colorscale=my.cpt
# automatically get the region
region=`minmax $input_file -I1`
projection=-JM7 # this is my projection
ann="-Ba5f1WesN "
#gmtset BASEMAP_TYPE fancy
# plot coastline
if [ $plot_topo -eq 0 ];then
	pscoast $region $projection -Df -A1000 -Sgray \
		-W0.5 -P -Y2 -K > $output_file
elif [ $plot_topo -eq 1 ];then # plot topo
	# plot topo
	grdimage $HOME/data/etopo2/etopo2.grd -Cmytopo.cpt \
		$region $projection -P -Y2 -K > $output_file
	pscoast $region $projection -Df -A1000 \
		-W0.5  -K -O  >> $output_file
else
	# plot topo and shading
	# cut
	grdcut $HOME/data/etopo2/etopo2.grd \
	    -G$tmpn.grd $region -fg
	# gradient
	grdgradient -Nt $tmpn.grd -G$tmpn.i.grd -A0 
	grdimage -Cmytopo.cpt -I$tmpn.i.grd $tmpn.grd \
		$region $projection -P -Y2 -K > $output_file
	pscoast $region $projection -Df -A1000 \
		-W0.5  -K -O  >> $output_file
fi
# plot plate
psxy $region $projection \
	$HOME/data/plate_boundaries/bird_PB2002/PB2002_tdiddy.gmt \
	-W2,blue -m -O -K -fg >> $output_file
#....
# plot earthquakes
# constant symb size
#gawk '{print($1,$2)}' $input_file | \
#	psxy  -fg $region $projection -Gred -W0.5 \
#            -Sc0.1 -K -O >> $output_file
# adjusting symbol sizes
#gawk -v s=$scale '{print($1,$2,$4/s)}' $input_file | \
#	psxy  -fg $region $projection -Gred -W0.5 \
#            -Sc -K -O >> $output_file
gawk -v s=$scale '{print($1,$2,$3,$4/s)}' $input_file | \
	psxy -fg  $region $projection -C$colorscale \
    -Gred -W0.5 -Sa -K -O >> $output_file
# scale
psscale -C$colorscale -O -K -Ef -D3.5/-.25/3/.2h -L -B:"depth of seismicity":/:"[km]": >> $output_file
#...
# end plot
psbasemap $ann $region $projection -O \
	-Lf9/47.5/`echo $region | gawk -f reg2midlat.awk`/400k+u >> $output_file
modifybb $output_file
gmtset BASEMAP_TYPE plain
echo $0: produced $output_file



