#!/bin/tcsh -f

#==========================================================
# gravmap_new.gmt - script for plotting color gravity image for a specified region
#
# Carl Tape, 16-Nov-2005
# This script is modified from gravmap.gmt.
#
# INPUT: file label, region, projection, output file resolution
#
# If you want the maximal resolution, then you can either enter
# a large number, or comment out the -E command in grdimage.
#
# If you want to create a NEW grdfile, then be sure to delete the
# old one: the program will only generate a new grdfile if there
# is not one in the local directory.
#
# If you want a ps file, simply comment out the line at the end that
# deletes the ps file.
#
# Examples:
#   gravmap_new.gmt jdf 225 245 35 55 -Jm0.28 100
#   gravmap_new.gmt chile 274 293 -35 -15 -Jm0.3 100
#   gravmap_new.gmt coc-naz-pac 240 270 -10 20 -Jm0.2 100
#
#==========================================================

if ($#argv != 7) then
   echo " "
   echo "Usage: gravmap_new.gmt output_name W E S N Jproj output_res"
   echo "For examples of usage please look at the text of this script"
   echo " "
   exit
endif

set origin = "-X1 -Y2"
set bounds = "-R$2/$3/$4/$5"
set azm = 300
set proj = $6
set grd_res = $7

gmtset HEADER_FONT_SIZE 12 LABEL_FONT_SIZE 10 MEASURE_UNIT inch

# input files
set gravimg = "/home/datalib/World_Grav/world_grav.img.11.1"
set cptfile = "/home/datalib/World_Grav/grav.cpt"
#set cptfile = "grav.cpt"
set plates = "/home/datalib/Plate_Boundaries/nuvel_1_plates"
echo $gravimg
echo $cptfile

# output files
set mercfile = "merc.grd"
set grdfile = "$1.grd"
set gradfile = "$1.grad"
set psfile = "$1.ps"
set jpgfile = "$1.jpg"

# generate the gravity grdfile from the gravity image file
if (! -f $grdfile) then
  echo "wait a moment...making the grd file $grdfile"
  img2grd -S0.1 -V $bounds $gravimg -T1 -m1.0 -G$grdfile
  grdedit $grdfile $bounds -A

  # take gradient of grdfile needed in creating lighting in color grdimage
  grdgradient $grdfile -A$azm -G$gradfile -Nt -V
endif

# generate colorpoint file from grdimage
#set cmax = 100
#grd2cpt -Chaxby -S-${cmax}/${cmax}/20 $grdfile > $cptfile

# plot the grd gravity file
grdimage $grdfile $proj $bounds -E${grd_res} -I$gradfile -C$cptfile -P -K -V $origin > $psfile

# coastlines
set cinfo = "-W1 -Di -G220"
pscoast $bounds $proj $cinfo -Ba4g4f4/a2g2f2:."Gravity $1": -K -O -P -V >> $psfile

# plate boundaries
#psxy $plates $proj $bounds -M -W2p -K -O -P -V >> $psfile

# colorbar
psscale -C$cptfile -B50f25:"Satellite Gravity Anomaly, mGals": -D3/-0.5/6/0.25h -P -O -V >> $psfile

convert $psfile $jpgfile
rm $psfile
xv $jpgfile &