#!/bin/csh
# Name:    c.surface
# Author:  Marvin Simkin
# Date:    2002-11-13
# Purpose: Generate surfaces from xyz data

# ARGUMENTS (required)
set WESTFEET = $1
set EASTFEET = $2
set SOUTHFEET = $3
set NORTHFEET = $4
set XYZFILE = $5



# calculate X and Y grid spacing
# values 232 and 78 are size of NED data downloaded from seamless.usgs.gov
# this produces a grid that is one pixel too big on both sides, so -1

set XSPACE = `/usr/bin/perl -e "print (($EASTFEET - $WESTFEET) / (232 - 1))"`

set YSPACE = `/usr/bin/perl -e "print (($NORTHFEET - $SOUTHFEET) / (78 - 1))"`

##echo "Calculated X/Y = $XSPACE/$YSPACE vs. 85.16/103.28"

# for geophiz...
source /home/guest/.login
source /home/guest/.cshrc

# ------------------------------------------------------------------------------
# average the known data points when more than one per pixel
echo
set MEAN = $XYZFILE.mean
/usr/bin/rm -f $MEAN
set INPUT = $XYZFILE.xyz
echo /export/apps/gmt/bin/blockmean $INPUT -I$XSPACE/$YSPACE -R$WESTFEET/$EASTFEET/$SOUTHFEET/$NORTHFEET -V '>' $MEAN
/export/apps/gmt/bin/blockmean $INPUT -I$XSPACE/$YSPACE -R$WESTFEET/$EASTFEET/$SOUTHFEET/$NORTHFEET -V > $MEAN


# ------------------------------------------------------------------------------
set GRID = $XYZFILE.grd
set SURF = $XYZFILE.surf
set LINE = $XYZFILE.lin
# Loop through 11 possible tension settings to see what best matches actual
foreach TENSION (0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0)


  # ------------------------------------------------------------------------------
  # cast a surface over the known points
  echo
  /usr/bin/rm -f $TENSION$GRID
  echo /export/apps/gmt/bin/surface $MEAN -I$XSPACE/$YSPACE -R$WESTFEET/$EASTFEET/$SOUTHFEET/$NORTHFEET -V -N1000 -T$TENSION -G$TENSION$GRID
  /export/apps/gmt/bin/surface $MEAN -I$XSPACE/$YSPACE -R$WESTFEET/$EASTFEET/$SOUTHFEET/$NORTHFEET -V -N1000 -T$TENSION -G$TENSION$GRID


  # ------------------------------------------------------------------------------
  # turn that binary .grd file back into ascii so we can work with it
  echo
  /usr/bin/rm -f $TENSION$SURF
  echo /export/apps/gmt/bin/grd2xyz $TENSION$GRID -R$WESTFEET/$EASTFEET/$SOUTHFEET/$NORTHFEET -V '>' $TENSION$SURF
  /export/apps/gmt/bin/grd2xyz $TENSION$GRID -R$WESTFEET/$EASTFEET/$SOUTHFEET/$NORTHFEET -V > $TENSION$SURF


  # ------------------------------------------------------------------------------
  # make a mesh
  /usr/bin/rm -f $TENSION$LINE
  ./surf2lin.pl < $TENSION$SURF > $TENSION$LINE


  # ------------------------------------------------------------------------------
  ## This is sh syntax and lame anyhow:
  ##cut -f2-3 $TENSION$SURF | while read Y Z; do if test "x$Y" = "x$LASTY"; then D="`echo $D` $Z"; else LASTY=$Y; echo $D; D=; fi; done > $TENSION$XYZFILE.d


# ------------------------------------------------------------------------------
# repeat loop
end