#!/bin/tcsh -f

#==========================================================
# gravmap_agemap.gmt - script for plotting color gravity image for a specified region
#
# Carl Tape, 22-Nov-2006
# This script is modified from gravmap.gmt.
# It was further modified by Joann Stock in Nov. 2008 to also plot ages of the
# ocean floor.
#
# 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 cascadia 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
#   gravmap_new.gmt lesser-antilles -65 -54 5.5 21 -Jm0.5 100
#   gravmap_new.gmt sumatra 93 107 -13 7.5 -Jm0.5 100
#
#==========================================================

if ($#argv != 7) then
   echo " "
   echo "Usage: gravmap_agemap.gmt output_name W E S N Jproj output_res"
   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 PLOT_DEGREE_FORMAT D

# input files
set gravimg = "/home/datalib/World_Grav/world_grav.img.15.1"
set ageimg = "/home/datalib/Seafloor_Age/global_age_2.2/age.2.2.img"
set cptfile = "/home/datalib/World_Grav/grav.cpt"
set agecptfile = "age.cpt"
set cptfile1 = $1".cpt"
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"
set agegrd = "$1.age.grd"
set agefile = "$1.age.ps"

# 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 = 300
set cmax = ` grdinfo $1.grd | grep z_min | awk '{print int($5)}'`
set cmin = ` grdinfo $1.grd | grep z_min | awk '{print int($3)}'`
echo "Name of gravity grid file is " $1.grd 
echo " Range of gravity values in grid is " $cmin, $cmax
#grd2cpt -Chaxby -S-${cmax}/${cmax}/20 $grdfile > $cptfile
grd2cpt -Crainbow -S-${cmax}/${cmax}/20 $grdfile > $cptfile1
#makecpt -Crainbow -D -T-200/0/5 -Z > $cptfile

# plot the grd gravity file
grdimage $grdfile $proj $bounds -E${grd_res} -I$gradfile -C$cptfile1 -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
#set plates = "/home/datalib/Plate_Boundaries/Simons_Gravity/simons_grav_sub_zones_LUR"
#psxy $plates $proj $bounds -M -W2p -K -O -P -V >> $psfile

##set plates = "/home/datalib/Plate_Boundaries/nuvel1_boundaries"
#set plates = "/home/datalib/Plate_Boundaries/bird_boundaries"
#psxy $plates $proj $bounds -M -W2p,255/255/255,-- -K -O -P -V >> $psfile

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

display $psfile

# plot the grd agemap file
#
if (! -f $agegrd) then
  echo "wait a moment...making the Age Grid file $agegrd"
  img2grd -V $bounds $ageimg -S0.01 -T1 -G$agegrd -D-80.738/80.738
  grdedit $agegrd $bounds -A
endif
# generate colorpoint file from grdimage
set cmax = 200
#grd2cpt -Crainbow -I -S0/${cmax}/20 $agegrd > $agecptfile
makecpt -Crainbow -I -T0/200/10  > $agecptfile
grdimage $agegrd $proj $bounds -E${grd_res} -C$agecptfile -P -K -V $origin > $agefile
set cinfo = "-W1 -Di -G220"
pscoast $bounds $proj $cinfo -Ba4g4f4/a2g2f2:."Seafloor Age  $1": -K -O -P -V >> $agefile
# colorbar
psscale -C$agecptfile -B50f25:"Seafloor Age, Ma": -D3/-0.5/6/0.25h -P -O -V >> $agefile

display $agefile