#!/bin/csh -f
                                                                                                    
# DEFINE DEFAULTS:
set plot_sismi = n
set area = 1
set plot_topo = n
set plot = ""
set palette = etopo5.cpt

# USAGE:
if ($#argv < 2) then
  echo "Usage: run.noam -a area [-s plot_sismi mm mM dm dM] [-t plot_topo] [-o plot]"
  echo " -a  geographic area to plot:"
  echo "     1  --> Central + Eastern US"
  echo "     2  --> Western US"
  echo " -s  plot seismicity: 1 => size = f(mag)"
  echo "                      2 => equal size [$plot_sismi]"
  echo "     then give limits: magmin, magmax, depthmin, depthmax"
  echo " -p  plot topo, y/n [$plot_topo]"
  echo " -o  name of ps file [$plot]"
  exit
endif

# READ USER INPUT
foreach a ($argv)
switch ($a)
  case -o:
    set plot = $argv[2]
    breaksw
  case -t:
    set plot_topo = $argv[2]
    breaksw
  case -s:
    set plot_sismi = $argv[2]
    set mm = $argv[3]; set mM = $argv[4]
    set dm = $argv[5]; set dM = $argv[6]
    breaksw
  case -a:
    set area = $argv[2]
        switch ($area)
        case 1:
          set latm=24; set latM=50; set lonm=-107; set lonM=-65
          set proj=-JM6.5; set grid=-B5a10
          set topo_file = noam_5m.grd
          set sis_file = central_eastern_us.neic
          breaksw
        case 2:
          set latm=30.0; set latM=45.0; set lonm=-128; set lonM=-102
          set proj=-JM6.5; set grid=-B1a5::WNes
          set topo_file = westus.grd
          set sis_file = west_us.neic
          breaksw
        endsw
    breaksw
endsw
shift
end

# DID WE FORGET THE PS FILE NAME?
if ($plot == "") then
  echo -n "Please define ps file name: "
  set plot = $<
  echo "Thank you, now plotting $plot"
endif

# DUMMY GMT COMMAND
set frame="-R$lonm/$lonM/$latm/$latM"
echo $lonm $latm | psxy $proj $frame -Sp -K >! $plot

# PLOT TOPO IMAGE IF REQUESTED
if ($plot_topo == 'y') then
  echo "Plotting topo..."
  grdimage $topo_file $frame $proj -C$palette -O -K >> $plot
endif

# PLOT COASTLINE
pscoast $frame $proj -W1/0 -Di -A400 -N1/1 -N2/0.25tap -K -O >> $plot

# PLOT SEISMICITY
if ($plot_sismi != 'n') then
  echo "Plotting seismicity..."
  awk 'substr($0,50,4)!="    " && NR>18 { \
       lon=substr($0,38,7); lat=substr($0,31,6); mag=substr($0,50,4); dep=substr($0,46,3); \
       lon=strtonum(lon); lat=strtonum(lat); mag=strtonum(mag); dep=strtonum(dep); \
       if (lon>lonmin && lon<lonmax && lat>latmin && lat<latmax && \
       dep>dm && dep<dM && mag>mm && mag<mM) {print lon,lat,mag,mag*mag/300}}' \
       latmin=$latm latmax=$latM lonmin=$lonm lonmax=$lonM \
       dm=$dm dM=$dM mm=$mm mM=$mM $sis_file >! tmp.sis
  if ($plot_sismi == 1) then
    makecpt -T4/9/0.5 -Z -Crainbow >! mag.cpt
    psxy tmp.sis $frame $proj -Sc -W1 -Cmag.cpt -O -K >> $plot
  else if ($plot_sismi == 2) then
    psxy tmp.sis $frame $proj -Sc0.005 -G255/0/0 -O -K >> $plot
  endif
endif

# END MAP
psbasemap $frame $proj $grid -O >> $plot

# ASK USER ABOUT VIEWING
echo -n "Do you want to view $plot? "
set answer = $<
if ($answer == 'y') gv $plot&