#!/bin/tcsh -x

#### gmtset will change the default settings.  type <gmtdefaults -D | more> to see what is possible

gmtset PLOT_DEGREE_FORMAT ddd:mm:ssF
gmtset PAPER_MEDIA Letter MEASURE_UNIT inch

#### set various things like the spatial extent, projections, annotation interval and origin of plot 

set bounds = "-R-122/-114/32/38"
set proj   = "-JM6.5"
set tics   = "-B1neSW"
set origin = "-X1 -Y3"

#### name of target postscript file
set psfile = socal.ps

#### name of file with topographic data (will be made here)
set topogrd   = "socal.grd"

#### name of colormap file that gives you height->color mapping, will be made here
set cpt    = "socal.cpt"

#### name (absolute path) of file with fault locations.  Note file format.
#### datalib is the default location of many useful things in seismolab
set faults = "/home/datalib/Faults/Cal/scitex.lin"
set gpsfile = "socal_vfield_4p0_3D.dat"

#### construct the colormap using a template colormap rescaled for our desired range
makecpt -Csealand -D -T-0/3000/1000 -Z >! $cpt

#### extract the topgraphic data from the relevant database we have already made
grdraster 48 $bounds -G$topogrd -V

#### construct a file of the topographic gradient that will be used to produce the illumination effect
grdgradient $topogrd -A45 -Gintens.grd -Nt -M -V

#### begin the actual plot.  
#### -P means portrait mode.  -K means more stuff will be added later.  -V means verbose output
#### make a default map frame.  Not strictly necessary.  But is a bit cleaner. Will overwrite file.

set misc   = "-K -V -P"
psbasemap $origin $proj $bounds $misc $tics >! $psfile

#### -O means it is following a previous plotting call
#### make the shaded relief image. color for elevation, gradient for illumination.
set misc   = "-K -O -V -P"
grdimage $topogrd $proj $bounds -C$cpt $misc -Iintens.grd  >> $psfile

#### add coastline, color lakes, add national boundaries.
set cinfo  = "-Dh -Cblue -W1pt "
pscoast $proj $bounds $misc $cinfo >> $psfile
set cinfo  = "-N1 -N2 -W1pt -Dh"
pscoast $proj $bounds $misc $cinfo >> $psfile

#### plot the fault database.  -M tells it there are multiple line segments in the file.
psxy $faults $proj $bounds $misc -W1pt -M >> $psfile

#### plot GPS site locations
gawk '{printf($1 " " $2"\n")}' $gpsfile | psxy $proj $bounds $misc -Sc0.05 -Gred -W1 >> $psfile 

#### -P says this is the last command.
#### plot the colorbar
set misc = "-O -V -P"
set scale  = "-D2.5/-0.5/4/0.2h -E -I"
psscale -C$cpt $scale -B500 $misc >> $psfile

#### call something to show the plot.  Could have used <gs> instead
display $psfile &