#!/bin/csh

/bin/rm tmp.ps tmp.cpt

## Make a basemap of western United States:
psbasemap -E180/45 -R-126.000/-100.000/30.500/50.000 \
-JM5 -Ba5f1 -P -X2.0 -Y2 -K >tmp.ps

## Convert the xyz file into GRD file using xyz2grd
awk '{print $1, $2, $4}' data.dat | xyz2grd -R -I0.5 -Gtmp.grd 

## Resample tmp.grd from coarse (0.5 deg.) into finer resolution (0.2 degree) using grdsample.
#  The output file is tmph.grd.
grdsample tmp.grd -Gtmph.grd -I0.2 -Ql/0.1 -Lx

## Filter tmph.grd using a Gaussian filter of 4 degree wide 
## This is used to get rid of short wavelength features
#  The output file is tmpf.grd. 
grdfilter tmph.grd -D0 -Fg4 -Gtmpf.grd

## Make a color table
makecpt -Crainbow -T410/450/20 -Z >tmp.cpt

## Produce the image using the resampled grd file
grdview tmpf.grd -E180/45 -Ctmp.cpt -R -JM -Jz0.05 -O -K >>tmp.ps


## Plot the scale bar at the bottom of the plot
#psscale -D2.5/-0.4/5.0/0.1h -B10 -Ctmp.cpt -O -K >>tmp.ps


## Overlay with coastlines and political boudnaries
pscoast -R -J -W3 -E180/45 -Df -Na/1/0t10:5_0 -A1000 -O -K -B  >>tmp.ps


## Add a title
pstext -R -J -O -N  <<END>>tmp.ps
-113 51 15 0 0 6 Depth of the 410 km discontinuity: Preliminary results
END

## Convert tmp.ps into a JPG file:

ps2jpeg100 tmp.ps