Initial commit

19 files changed, 1650 insertions, 0 deletions

diff --git a/pkg/utilities/nttools/trebin/mkpkg b/pkg/utilities/nttools/trebin/mkpkg
new file mode 100644
index 00000000..7b1e83c3
--- /dev/null
+++ b/pkg/utilities/nttools/trebin/mkpkg
@@ -0,0 +1,27 @@
+# Update the trebin application code in the ttools package library
+# Author: Phil Hodge, 2-DEC-1988
+
+$checkout libpkg.a ../
+$update   libpkg.a
+$checkin  libpkg.a ../
+$exit
+
+libpkg.a:
+	tnamcls.x
+	tnamgio.x
+	tnaminit.x
+	trebin.x		<error.h> <tbset.h>
+	tucspl.f
+	tudcol.x		<tbset.h>
+	tugcol.x		<error.h> <tbset.h>
+	tugetput.x		<tbset.h>
+	tuhunt.f
+	tuiep3.f
+	tuifit.x		trebin.h
+	tuinterp.x		<error.h> <tbset.h>
+	tuiset.x		trebin.h
+	tuispl.f
+	tuival.x		trebin.h
+	tutrim.x
+	tuxget.x		<tbset.h>
+	;
diff --git a/pkg/utilities/nttools/trebin/tnamcls.x b/pkg/utilities/nttools/trebin/tnamcls.x
new file mode 100644
index 00000000..80ded94f
--- /dev/null
+++ b/pkg/utilities/nttools/trebin/tnamcls.x
@@ -0,0 +1,24 @@
+# tnam_cls -- close input & output fnt
+# Close the file name templates.
+#
+# Phil Hodge, 15-Apr-1988  Subroutine created.
+# Phil Hodge,  4-Oct-1995  Modify to use tbn instead of fnt.
+# Phil Hodge, 25-Apr-2000  Add xin_t to calling sequence, and remove dir_only.
+
+procedure tnam_cls (in_t, xin_t, out_t)
+
+pointer in_t		# io: fnt pointer for input tables
+pointer xin_t		# io: fnt pointer for tables of output indep var
+pointer out_t		# io: fnt pointer for output tables
+#--
+
+begin
+	if (in_t != NULL)
+	    call tbnclose (in_t)
+
+	if (xin_t != NULL)
+	    call tbnclose (xin_t)
+
+	if (out_t != NULL)
+	    call tbnclose (out_t)
+end
diff --git a/pkg/utilities/nttools/trebin/tnamgio.x b/pkg/utilities/nttools/trebin/tnamgio.x
new file mode 100644
index 00000000..a3d29af5
--- /dev/null
+++ b/pkg/utilities/nttools/trebin/tnamgio.x
@@ -0,0 +1,79 @@
+# tnam_gio -- get input & output names
+# Get the next input and output table names.
+#
+# Phil Hodge, 15-Apr-1988  Task created.
+# Phil Hodge,  4-Oct-1995  Modify to use tbn instead of fnt; use tbparse.
+# Phil hodge, 16-Apr-1999  Remove ttype from calling sequence of tbparse.
+# Phil Hodge, 25-Apr-2000  Add xin_t, xtable, maxch to calling sequence,
+#			and remove dir_only.
+
+int procedure tnam_gio (in_t, xin_t, out_t, outdir,
+			intable, xtable, outtable, maxch)
+
+pointer in_t			# i: fnt pointer for input tables
+pointer xin_t			# i: fnt pointer for tables of output X
+pointer out_t			# i: fnt pointer for output tables
+char	outdir[ARB]		# i: name of output directory
+char	intable[ARB]		# o: name of next input table
+char	xtable[ARB]		# o: name of next table for output X
+char	outtable[ARB]		# o: name of next output table
+int	maxch			# i: size of table name strings
+#--
+pointer sp
+pointer filename		# name of table file, without brackets
+pointer scratch
+int	nchar			# value of tbnget for intable
+int	junk, hdu, tbparse()
+int	dir_len			# length of root portion of table name
+int	tbnget(), tbnlen(), fnldir()
+errchk	tbparse
+
+begin
+	# Get the next input table name.
+	nchar = tbnget (in_t, intable, SZ_LINE)
+	if (nchar == EOF)
+	    return (EOF)
+
+	# Get the next table for output independent variable values.
+	if (xin_t != NULL) {
+	    if (tbnlen (xin_t) == 1)		# only one xtable?
+		call tbnrew (xin_t)
+	    if (tbnget (xin_t, xtable, SZ_LINE) == EOF)
+		return (EOF)
+	} else {
+	    xtable[1] = EOS
+	}
+
+	if (out_t != NULL) {
+
+	    # Get the next output table name.
+	    if (tbnget (out_t, outtable, SZ_LINE) == EOF)
+		return (EOF)
+
+	} else {		# output is a directory name
+
+	    call smark (sp)
+	    call salloc (filename, SZ_LINE, TY_CHAR)
+	    call salloc (scratch, SZ_LINE, TY_CHAR)
+
+	    # Copy the portion of the table name without brackets to
+	    # Memc[filename]; we need to get rid of the brackets because
+	    # they confuse fnldir.
+	    junk = tbparse (intable, Memc[filename], Memc[scratch],
+			SZ_LINE, hdu)
+
+	    # Get the length of the directory prefix.
+	    dir_len = fnldir (Memc[filename], Memc[scratch], SZ_LINE)
+
+	    # Copy the output directory name to outtable.
+	    call strcpy (outdir, outtable, SZ_LINE)
+
+	    # Append the name of the input file (without directory prefix
+	    # and without any bracket suffix) to the output directory.
+	    call strcat (Memc[filename+dir_len], outtable, SZ_LINE)
+
+	    call sfree (sp)
+	}
+
+	return (nchar)
+end
diff --git a/pkg/utilities/nttools/trebin/tnaminit.x b/pkg/utilities/nttools/trebin/tnaminit.x
new file mode 100644
index 00000000..9e5ef686
--- /dev/null
+++ b/pkg/utilities/nttools/trebin/tnaminit.x
@@ -0,0 +1,75 @@
+include <fset.h>		# to check whether input or output is redirected
+
+# tnam_init -- initialize for input & output names
+# Get the input and output table name lists.  If the output is just a
+# directory name, the name will be copied to outdir; otherwise, the
+# number of names in the input and output lists must be the same.
+#
+# Phil Hodge, 14-Apr-1988  Task created.
+# Phil Hodge, 17-Jun-1993  Change YES to NO in calls to fntopnb.
+# Phil Hodge,  4-Oct-1995  Modify to use tbn instead of fnt.
+# Phil Hodge, 22-Apr-1999  Include explicit test for STDOUT, since
+#		isdirectory thinks STDOUT is a directory.
+# Phil Hodge,  8-Jun-1999  Set input/output to STDIN/STDOUT if redirected.
+# Phil Hodge, 25-Apr-2000  Get inlist, xlist, outlist in trebin, and add
+#		those three and xin_t to the calling sequence.
+
+procedure tnam_init (inlist, xlist, outlist,
+		in_t, xin_t, out_t, outdir, maxch)
+
+char	inlist[ARB]		# i: list of input table names
+char	xlist[ARB]		# i: list of table names for output indep var
+char	outlist[ARB]		# i: list of output table names
+pointer in_t			# o: fnt pointer for input tables
+pointer xin_t			# o: fnt pointer for tables of output X
+pointer out_t			# o: fnt pointer for output tables
+char	outdir[ARB]		# o: if dir_only, name of output directory
+int	maxch			# i: size of outdir string
+#--
+int	n_in, n_xin, n_out	# number of tables in each list
+bool	dir_only		# output just a directory name?
+pointer tbnopen()
+int	isdirectory(), tbnlen()
+bool	strne()
+
+begin
+	dir_only = false
+	if (isdirectory (outlist, outdir, SZ_LINE) > 0 &&
+		strne (outlist, "STDOUT"))
+	    dir_only = true
+
+	in_t = tbnopen (inlist)
+	xin_t = tbnopen (xlist)
+
+	n_in = tbnlen (in_t)
+	n_xin = tbnlen (xin_t)
+	if (n_xin < 1) {
+	    call tbnclose (xin_t)
+	    xin_t = NULL
+	}
+
+	if (dir_only) {
+	    out_t = NULL
+	    n_out = 0
+	} else {
+	    out_t = tbnopen (outlist)
+	    n_out = tbnlen (out_t)
+	}
+
+	if (xin_t != NULL) {
+	    # It's OK to have just one xtable for all intables.
+	    if (n_in != n_xin && n_xin != 1) {
+		call tnam_cls (in_t, xin_t, out_t)
+		call error (1,
+			"intable and xtable lists are not the same length")
+	    }
+	}
+
+	if (out_t != NULL) {
+	    if (n_in != n_out) {
+		call tnam_cls (in_t, xin_t, out_t)
+		call error (1,
+			"intable and outtable lists are not the same length")
+	    }
+	}
+end
diff --git a/pkg/utilities/nttools/trebin/trebin.h b/pkg/utilities/nttools/trebin/trebin.h
new file mode 100644
index 00000000..0164ab3e
--- /dev/null
+++ b/pkg/utilities/nttools/trebin/trebin.h
@@ -0,0 +1,5 @@
+# Codes for interpolating functions.
+define	I_NEAREST	1	# nearest neighbor
+define	I_LINEAR	2	# linear interpolation
+define	I_POLY3		3	# four-point polynomial interpolation
+define	I_SPLINE	4	# cubic spline interpolation
diff --git a/pkg/utilities/nttools/trebin/trebin.x b/pkg/utilities/nttools/trebin/trebin.x
new file mode 100644
index 00000000..873f3b05
--- /dev/null
+++ b/pkg/utilities/nttools/trebin/trebin.x
@@ -0,0 +1,136 @@
+include <error.h>
+include <tbset.h>
+include "trebin.h"
+
+# trebin -- resample to uniform spacing
+# This task resamples a table or list of tables to uniformly spaced
+# values of the independent variable.
+#
+# Phil Hodge, 14-Apr-1988  Task created.
+# Phil Hodge, 30-Jan-1992  Delete inlist, outlist.
+# Phil Hodge, 16-Jun-1993  Set the sign of 'step' based on 'start' and 'end'.
+# Phil Hodge,  4-Oct-1995  Modify to use tbn instead of fnt.
+# Phil Hodge, 21-May-1996  Include extrapolate and ext_value.
+# Phil Hodge, 22-Apr-1999  Get 'step' even if 'start' and 'end' are the same.
+# Phil Hodge, 25-Apr-2000  Get inlist, outlist, xlist in this routine
+#			instead of in tnam_init; also get padvalue.
+# Phil Hodge,  4-Nov-2000  It is an error if step = 0, unless start = end
+
+procedure trebin()
+
+pointer sp
+pointer inlist			# scratch for list of input table names
+pointer outlist			# scratch for list of output table names
+pointer xlist			# scratch for list of table names for X
+pointer intable			# scratch for name of input table
+pointer outtable		# scratch for name of output table
+pointer outdir			# scratch for name of output directory
+pointer xtable			# scratch for name of indep var table
+double	iv_start		# starting value of independent variable
+double	iv_end			# ending value of independent variable
+double	iv_step			# increment in independent variable
+bool	extrapolate		# true means extrapolate if out of bounds
+double	ext_value		# value to use when out of bounds
+double	padvalue		# value at end of input indep. var. to ignore
+char	iv_col[SZ_COLNAME]	# name of independent variable column
+char	func[SZ_FNAME]		# interpolation function
+int	i_func			# interpolation function
+pointer in_t, xin_t, out_t	# fn template pointers for input & output lists
+bool	verbose			# print file names?
+double	clgetd()
+bool	clgetb()
+int	tnam_gio()
+
+begin
+	# Get input and output table template lists.
+	call smark (sp)
+	call salloc (inlist, SZ_LINE, TY_CHAR)
+	call salloc (outlist, SZ_LINE, TY_CHAR)
+	call salloc (xlist, SZ_LINE, TY_CHAR)
+	call salloc (intable, SZ_FNAME, TY_CHAR)
+	call salloc (outtable, SZ_FNAME, TY_CHAR)
+	call salloc (xtable, SZ_FNAME, TY_CHAR)
+	call salloc (outdir, SZ_FNAME, TY_CHAR)
+
+	call clgstr ("intable", Memc[inlist], SZ_LINE)
+	call clgstr ("outtable", Memc[outlist], SZ_LINE)
+	call clgstr ("column", iv_col, SZ_COLNAME)
+	call clgstr ("xtable", Memc[xlist], SZ_LINE)
+
+	# Open the input & output lists of table names.
+	call tnam_init (Memc[inlist], Memc[xlist], Memc[outlist],
+		in_t, xin_t, out_t, Memc[outdir], SZ_FNAME)
+
+	if (xin_t == NULL) {
+
+	    # Get parameters for linearly spaced output independent variable.
+	    iv_start = clgetd ("start")
+	    iv_end = clgetd ("end")
+	    iv_step = clgetd ("step")
+	    if (iv_step == 0.d0 && iv_start != iv_end)
+		call error (1, "step = 0 is invalid")
+
+	    # Set the sign of 'step', rather than expecting the user
+	    # to set it correctly.
+	    if (iv_start < iv_end)
+		iv_step = abs (iv_step)
+	    else if (iv_start > iv_end)
+		iv_step = -abs (iv_step)
+
+	} else {
+
+	    iv_start = 0.d0
+	    iv_end = 0.d0
+	    iv_step = 0.d0
+	}
+
+	call clgstr ("function", func, SZ_FNAME)
+	extrapolate = clgetb ("extrapolate")
+	if (extrapolate)
+	    ext_value = INDEFD			# not used
+	else
+	    ext_value = clgetd ("value")
+
+	padvalue = clgetd ("padvalue")
+
+	verbose = clgetb ("verbose")
+
+	call tuiset (func, i_func)		# set interpolator type
+
+	# Process each table.
+	while (tnam_gio (in_t, xin_t, out_t, Memc[outdir],
+		Memc[intable], Memc[xtable], Memc[outtable], SZ_FNAME) != EOF) {
+
+	    if (verbose) {
+		if (Memc[xtable] != EOS) {
+		    call printf ("%s, %s --> %s\n")
+			call pargstr (Memc[intable])
+			call pargstr (Memc[xtable])
+			call pargstr (Memc[outtable])
+		} else {
+		    call printf ("%s --> %s\n")
+			call pargstr (Memc[intable])
+			call pargstr (Memc[outtable])
+		}
+		call flush (STDOUT)
+	    }
+
+	    iferr {
+		call tuinterp (Memc[intable], Memc[xtable], Memc[outtable],
+			i_func, iv_col, iv_start, iv_end, iv_step,
+			extrapolate, ext_value, padvalue, verbose)
+	    } then {
+		call erract (EA_WARN)
+		if (verbose) {
+		    call eprintf ("This table will be skipped.\n")
+		} else {
+		    call eprintf ("Table %s will be skipped.\n")
+			call pargstr (Memc[intable])
+		}
+		next
+	    }
+	}
+
+	call tnam_cls (in_t, xin_t, out_t)
+	call sfree (sp)
+end
diff --git a/pkg/utilities/nttools/trebin/tucspl.f b/pkg/utilities/nttools/trebin/tucspl.f
new file mode 100644
index 00000000..511211b8
--- /dev/null
+++ b/pkg/utilities/nttools/trebin/tucspl.f
@@ -0,0 +1,52 @@
+	subroutine tucspl (xa, ya, n, work, y2)
+C
+C  Compute the second derivative of YA at each point in the array.
+C  This is the initialization needed in preparation for interpolating
+C  using cubic splines by the subroutine TUISPL.  Input and output
+C  are all double precision.
+C
+C  This routine was copied with slight modifications from the SPLINE
+C  subroutine in Numerical Recipes by Press, Flannery, Teukolsky and
+C  Vetterling.
+C
+C  N		i: number of elements in each array
+C  XA		i: array of independent-variable values
+C  YA		i: array of dependent-variable values
+C  WORK		io: scratch array used for work space
+C  Y2		o: second derivative of YA at each point
+C
+CH  Phil Hodge, 14-Apr-1988  Subroutine copied from Numerical Recipes SPLINE.
+C
+	integer n
+	double precision xa(n), ya(n), work(n), y2(n)
+C--
+	integer i
+	double precision p, sig
+
+C  These values (and y2(n) = 0) are for a "natural" spline.
+	y2(1) = 0.
+	work(1) = 0.
+C
+C  This is the decomposition loop of the tridiagonal algorithm.
+C  Y2 and WORK are used for temporary storage of the decomposed factors.
+C
+	do 10 i = 2, n-1
+	    sig = (xa(i) - xa(i-1)) / (xa(i+1) - xa(i-1))
+	    p = sig * y2(i-1) + 2.
+	    y2(i) = (sig - 1.) / p
+	    work(i) = (6. * ((ya(i+1) - ya(i)) / (xa(i+1) - xa(i))
+     +			- (ya(i) - ya(i-1)) / (xa(i) - xa(i-1))) /
+     +			(xa(i+1) - xa(i-1)) - sig * work(i-1)) / p
+ 10	continue
+
+C					"natural" spline
+	y2(n) = 0.
+C
+C  This is the backsubstitution loop of the tridiagonal algorithm.
+C
+	do 20 i = n-1, 1, -1
+	    y2(i) = y2(i) * y2(i+1) + work(i)
+ 20	continue
+
+	return
+	end
diff --git a/pkg/utilities/nttools/trebin/tudcol.x b/pkg/utilities/nttools/trebin/tudcol.x
new file mode 100644
index 00000000..8a6d1300
--- /dev/null
+++ b/pkg/utilities/nttools/trebin/tudcol.x
@@ -0,0 +1,140 @@
+include <tbset.h>
+
+# tudcol -- get column pointers
+# Get pointers to input and output dependent variable columns, define
+# columns in output table.
+# The arrays icp & ocp of column pointers are of the same length,
+# which will be less than the total number of columns because they
+# do not include the independent variable column.
+#
+# Columns of type text or boolean will not be copied to the output table.
+# If the independent variable column contains arrays, scalar columns will
+# be copied to output without interpolation; array columns must be the same
+# length as the independent variable column.
+# If the independent variable column contains scalars, array columns will
+# not be copied to output.
+# If verbose is true, a message will be printed regarding skipped columns.
+#
+# Phil Hodge, 26-Apr-2000  Subroutine created, based on previous tugcol.
+
+procedure tudcol (itp, otp, iv_colname, outrows,
+		iv_icp, iv_ocp, icp, ocp, ncols, array, verbose)
+
+pointer itp			# i: pointer to input table descriptor
+pointer otp			# i: pointer to output table descriptor
+char	iv_colname[ARB]		# i: name of indep var column
+int	outrows			# i: array length for output array columns
+pointer iv_icp			# o: ptr to descr for input indep var column
+pointer iv_ocp			# o: ptr to descr for output indep var column
+pointer icp[ARB]		# o: ptr to column descr for input table
+pointer ocp[ARB]		# o: ptr to column descr for output table
+int	ncols			# o: number of dependent variable columns
+bool	array			# o: true if indep var column contains arrays
+bool	verbose			# i: print info?
+#--
+pointer sp
+pointer why			# note regarding why a column is skipped
+bool	skip_this		# true if column will not be copied to output
+pointer cp			# a column pointer
+char	cname[SZ_COLNAME]	# a column name
+char	cunits[SZ_COLUNITS]	# units for a column
+char	cfmt[SZ_COLFMT]		# print format for a column
+int	dtype			# data type of a column
+int	xnelem			# number of input elements for indep var column
+int	iv_nelem		# number of output elements for indep var column
+int	nelem			# number of elements
+int	lenfmt			# length of print format
+int	incols			# number of columns in input table
+int	k			# loop index
+int	cnum			# column number (ignored)
+pointer tbcnum()
+int	tbpsta()
+
+begin
+	call smark (sp)
+	call salloc (why, SZ_FNAME, TY_CHAR)
+
+	incols = tbpsta (itp, TBL_NCOLS)
+
+	call tbcfnd1 (itp, iv_colname, iv_icp)
+	if (iv_icp == NULL)
+	    call error (1, "independent variable column not found")
+
+	# Get info about indep var column in input table.
+	call tbcinf (iv_icp, cnum, cname, cunits, cfmt, dtype, xnelem, lenfmt)
+
+	# Note that this test is based on the independent variable column.
+	array = (xnelem > 1)
+
+	# The indep var column in the output table may contain arrays;
+	# iv_nelem will be used in the loop below when defining the output
+	# column of independent variable values.
+	if (array)
+	    iv_nelem = outrows
+	else
+	    iv_nelem = 1
+
+	if (verbose && array) {
+	    call printf ("note:  array columns in each row will be rebinned\n")
+	    call flush (STDOUT)
+	}
+
+	# Define the columns in the output table.
+	ncols = 0
+	do k = 1, incols {
+
+	    skip_this = false		# initial value
+	    cp = tbcnum (itp, k)
+
+	    call tbcinf (cp, cnum, cname, cunits, cfmt, dtype, nelem, lenfmt)
+
+	    # Indep var column.
+	    if (cp == iv_icp) {
+		call tbcdef1 (otp, iv_ocp, cname, cunits, cfmt, dtype, iv_nelem)
+		next
+	    }
+
+	    if (array) {
+		if (nelem > 1 && nelem != xnelem) {	# not the same size
+		    skip_this = true
+		    call strcpy ("array size is not the same",
+				Memc[why], SZ_FNAME)
+		}
+	    } else {
+		if (nelem > 1) {			# skip array columns
+		    skip_this = true
+		    call strcpy ("column contains arrays", Memc[why], SZ_FNAME)
+		}
+		if (dtype <= 0 || dtype == TY_CHAR || dtype == TY_BOOL) {
+		    skip_this = true
+		    if (dtype == TY_BOOL) {
+			call strcpy ("data type is boolean",
+				Memc[why], SZ_FNAME)
+		    } else {
+			call strcpy ("data type is text string",
+				Memc[why], SZ_FNAME)
+		    }
+		}
+	    }
+
+	    if (skip_this) {
+		if (verbose) {
+		    call printf ("  skipping column `%s' (%s)\n")
+			call pargstr (cname)
+			call pargstr (Memc[why])
+		    call flush (STDOUT)
+		}
+		next
+	    }
+
+	    # Define output column; save pointers for input & output.
+	    ncols = ncols + 1
+	    if (array && nelem > 1)
+		nelem = outrows
+	    icp[ncols] = cp
+	    call tbcdef1 (otp, ocp[ncols],
+				cname, cunits, cfmt, dtype, nelem)
+	}
+
+	call sfree (sp)
+end
diff --git a/pkg/utilities/nttools/trebin/tugcol.x b/pkg/utilities/nttools/trebin/tugcol.x
new file mode 100644
index 00000000..6d9dd10d
--- /dev/null
+++ b/pkg/utilities/nttools/trebin/tugcol.x
@@ -0,0 +1,87 @@
+include <error.h>
+include <tbset.h>
+
+# tugcol -- get input X values
+# Get input independent variable column and check it to make
+# sure it is either monotonically increasing or decreasing.
+#
+# Phil Hodge, 18-Apr-1988  Subroutine created
+# Phil Hodge, 30-Jan-1992  Check independent variables more carefully.
+# Phil Hodge, 27-Apr-2000  Move most of this routine to tudcol;
+#			rewrite to allow either array or scalar column.
+
+procedure tugcol (itp, iv_icp, row, xin, xnelem, padvalue, array)
+
+pointer itp			# i: pointer to input table descriptor
+pointer iv_icp			# i: ptr to descr for input indep var column
+int	row			# i: row number, if input column contains arrays
+double	xin[ARB]		# o: input independent variable values
+int	xnelem			# o: actual number of elements in xin array
+double	padvalue		# i: ignore this value at end of xin array
+bool	array			# i: true if input column contains arrays
+#--
+pointer sp
+pointer temp			# scratch for checking indep var for duplicates
+int	nelem			# array size
+int	nvals			# number of elements actually gotten
+int	nrows			# number of rows in input table
+int	i			# loop index
+int	op			# index in temp
+int	tbcigi(), tbpsta(), tbagtd()
+string	NOT_MONOTONIC	"input independent variable is not monotonic"
+
+begin
+	if (array) {
+
+	    nelem = tbcigi (iv_icp, TBL_COL_LENDATA)
+	    nvals = tbagtd (itp, iv_icp, row, xin, 1, nelem)
+	    if (nvals != nelem) {
+		call eprintf (
+	"Not all input independent variable data were gotten from row %d\n")
+		    call pargi (row)
+		call error (1, "")
+	    }
+	    xnelem = nvals
+
+	} else {
+
+	    nrows = tbpsta (itp, TBL_NROWS)
+	    do i = 1, nrows
+		call tbegtd (itp, iv_icp, i, xin[i])
+	    xnelem = nrows
+	}
+
+	# Trim trailing INDEF and pad values by reducing xnelem.
+	call tu_trim (xin, xnelem, padvalue)
+
+	call smark (sp)
+	call salloc (temp, xnelem, TY_DOUBLE)
+
+	# Copy the independent variable data to scratch, skipping embedded
+	# INDEF values.
+	op = 0
+	do i = 1, xnelem {
+	    if (!IS_INDEFD(xin[i])) {
+		Memd[temp+op] = xin[i]	# op is zero indexed at this point
+		op = op + 1
+	    }
+	}
+
+	if (op > 1) {
+	    # Check the independent variable values to make sure they're
+	    # monotonically increasing or decreasing.
+	    if (Memd[temp+1] > Memd[temp]) {		# increasing
+		do i = 2, op {				# one indexed
+		    if (Memd[temp+i-1] <= Memd[temp+i-2])
+			call error (1, NOT_MONOTONIC)
+		}
+	    } else {					# decreasing
+		do i = 2, op {
+		    if (Memd[temp+i-1] >= Memd[temp+i-2])
+			call error (1, NOT_MONOTONIC)
+		}
+	    }
+	}
+
+	call sfree (sp)
+end
diff --git a/pkg/utilities/nttools/trebin/tugetput.x b/pkg/utilities/nttools/trebin/tugetput.x
new file mode 100644
index 00000000..e1359929
--- /dev/null
+++ b/pkg/utilities/nttools/trebin/tugetput.x
@@ -0,0 +1,142 @@
+include <tbset.h>
+
+# tu_getput -- do the interpolation
+# This routine reads the independent and dependent variable values,
+# does the interpolation, and writes the results to the output table.
+#
+# Phil Hodge, 24-April-2000  subroutine created
+
+procedure tu_getput (itp, otp, iv_icp, iv_ocp,
+		icp, ocp, ncols, 
+		xout, outrows, iv_step,
+		i_func, extrapolate, ext_value, padvalue, array, verbose)
+
+pointer itp, otp		# i: descr of input & output tables
+pointer iv_icp			# i: column descr for input indep var column
+pointer iv_ocp			# i: column descr for output indep var column
+pointer icp[ncols]		# i: column descriptors for input table
+pointer ocp[ncols]		# i: column descriptors for output table
+int	ncols			# i: number of columns to copy
+double	xout[ARB]		# i: output indep var values
+int	outrows			# i: size of xout
+double	iv_step			# i: increment in independent variable
+int	i_func			# i: interpolation function code
+bool	extrapolate		# i: true means don't use ext_value
+double	ext_value		# i: value to assign to extrapolated points
+double	padvalue		# i: value at end of input indep. var. to ignore
+bool	array			# i: true if indep var column contains arrays
+bool	verbose			# i: print info?
+#--
+pointer sp
+pointer xin		# scratch for input indep var values
+pointer yin		# scratch for input data values
+pointer y2		# scratch for second derivatives (spline only)
+pointer xa		# scratch for non-indef indep var values
+pointer ya		# scratch for non-indef dep var values
+pointer yout		# array of interpolated values
+double	dbuf		# one interpolated value
+int	inrows		# number of rows in input table
+int	xnelem		# number of elements in xin, INDEF & padvalue trimmed
+int	nelem		# allocated number of elements in array
+int	nvals		# number of array elements actually gotten
+int	n		# size of xa, ya, y2 (counting only non-INDEF values)
+int	row, col	# loop indices
+int	i
+int	tbpsta(), tbcigi(), tbagtd()
+errchk	tugcol, tbagtd, tbaptd, tbegtd, tbeptd, tbrcsc
+
+begin
+	call smark (sp)
+
+	inrows = tbpsta (itp, TBL_NROWS)
+
+	if (array)
+	    nelem = tbcigi (iv_icp, TBL_COL_LENDATA)
+	else
+	    nelem = inrows
+	call salloc (xin, nelem, TY_DOUBLE)
+	call salloc (yin, nelem, TY_DOUBLE)
+	call salloc (y2, nelem, TY_DOUBLE)
+	call salloc (xa, nelem, TY_DOUBLE)
+	call salloc (ya, nelem, TY_DOUBLE)
+
+	if (array) {
+
+	    call salloc (yout, outrows, TY_DOUBLE)
+
+	    # for each row in the input table ...
+	    do row = 1, inrows {
+
+		# Get input independent variable array from current row.
+		call tugcol (itp, iv_icp, row, Memd[xin], xnelem,
+			padvalue, array)
+
+		# Put output indep var values into current row.
+		call tbaptd (otp, iv_ocp, row, xout, 1, outrows)
+
+		# for each column to be interpolated ...
+		do col = 1, ncols {
+		    nelem = tbcigi (icp[col], TBL_COL_LENDATA)
+		    if (nelem == 1) {
+			# just copy scalar column
+			call tbrcsc (itp, otp, icp[col], ocp[col], row, row, 1)
+		    } else {
+			nvals = tbagtd (itp, icp[col], row, Memd[yin], 1, nelem)
+			call tuifit (i_func, Memd[xin], Memd[yin], xnelem,
+				Memd[xa], Memd[ya], Memd[y2], n)
+			if (n > 0) {
+			    do i = 1, outrows {
+				# interpolate
+				call tuival (i_func,
+					Memd[xa], Memd[ya], Memd[y2], n,
+					extrapolate, ext_value, xout,
+					i, outrows, iv_step, Memd[yout+i-1])
+			    }
+			    # put the result
+			    call tbaptd (otp, ocp[col], row,
+					Memd[yout], 1, outrows)
+			} else if (verbose && row == 1) {
+			    call printf ("not enough data for interpolation\n")
+			    call flush (STDOUT)
+			}
+		    }
+		}
+	    }
+
+	} else {
+
+	    # Get input independent variable column.
+	    row = 1
+	    call tugcol (itp, iv_icp, row, Memd[xin], xnelem, padvalue, array)
+
+	    # Put output independent variable values into output table column.
+	    do row = 1, outrows
+		call tbeptd (otp, iv_ocp, row, xout[row])
+
+	    # for each column to be interpolated ...
+	    do col = 1, ncols {
+
+		do row = 1, inrows
+		    call tbegtd (itp, icp[col], row, Memd[yin+row-1])
+
+		# xnelem will be less than or equal to inrows.
+		call tuifit (i_func, Memd[xin], Memd[yin], xnelem,
+			Memd[xa], Memd[ya], Memd[y2], n)
+
+		if (n > 0) {
+		    do row = 1, outrows {
+			# interpolate, and put the result
+			call tuival (i_func, Memd[xa], Memd[ya], Memd[y2], n,
+				extrapolate, ext_value,
+				xout, row, outrows, iv_step, dbuf)
+			call tbeptd (otp, ocp[col], row, dbuf)
+		    }
+		} else if (verbose) {
+		    call printf ("not enough data for interpolation\n")
+		    call flush (STDOUT)
+		}
+	    }
+	}
+
+	call sfree (sp)
+end
diff --git a/pkg/utilities/nttools/trebin/tuhunt.f b/pkg/utilities/nttools/trebin/tuhunt.f
new file mode 100644
index 00000000..3d8df648
--- /dev/null
+++ b/pkg/utilities/nttools/trebin/tuhunt.f
@@ -0,0 +1,103 @@
+	subroutine tuhunt (xa, n, x, klo)
+C
+C  The array XA is searched for an element KLO such that
+C       xa(klo) <= x <= xa(klo+1)
+C  If X < XA(1) then KLO is set to zero; if X > XA(N) then KLO is set to N.
+C  That is, KLO = 0 or N is a flag indicating X is out of bounds.
+C
+C  KLO must be set to an initial guess on input; it will then be replaced
+C  by the correct value on output.
+C
+C  This routine was copied with some modifications from the HUNT
+C  subroutine in Numerical Recipes by Press, Flannery, Teukolsky and
+C  Vetterling.
+C
+C  N		i: number of elements in each array
+C  XA		i: array of independent-variable values
+C  X		i: the value to be bracketed by elements in XA
+C  KLO		io: the lower index in XA that brackets X
+C
+CH  Phil Hodge, 14-Apr-1988  Subroutine copied from Numerical Recipes HUNT.
+CH  Phil Hodge, 21-May-1996  Don't flag endpoints of XA as out of bounds.
+C
+	integer n, klo
+	double precision xa(n), x
+C--
+	integer inc, km, khi
+	logical ascnd
+
+C	Set ASCND, and check for X out of bounds.
+
+	if (xa(n).gt.xa(1)) then
+	    ascnd = .true.
+	    if (x .lt. xa(1)) then
+		klo = 0
+		return
+	    else if (x .gt. xa(n)) then
+		klo = n
+		return
+	    end if
+	else
+	    ascnd = .false.
+	    if (x .gt. xa(1)) then
+		klo = 0
+		return
+	    else if (x .lt. xa(n)) then
+		klo = n
+		return
+	    end if
+	end if
+
+	if ((klo .le. 0) .or. (klo .gt. n)) then
+	    klo = 1
+	    khi = n
+	    go to 3
+	endif
+
+	inc = 1
+	if ((x.ge.xa(klo) .and. ascnd) .or.
+     +      (x.lt.xa(klo) .and. .not.ascnd)) then
+1	    khi = klo + inc
+	    if (khi .gt. n) then
+		khi = n + 1
+	    else if ((x.ge.xa(khi) .and. ascnd) .or.
+     +		     (x.lt.xa(khi) .and. .not.ascnd)) then
+		klo = khi
+		inc = inc + inc
+		go to 1
+	    endif
+	else
+	    khi = klo
+2	    klo = khi - inc
+	    if (klo .lt. 1) then
+		klo = 0
+	    else if ((x.lt.xa(klo) .and. ascnd) .or.
+     +		     (x.ge.xa(klo) .and. .not.ascnd)) then
+		khi = klo
+		inc = inc + inc
+		go to 2
+	    endif
+	endif
+
+3	continue
+C	Before we return, make sure we don't return a value of KLO that
+C	implies X is out of bounds.  We know it isn't because we checked
+C	at the beginning.
+	if (khi-klo .eq. 1) then
+	    klo = max (klo, 1)
+	    klo = min (klo, n-1)
+	    return
+	end if
+
+	km = (khi + klo) / 2
+
+	if ((x .gt. xa(km) .and. ascnd) .or.
+     +      (x .le. xa(km) .and. .not.ascnd)) then
+	    klo = km
+	else
+	    khi = km
+	endif
+
+	go to 3
+
+	end
diff --git a/pkg/utilities/nttools/trebin/tuiep3.f b/pkg/utilities/nttools/trebin/tuiep3.f
new file mode 100644
index 00000000..58f81a99
--- /dev/null
+++ b/pkg/utilities/nttools/trebin/tuiep3.f
@@ -0,0 +1,71 @@
+	subroutine tuiep3 (xa, ya, x, y, dy)
+C
+C  Evaluate a cubic polynomial interpolation function at X.
+C       xa(klo) <= x <= xa(klo+1)
+C  This routine was copied with slight modifications from the POLINT
+C  subroutine in Numerical Recipes by Press, Flannery, Teukolsky and
+C  Vetterling.
+C
+C  XA		i: array of four independent-variable values
+C  YA		i: array of four dependent-variable values
+C  X		i: the independent variable
+C  Y		o: the result of interpolations
+C  DY		o: an estimate of the error of interpolation
+C
+CH  Phil Hodge, 18-Apr-1988  Subroutine copied from Numerical Recipes POLINT.
+C
+	double precision xa(4), ya(4), x, y, dy
+C--
+	integer n
+C					four-point interpolation
+	parameter (n = 4)
+
+	double precision c(n), d(n), dif, dift, den
+	double precision ho, hp, w
+	integer i, m, ns
+
+	do 10 i = 1, n
+	    if (x .eq. xa(i)) then
+		y = ya(i)
+		return
+	    endif
+ 10	continue
+
+	ns = 1
+	dif = abs (x - xa(1))
+
+	do 20 i = 1, n 
+	    dift = abs (x - xa(i))
+	    if (dift .lt. dif) then
+		ns = i
+		dif = dift
+	    endif
+	    c(i) = ya(i)
+	    d(i) = ya(i)
+ 20	continue
+
+	y = ya(ns)
+	ns = ns - 1
+
+	do 40 m = 1, n-1
+
+	    do 30 i = 1, n-m
+		ho = xa(i) - x
+		hp = xa(i+m) - x
+		w = c(i+1) - d(i)
+		den = w / (ho - hp)
+		d(i) = hp * den
+		c(i) = ho * den
+ 30	    continue
+
+	    if (2*ns .lt. n-m) then
+		dy = c(ns+1)
+	    else
+		dy = d(ns)
+		ns = ns - 1
+	    endif
+	    y = y + dy
+ 40	continue
+
+	return
+	end
diff --git a/pkg/utilities/nttools/trebin/tuifit.x b/pkg/utilities/nttools/trebin/tuifit.x
new file mode 100644
index 00000000..2f3a2d84
--- /dev/null
+++ b/pkg/utilities/nttools/trebin/tuifit.x
@@ -0,0 +1,63 @@
+include "trebin.h"
+
+# tuifit -- initialize or fit
+# The input arrays of independent and dependent variable values are copied
+# to output arrays, skipping indef values.  The number of good values
+# (i.e. not indef) is checked to make sure there are enough for the type
+# of interpolation specified.  In the case of spline interpolation, the
+# array y2 of second derivatives is computed.
+#
+# Phil Hodge, 18-Apr-1988  Subroutine created
+# Phil Hodge, 17-Jun-1993  Require at least two values for linear interpolation.
+
+procedure tuifit (i_func, xin, yin, inrows,
+		xa, ya, y2, n)
+
+int	i_func			# i: interpolation function code
+double	xin[ARB]		# i: array of independent-variable values
+double	yin[ARB]		# i: array of dependent-variable values
+int	inrows			# i: size of xin, yin arrays
+double	xa[ARB]			# o: array of independent-variable values
+double	ya[ARB]			# o: array of dependent-variable values
+double	y2[ARB]			# o: used only by spline interpolation
+int	n			# o: size of xa, ya, y2 arrays
+#--
+pointer sp
+pointer work			# scratch for spline fit
+int	k			# loop index
+
+begin
+	n = 0					# initial value
+
+	do k = 1, inrows {
+	    if ( ! IS_INDEFD(xin[k]) && ! IS_INDEFD(yin[k])) {
+		n = n + 1
+		xa[n] = xin[k]
+		ya[n] = yin[k]
+	    }
+	}
+
+	switch (i_func) {
+	case I_NEAREST:
+	    if (n < 1)
+		n = 0				# flag it as bad
+
+	case I_LINEAR:
+	    if (n < 2)
+		n = 0
+
+	case I_POLY3:
+	    if (n < 4)
+		n = 0
+
+	case I_SPLINE:
+	    if (n >= 4) {
+		call smark (sp)
+		call salloc (work, n, TY_DOUBLE)
+		call tucspl (xa, ya, n, Memd[work], y2)
+		call sfree (sp)
+	    } else {
+		n = 0
+	    }
+	}
+end
diff --git a/pkg/utilities/nttools/trebin/tuinterp.x b/pkg/utilities/nttools/trebin/tuinterp.x
new file mode 100644
index 00000000..3a3b1d89
--- /dev/null
+++ b/pkg/utilities/nttools/trebin/tuinterp.x
@@ -0,0 +1,139 @@
+include <error.h>
+include <tbset.h>
+include "trebin.h"
+
+# tuinterp -- interpolate to regrid a table
+# Open the input & output tables, interpolate to uniformly spaced values
+# of the independent variable, and close the tables.
+#
+# Phil Hodge, 15-Apr-1988  Subroutine created
+# Phil Hodge, 12-May-1989  Include check for not enough data to interpolate.
+# Phil Hodge, 12-Jun-1989  Also copy header parameters.
+# Phil Hodge, 30-Jan-1992  Call tbtclo instead of close.
+# Phil Hodge, 16-Jun-1993  Check number of rows for interpolation function.
+# Phil Hodge,  4-Apr-1994  Errchk tbtopn, and use iferr for tbtcre.
+# Phil Hodge, 20-May-1996  Pass extrapolate and ext_value to tuival.
+# Phil Hodge, 29-Jul-1998  Add iv_step to calling sequence of tuival.
+# Phil Hodge,  8-Apr-1999  Call tbfpri.
+# Phil Hodge, 22-Apr-1999  Don't set output table type if outtable = STDOUT.
+# Phil Hodge, 25-Apr-2000  Add xtable, padvalue, verbose to calling sequence.
+# Phil Hodge, 30-Oct-2001  Delete just the output table, not the whole file,
+#			if there's an error.
+# Phil Hodge,  2-Jan-2002  Remove the statements to delete the output table
+#			in case the input table is not monotonic (because
+#			calling tbtclo for a text table caused the error
+#			message to be replaced by a misleading message).
+
+procedure tuinterp (intable, xtable, outtable,
+		i_func, iv_colname, iv_start, iv_end, iv_step,
+		extrapolate, ext_value, padvalue, verbose)
+
+char	intable[ARB]		# i: name of input table
+char	xtable[ARB]		# i: table of output indep var values
+char	outtable[ARB]		# i: name of output table
+int	i_func			# i: interpolation function code
+char	iv_colname[SZ_COLNAME]	# i: name of independent variable column
+double	iv_start		# i: starting value of independent variable
+double	iv_end			# i: ending value of independent variable
+double	iv_step			# i: increment in independent variable
+bool	extrapolate		# i: true means don't use ext_value
+double	ext_value		# i: value to assign to extrapolated points
+double	padvalue		# i: value at end of input indep. var. to ignore
+bool	verbose			# i: print info?
+#--
+pointer sp
+pointer itp, otp		# descr of input & output tables
+pointer iv_icp			# descr for input indep var column
+pointer iv_ocp			# descr for output indep var column
+pointer icpp, ocpp		# column descr for i & o tables
+pointer xout			# scratch for output indep var values
+int	ttype			# indicates row- or column-ordered
+int	ncols			# number of dependent variable columns
+int	incols			# total number of input columns
+int	outrows			# number of rows in output table
+int	phu_copied		# set by tbfpri and ignored
+bool	array			# true if indep var column contains arrays
+pointer tbtopn()
+int	tbpsta()
+bool	strne()
+errchk	tbfpri, tbtopn, tbtdel, tudcol, tuxget, tu_getput
+
+begin
+	call smark (sp)
+
+	itp = tbtopn (intable, READ_ONLY, 0)
+	incols = tbpsta (itp, TBL_NCOLS)
+
+	# Either read or compute the values of the independent variable
+	# at which we will interpolate the values from the input table.
+	iferr {
+	    call tuxget (xtable, iv_start, iv_end, iv_step, padvalue,
+		xout, outrows)
+	} then {
+	    call tbtclo (itp)
+	    call erract (EA_ERROR)
+	}
+
+	iferr {
+	    call tbfpri (intable, outtable, phu_copied)
+	    otp = tbtopn (outtable, NEW_FILE, NULL)
+	} then {
+	    call mfree (xout, TY_DOUBLE)
+	    call tbtclo (itp)
+	    call erract (EA_ERROR)
+	}
+
+	call salloc (icpp, incols, TY_POINTER)
+	call salloc (ocpp, incols, TY_POINTER)
+
+	# Define output columns, and get column pointers.
+	iferr {
+	    call tudcol (itp, otp, iv_colname, outrows,
+		iv_icp, iv_ocp, Memi[icpp], Memi[ocpp], ncols, array, verbose)
+	} then {
+	    call mfree (xout, TY_DOUBLE)
+	    call tbtclo (itp)
+	    call tbtclo (otp)
+	    call erract (EA_ERROR)
+	}
+
+	# Output table should be same type as input table.
+	if (strne (outtable, "STDOUT")) {
+	    ttype = tbpsta (itp, TBL_WHTYPE)
+	    call tbpset (otp, TBL_WHTYPE, ttype)
+	    if (ttype == TBL_TYPE_S_COL) {
+		call tbpset (otp, TBL_ALLROWS,
+			max (outrows, tbpsta (itp, TBL_ALLROWS)))
+	    }
+	}
+
+	iferr {
+	    call tbtcre (otp)			# create output table
+	} then {
+	    call mfree (xout, TY_DOUBLE)
+	    call tbtclo (itp)
+	    call tbtclo (otp)
+	    call erract (EA_ERROR)
+	}
+
+	call tbhcal (itp, otp)			# copy all header parameters
+
+	# For each column, get the data, do the interpolation,
+	# write the results.
+	iferr {
+	    call tu_getput (itp, otp, iv_icp, iv_ocp,
+		Memi[icpp], Memi[ocpp], ncols, 
+		Memd[xout], outrows, iv_step,
+		i_func, extrapolate, ext_value, padvalue, array, verbose)
+	} then {
+	    call mfree (xout, TY_DOUBLE)
+	    call tbtclo (itp)
+	    call erract (EA_ERROR)
+	}
+
+	call mfree (xout, TY_DOUBLE)
+	call sfree (sp)
+
+	call tbtclo (itp)
+	call tbtclo (otp)
+end
diff --git a/pkg/utilities/nttools/trebin/tuiset.x b/pkg/utilities/nttools/trebin/tuiset.x
new file mode 100644
index 00000000..427df021
--- /dev/null
+++ b/pkg/utilities/nttools/trebin/tuiset.x
@@ -0,0 +1,26 @@
+include "trebin.h"
+
+# tuiset -- set interpolation type
+# 
+#
+# P.E. Hodge, 18-Apr-88  Subroutine created
+
+procedure tuiset (func, i_func)
+
+char	func[ARB]		# i: interpolation function
+int	i_func			# o: interpolation function code
+#--
+int	strncmp()
+
+begin
+	if (func[1] == 'n')
+	    i_func = I_NEAREST
+	else if (func[1] == 'l')
+	    i_func = I_LINEAR
+	else if (strncmp (func, "poly3", 5) == 0)
+	    i_func = I_POLY3
+	else if (func[1] == 's')
+	    i_func = I_SPLINE
+	else
+	    call error (1, "unknown interpolation function")
+end
diff --git a/pkg/utilities/nttools/trebin/tuispl.f b/pkg/utilities/nttools/trebin/tuispl.f
new file mode 100644
index 00000000..2f4c68b8
--- /dev/null
+++ b/pkg/utilities/nttools/trebin/tuispl.f
@@ -0,0 +1,32 @@
+	subroutine tuispl (xa, ya, y2, x, y)
+C
+C  Interpolate at point X using cubic splines.  The array Y2 must have
+C  previously been computed by calling TUCSPL.  Note that XA, YA, Y2
+C  are two-element subarrays of the arrays with the same names elsewhere.
+C  Input and output are all double precision.
+C  This routine was copied with slight modifications from the SPLINT
+C  subroutine in Numerical Recipes by Press, Flannery, Teukolsky and
+C  Vetterling.
+C
+C  XA		i: pair of independent-variable values
+C  YA		i: pair of dependent-variable values
+C  Y2		i: second derivatives of YA at each point
+C  X		i: value at which spline is to be computed
+C  Y		o: interpolated value at X
+C
+CH  Phil Hodge, 14-Apr-1988  Subroutine copied from Numerical Recipes SPLINT.
+C
+	double precision xa(2), ya(2), y2(2), x, y
+C--
+	double precision h, a, b
+
+	h = xa(2) - xa(1)
+
+	a = (xa(2) - x) / h
+	b = (x - xa(1)) / h
+	y = a * ya(1) + b * ya(2) +
+     +		((a**3 - a) * y2(1) + (b**3 - b) * y2(2))
+     +		* h * h / 6.
+
+	return
+	end
diff --git a/pkg/utilities/nttools/trebin/tuival.x b/pkg/utilities/nttools/trebin/tuival.x
new file mode 100644
index 00000000..20354a5d
--- /dev/null
+++ b/pkg/utilities/nttools/trebin/tuival.x
@@ -0,0 +1,272 @@
+include "trebin.h"
+
+# tuival -- interpolate
+# Interpolate to get one value.
+# A zero value of iv_step means that the output independent-variable values
+# may not be uniformly spaced, so the interval x1,x2 must be gotten from
+# the xout array.  (This is relevant for linear interpolation only.)
+#
+# Phil Hodge, 18-Apr-1988  Subroutine created
+# Phil Hodge, 20-May-1996  Include extrapolate and ext_value for extrapolation.
+# Phil Hodge, 29-Jul-1998  For 1-D, fit a line instead of interpolating,
+#			add iv_step to calling sequence, add subroutines
+#			tu_range and tu_fit1.
+# Phil Hodge, 22-Apr-1999  The test on whether x is outside the range of
+#		the data did not consider that xa could be decreasing.
+# Phil Hodge, 25-Apr-2000  Change calling sequence:  x --> xout, j, outrows;
+#		this is to support the option of a nonuniformly spaced output
+#		independent variable array.
+
+procedure tuival (i_func, xa, ya, y2, n,
+	extrapolate, ext_value, xout, j, outrows, iv_step, y)
+
+int	i_func			# i: interpolation function code
+double	xa[n]			# i: input independent-variable values
+double	ya[n]			# i: input dependent-variable values
+double	y2[n]			# i: used only by spline interpolation
+int	n			# i: size of xa, ya, y2 arrays
+bool	extrapolate		# i: extrapolate rather than use ext_value?
+double	ext_value		# i: value to assign if x is out of bounds
+double	xout[ARB]		# i: output independent variable array
+int	j			# i: current element of xout
+int	outrows			# i: size of xout array
+double	iv_step			# i: spacing of x values (can be negative)
+double	y			# o: the interpolated value
+#--
+double	x			# independent variable
+double	x1, x2			# beginning and end of output pixel
+double	dy			# an estimate of the error of interpolation
+int	klo, khi		# indexes within xa that bracket x
+int	npts			# number of points in interval
+
+begin
+	x = xout[j]
+
+	# Initial guess for klo; we want xa[klo] <= x <= xa[klo+1].
+	klo = (n + 1) / 2
+	khi = klo
+
+	# If x falls outside the range of the data, either extrapolate
+	# (below) or assign a fixed value ext_value.
+	if (!extrapolate && n > 1) {
+	    if (xa[1] < xa[2] && (x < xa[1] || x > xa[n])) {	# increasing
+		y = ext_value
+		return
+	    }
+	    if (xa[1] > xa[2] && (x > xa[1] || x < xa[n])) {	# decreasing
+		y = ext_value
+		return
+	    }
+	}
+
+	switch (i_func) {
+	case I_NEAREST:
+
+	    # Find klo such that xa[klo] <= x <= xa[klo+1], starting from
+	    # current klo.
+	    call tuhunt (xa, n, x, klo)
+	    klo = max (klo, 1)
+	    klo = min (klo, n)
+	    khi = min (klo+1, n)
+	    if (abs (x - xa[klo]) < abs (x - xa[khi]))
+		y = ya[klo]
+	    else
+		y = ya[khi]
+
+	case I_LINEAR:
+
+	    # Get x1 and x2.
+	    call tu_x1x2 (xout, j, outrows, iv_step, x1, x2)
+
+	    # Get klo, khi, and npts.
+	    call tu_range (xa, ya, n, x1, x2, klo, khi, npts)
+	    if (npts > 1) {
+		# Fit a line to the data from klo to khi inclusive.
+		call tu_fit1 (xa, ya, x, klo, khi, y)
+	    } else {
+		# Interpolate.
+		klo = klo - 1
+		call tuhunt (xa, n, x, klo)
+		klo = max (klo, 1)
+		klo = min (klo, n-1)
+		call tuiep1 (xa[klo], ya[klo], x, y)
+	    }
+
+	case I_POLY3:
+
+	    call tuhunt (xa, n, x, klo)
+	    klo = max (klo, 2)
+	    klo = min (klo, n-2)
+	    # Pass tuiep3 the four points at klo-1, klo, klo+1, klo+2.
+	    call tuiep3 (xa[klo-1], ya[klo-1], x, y, dy)
+
+	case I_SPLINE:
+
+	    call tuhunt (xa, n, x, klo)
+	    klo = max (klo, 1)
+	    klo = min (klo, n-1)
+	    call tuispl (xa[klo], ya[klo], y2[klo], x, y)
+	}
+end
+
+# This routine gets the endpoints x1,x2 of the current output pixel.
+# If the output independent variable array is uniformly spaced, x1 and x2
+# will just be the current pixel xout[j] - or + iv_step/2.  If the output
+# independent variable array is (or may be) nonuniformly spaced, indicated
+# by iv_step being zero, then x1 and x2 will be the midpoints of intervals
+# adjacent to xout[j] on the left and right.
+#
+# x1 may be less than, equal to, or greater than x2.
+
+procedure tu_x1x2 (xout, j, outrows, iv_step, x1, x2)
+
+double	xout[ARB]	# i: output independent variable array
+int	j		# i: current element of xout
+int	outrows		# i: size of xout array
+double	iv_step		# i: spacing of x values (can be negative)
+double	x1, x2		# o: endpoints of output pixel
+
+begin
+	if (iv_step != 0) {
+	    # output is uniformly spaced
+	    x1 = xout[j] - iv_step / 2.d0
+	    x2 = xout[j] + iv_step / 2.d0
+	} else {
+	    if (outrows == 1) {
+		x1 = xout[1]
+		x2 = xout[1]
+	    } else if (j == 1) {
+		x2 = (xout[1] + xout[2]) / 2.d0
+		x1 = xout[1] - (x2 - xout[1])
+	    } else if (j == outrows) {
+		x1 = (xout[outrows-1] + xout[outrows]) / 2.d0
+		x2 = xout[outrows] + (xout[outrows] - x1)
+	    } else {
+		x1 = (xout[j-1] + xout[j]) / 2.d0
+		x2 = (xout[j] + xout[j+1]) / 2.d0
+	    }
+	}
+end
+
+# tu_range -- find range of indexes
+# This routine finds the range of indexes in xa and ya corresponding to
+# the interval x1 to x2.  The range is klo to khi; klo will be less than
+# or equal to khi, even though xa can be either increasing or decreasing.
+# klo and khi will also be within the range 1 to n inclusive, even if
+# x1 to x2 is outside the range xa[1] to xa[n].  npts will be set to
+# the actual number of elements of xa within the interval x1 to x2;
+# thus, npts can be zero.
+#
+# x1 and x2 will be interchanged, if necessary, so that the array index
+# in xa corresponding to x1 will be smaller than the array index corresponding
+# to x2, i.e. so klo will be smaller than khi.
+#
+# Note that klo is used as an initial value when hunting for a value in xa,
+# so klo must have been initialized to a reasonable value before calling
+# this routine.
+
+procedure tu_range (xa, ya, n, x1, x2, klo, khi, npts)
+
+double	xa[n]		# i: array of independent-variable values
+double	ya[n]		# i: array of dependent-variable values
+int	n		# i: size of xa, ya, y2 arrays
+double	x1, x2		# io: endpoints of output pixel
+int	klo, khi	# io: range of indices in xa within x1,x2
+int	npts		# o: number of elements in xa within x1,x2
+#--
+double	temp		# for swapping x1 and x2, if appropriate
+int	k1, k2
+
+begin
+	# Swap x1 and x2 if necessary, so that klo will be less than or
+	# equal to khi.
+	if (xa[1] <= xa[2]) {	# input X is increasing
+	    if (x1 > x2) {
+		temp = x1; x1 = x2; x2 = temp
+	    }
+	} else {		# input X is decreasing
+	    if (x1 < x2) {
+		temp = x1; x1 = x2; x2 = temp
+	    }
+	}
+
+	k1 = klo
+	call tuhunt (xa, n, x1, k1)
+	k1 = k1 + 1			# next point
+	klo = k1
+	klo = min (klo, n)
+
+	k2 = k1
+	call tuhunt (xa, n, x2, k2)
+	khi = max (k2, 1)
+	khi = min (khi, n)
+
+	# npts can be different from khi - klo + 1, and it can be zero.
+	npts = k2 - k1 + 1
+end
+
+# tu_fit1 -- fit a line (1-D) to data
+# This routine fits a straight line to (xa[k],ya[k]) for k = klo to khi
+# inclusive.  The fit is then evaluated at x, and the result is returned
+# as y.  There must be at least two points, i.e. khi > klo.
+
+procedure tu_fit1 (xa, ya, x, klo, khi, y)
+
+double	xa[ARB]		# i: array of independent-variable values
+double	ya[ARB]		# i: array of dependent-variable values
+double	x		# i: independent variable
+int	klo, khi	# i: range of indices in xa covered by iv_step
+double	y		# o: value of fit at x
+#--
+double	sumx, sumy, sumxy, sumx2
+double	xmean, ymean
+double	slope, intercept
+int	k
+
+begin
+	sumx = 0.d0
+	sumy = 0.d0
+	do k = klo, khi {
+	    sumx = sumx + xa[k]
+	    sumy = sumy + ya[k]
+	}
+	xmean = sumx / (khi - klo + 1)
+	ymean = sumy / (khi - klo + 1)
+
+	sumxy = 0.d0
+	sumx2 = 0.d0
+	do k = klo, khi {
+	    sumxy = sumxy + (xa[k] - xmean) * (ya[k] - ymean)
+	    sumx2 = sumx2 + (xa[k] - xmean)**2
+	}
+	slope = sumxy / sumx2
+	intercept = (ymean - slope * xmean)
+
+	y = intercept + slope * x
+end
+
+
+# tuiep1 -- linear interpolation
+# Interpolate to get one value.  X is supposed to be between xa[1] and
+# xa[2], but it is not an error if x is outside that interval.  Note
+# that xa, ya are subarrays of the arrays with the same names elsewhere.
+
+procedure tuiep1 (xa, ya, x, y)
+
+double	xa[2]			# i: pair of independent-variable values
+double	ya[2]			# i: pair of dependent-variable values
+double	x			# i: independent variable
+double	y			# o: the interpolated value
+#--
+double	p			# fraction of distance between indep var val
+
+begin
+	if (x == xa[1]) {
+	    y = ya[1]
+	} else if (x == xa[2]) {
+	    y = ya[2]
+	} else {
+	    p = (x - xa[1]) / (xa[2] - xa[1])
+	    y = p * ya[2] + (1.-p) * ya[1]
+	}
+end
diff --git a/pkg/utilities/nttools/trebin/tutrim.x b/pkg/utilities/nttools/trebin/tutrim.x
new file mode 100644
index 00000000..e792c37a
--- /dev/null
+++ b/pkg/utilities/nttools/trebin/tutrim.x
@@ -0,0 +1,43 @@
+# This routine trims elements from the end of the xout array by
+# decrementing nelem.  Values are trimmed if they are INDEF or
+# equal to padvalue (and padvalue itself is not INDEF).  Trimming
+# starts at the end and stops with the first value that is not
+# INDEF and not equal to padvalue.
+#
+# Phil Hodge, 26-Apr-2000
+
+procedure tu_trim (xout, nelem, padvalue)
+
+double	xout[ARB]	# i: independent variable values
+int	nelem		# io: size of xout array
+double	padvalue	# i: trim these values at end of xout array
+#--
+int	curr_nelem	# current value of nelem
+int	i
+
+begin
+	curr_nelem = nelem
+
+	if (!IS_INDEFD(padvalue)) {
+
+	    # Check for either INDEF or padvalue.
+	    do i = curr_nelem, 1, -1 {
+		if (IS_INDEFD(xout[i]))
+		    nelem = nelem - 1
+		else if (xout[i] == padvalue)
+		    nelem = nelem - 1
+		else			# neither INDEF nor a pad value
+		    break
+	    }
+
+	} else {
+
+	    # Just check for INDEF.
+	    do i = curr_nelem, 1, -1 {
+		if (IS_INDEFD(xout[i]))
+		    nelem = nelem - 1
+		else			# not INDEF
+		    break
+	    }
+	}
+end
diff --git a/pkg/utilities/nttools/trebin/tuxget.x b/pkg/utilities/nttools/trebin/tuxget.x
new file mode 100644
index 00000000..9b890bb4
--- /dev/null
+++ b/pkg/utilities/nttools/trebin/tuxget.x
@@ -0,0 +1,134 @@
+include <tbset.h>
+
+# The values of the independent variable for the output table can be
+# read either from an input table (xtable), or they can be assigned
+# from the start, increment, and end values.
+#
+# Phil Hodge, 25-Apr-2000  Subroutine created.
+# Phil Hodge, 12-May-2004  errchk the procedures that are called.
+
+procedure tuxget (xtable, iv_start, iv_end, iv_step, padvalue,
+		xout, outrows)
+
+char	xtable[ARB]	# i: table of output indep var values
+double	iv_start	# i: starting value of independent variable
+double	iv_end		# i: ending value of independent variable
+double	iv_step		# i: increment in independent variable
+double	padvalue	# i: value at end of indep var array to ignore
+pointer xout		# o: pointer to output indep var values
+int	outrows		# o: number of output rows (size of xout)
+#--
+pointer tp, cp
+int	row, nrows
+int	i
+int	nelem		# array size, or number of table rows
+int	nvals		# number of elements read from column if it's an array
+double	dbuf
+double	direction	# indicates increasing or decreasing data
+double	previous	# for comparing current and previous values
+pointer tbtopn(), tbcnum()
+int	tbpsta(), tbcigi(), tbagtd()
+errchk	tbtopn, tbpsta, tbcnum, tbcigi, tbtclo, tbagtd, tbegtd, tu_trim
+
+begin
+	if (xtable[1] != EOS) {
+
+	    tp = tbtopn (xtable, READ_ONLY, NULL)
+	    if (tbpsta (tp, TBL_NCOLS) < 1) {
+		call tbtclo (tp)
+		call error (1, "No data in xtable")
+	    }
+	    if (tbpsta (tp, TBL_NCOLS) > 1) {
+		call tbtclo (tp)
+		call eprintf ("xtable %s contains more than one column;\n")
+		    call pargstr (xtable)
+		call eprintf (
+	"use a column selector [c:<colname>] to specify which column.\n")
+		call error (1, "")
+	    }
+
+	    nrows = tbpsta (tp, TBL_NROWS)
+	    cp = tbcnum (tp, 1)
+	    nelem = tbcigi (cp, TBL_COL_LENDATA)
+
+	    if (nelem > 1 && nrows > 1) {
+		call tbtclo (tp)
+		call eprintf ("xtable %s contains more than one row,\n")
+		    call pargstr (xtable)
+		call eprintf ("and the column contains arrays;\n")
+		call eprintf (
+	"use a row selector [c:row=<rownum>] to specify which row.\n")
+		call error (1, "")
+	    }
+
+	    if (nelem == 1)
+		nelem = nrows
+	    call malloc (xout, nelem, TY_DOUBLE)
+
+	    # Read the data from the table.
+	    if (nrows == 1) {
+		row = 1
+		nvals = tbagtd (tp, cp, row, Memd[xout], 1, nelem)
+		if (nvals < nelem) {
+		    call tbtclo (tp)
+		    call error (1, "not all elements read from xtable")
+		}
+	    } else {
+		do row = 1, nrows
+		    call tbegtd (tp, cp, row, Memd[xout+row-1])
+	    }
+
+	    call tbtclo (tp)
+
+	    # Trim trailing garbage by decrementing outrows.
+	    outrows = nelem
+	    call tu_trim (Memd[xout], outrows, padvalue)
+
+	    # Check for embedded INDEF values, and make sure the values
+	    # are monotonically increasing or decreasing.
+	    if (outrows > 1) {
+
+		do i = 1, outrows {
+		    if (IS_INDEFD(Memd[xout+i-1])) {
+			call eprintf (
+			"xtable %s contains embedded INDEF values\n")
+			    call pargstr (xtable)
+			call eprintf ("(i.e. not just trailing INDEFs)\n")
+			call error (1, "")
+		    }
+		}
+
+		if (Memd[xout+1] >= Memd[xout])
+		    direction = 1.d0
+		else
+		    direction = -1.d0
+		previous = Memd[xout] - direction
+		do i = 1, outrows {
+		    if (direction * (Memd[xout+i-1] - previous) <= 0.d0) {
+			call eprintf (
+			"Values in xtable %s are not monotonic\n")
+			    call pargstr (xtable)
+			call error (1, "")
+		    }
+		    previous = Memd[xout+i-1]
+		}
+	    }
+
+	} else {			# no xtable
+
+	    # Find out how many rows the output table should have.
+	    if (iv_start == iv_end)
+		outrows = 1
+	    else
+		outrows = nint ((iv_end - iv_start) / iv_step + 1.0)
+
+	    call malloc (xout, outrows, TY_DOUBLE)
+
+	    # Compute the independent variable values for the output table.
+	    dbuf = iv_start
+	    do i = 1, outrows {
+		Memd[xout+i-1] = dbuf
+		dbuf = dbuf + iv_step
+	    }
+	}
+end