#!/bin/csh

/bin/rm tmp.ps tmp.cpt


## 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.1 degree) using grdsample. 
#  The output file is tmph.grd. 

grdsample tmp.grd -Gtmph.grd -I0.1 -Ql/0.1 -Lx

## Get the coordinates of points of 0.1 degree apart along the great circle arc from two points:

sample1d -I0.1 <<END>! profile.xy
-120.0 35.0
-105.0 45.0
END

## Get the data value at the points using grdtrack:

grdtrack profile.xy -Gtmph.grd >! profile.xyz

## Plot the data against logitude (first column):

awk '{print $1, $3}' profile.xyz | \
psxy -R-125/-100/410/440 -JX5/2 -K -X2.0 -Y7.5 -P \
-Ba5f1/a10f1 -W4/255/0/0 >tmp.ps


## For the bottom plot, make a basemap of western United States:
# It shifts by 5.8 inches downward

psbasemap -R-126.000/-100.000/30.500/50.000 \
-JM5 -Ba5f1 -P -O -Y-5.8 -K >>tmp.ps

## Make a color table:

makecpt -Crainbow -T410/450/20 -Z >tmp.cpt

## Produce the image using the resampled grd file:

grdimage tmph.grd -Ctmp.cpt -R -J -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 -Df -Na/1/0t10:5_0 -A1000 -O -K -B  >>tmp.ps

## Plot the profile on the map as dashed line:

awk '{print $1, $2}' profile.xyz | \
psxy -R -J -W20/255/0/0t150:100_0 -O -K >>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