Repatch (from linux) of OSX IRAF

11 files changed, 917 insertions, 0 deletions

diff --git a/pkg/utilities/nttools/tmatch/getmatch.x b/pkg/utilities/nttools/tmatch/getmatch.x
new file mode 100644
index 00000000..083d4156
--- /dev/null
+++ b/pkg/utilities/nttools/tmatch/getmatch.x
@@ -0,0 +1,101 @@
+#* HISTORY *
+#* B.Simon	24-Aug-94	Original
+#* B.Simon	18-Sep-00	Revised search termination criterion
+
+# GETMATCH -- Find rows in second table witch match rows in first
+
+procedure getmatch (in1, in2, ncol, col1, col2, weight, nrow1, index1, 
+		    nrow2, index2, maxnorm, sphere, closest, dist)
+
+pointer	in1		# i: first table descriptor
+pointer	in2		# i: second table descriptor
+int	ncol		# i: number of match columns
+pointer	col1[ARB]	# i: match columns in first table
+pointer	col2[ARB]	# i: match columns in second table
+double	weight[ARB]	# i: weights used in computing norm
+int	nrow1		# i: number of rows in first table
+int	index1[ARB]	# i: sorted row indices for first table
+int	nrow2		# i: number of rows in second table
+int	index2[ARB]	# i: sorted row indices for second table
+double	maxnorm		# i: maximum norm used in match
+bool	sphere		# i: apply spherical correction to first column?
+int	closest[ARB]	# o: closest match in second table to first
+double	dist[ARB]	# o: distance between matched rows
+#--
+double	sqnorm,  proj, abnorm, norm, minabnorm, minnorm
+int	idx, jdx, irow, jrow, krow, jlast
+
+begin
+	jlast = 1
+	sqnorm = maxnorm * maxnorm
+
+	# Find the row in the second table which minimizes the norm for
+	# each row of the first table
+
+
+	do idx = 1, nrow1 {
+	    irow = index1[idx]
+	    jrow = index2[jlast]
+
+	    # The initial guess is the row which matched last time
+
+	    call getnorm (in1, in2, ncol, col1, col2, irow, jrow, 
+			  weight, sphere, proj, abnorm, norm)
+
+	    minabnorm = abnorm
+	    minnorm = norm
+	    krow = jrow
+
+	    # Search backwards for a row which minimizes the norm
+	    # Terminate the search when the first dimension of the norm (proj)
+	    # is greater than the minimum norm, as all subsequent rows have
+	    # norms that must be greater than the minimum
+
+	    do jdx = jlast-1, 1, -1 {
+		jrow = index2[jdx]
+
+		call getnorm (in1, in2, ncol, col1, col2, irow, jrow,
+			      weight, sphere, proj, abnorm, norm)
+
+		if (proj > minabnorm)
+		    break
+
+		if (norm < minnorm) {
+		    minabnorm = abnorm
+		    minnorm = norm
+		    krow = jrow
+		}
+	    }
+
+	    # Search forwards for a row that minimizes the norm
+	    # Use the same termination condition as for the forwards search
+
+	    do jdx = jlast+1, nrow2 {
+		jrow = index2[jdx]
+
+		call getnorm (in1, in2, ncol, col1, col2, irow, jrow,
+			      weight, sphere, proj, abnorm, norm)
+
+		if (proj > minabnorm)
+		    break
+
+		if (norm < minnorm) {
+		    minabnorm = abnorm
+		    minnorm = norm
+		    krow = jrow
+		}
+	    }
+
+	    if (minnorm > sqnorm) {
+		dist[irow] = maxnorm
+		closest[irow] = 0
+	    } else {
+		dist[irow] = sqrt (minnorm)
+		closest[irow] = krow
+	    }
+
+	    jlast = krow
+	}
+
+end
+
diff --git a/pkg/utilities/nttools/tmatch/getnorm.x b/pkg/utilities/nttools/tmatch/getnorm.x
new file mode 100644
index 00000000..fed277c6
--- /dev/null
+++ b/pkg/utilities/nttools/tmatch/getnorm.x
@@ -0,0 +1,67 @@
+include	<math.h>
+
+#* HISTORY *
+#* B.Simon	24-Aug-94	original
+#* B.Simon	18-Sep-00	Revised computation of proj and added abnorm
+
+# GETNORM -- Compute the squared norm between two table rows
+
+procedure getnorm (in1, in2, ncol, col1, col2, row1, row2, weight, sphere, 
+		   proj, abnorm, norm)
+
+pointer	in1		# i: first table descriptor
+pointer	in2		# i: second table descriptor
+int	ncol		# i: number of match columns
+pointer	col1[ARB]	# i: match columns in first table
+pointer	col2[ARB]	# i: match columns in second table
+int	row1		# i: row number in first table
+int	row2		# i: row number in second table
+double	weight[ARB]	# i: weights used in computing norm
+bool	sphere		# i: apply spherical correction to first column?
+double	proj		# o: projection of norm on first axis
+double	abnorm		# o: norm, possibly without spherical correction
+double	norm		# o: norm (distance) between rows in two tables
+#--
+int	i
+double	val1, val2, dif
+
+begin
+	# Calculate first component of norm
+
+	call tbegtd (in1, col1[1], row1, val1)
+	call tbegtd (in2, col2[1], row2, val2)
+
+	dif = weight[1] * (val1 - val2)
+	proj = dif * dif
+	abnorm = proj
+
+	# Apply correction for spherical coordinates
+
+	if (sphere) {
+	    if (dif > 180) {
+		dif = dif - 360
+	    } else if (dif < -180) {
+		dif = dif + 360
+	    }
+
+	    call tbegtd (in1, col1[2], row1, val1)
+	    call tbegtd (in2, col2[2], row2, val2)
+
+	    val1 = 0.5 * weight[2] * (val1 + val2)
+	    dif = dif * cos (DEGTORAD(val1))
+	}
+
+	# Compute remaining components
+
+	norm = dif * dif
+	do i = 2, ncol {
+	    call tbegtd (in1, col1[i], row1, val1)
+	    call tbegtd (in2, col2[i], row2, val2)
+
+	    dif = weight[i] * (val1 - val2)
+	    abnorm = abnorm + dif * dif
+	    norm = norm + dif * dif
+	}
+
+end
+
diff --git a/pkg/utilities/nttools/tmatch/getweight.x b/pkg/utilities/nttools/tmatch/getweight.x
new file mode 100644
index 00000000..ef8b9204
--- /dev/null
+++ b/pkg/utilities/nttools/tmatch/getweight.x
@@ -0,0 +1,96 @@
+include	<math.h>
+include	<tbset.h>
+
+#* HISTORY *
+#* B.Simon	24-Aug-94	original
+
+# GETWEIGHT -- Get array of weights from list of factors or tables
+
+procedure getweight (ncol, col1, col2, factor, weight)
+
+int	ncol		# i: number of match columns
+pointer	col1[ARB]	# i: match columns from first table
+pointer	col2[ARB]	# i: match columns from second table
+char	factor[ARB]	# i: list of factors
+double	weight[ARB]	# o: array of weights
+#--
+double	unitval[6]
+int	invert[6]
+int	ic, jc, nc, icol, jcol, type1, type2, item
+pointer	sp, value, unit1, unit2, errmsg
+
+
+data	unitval	  / 1.0, 3600.0, 60.0, 1.0, 15.0, RADIAN /
+data	invert	  / NO,  YES,    YES,  NO,  NO,   NO     /
+
+string	unitlist  "|seconds|minutes|degrees|hours|radians|"
+string	badvalue  "Value in factor string is not a number (%s)"
+string	badunits  "Units mismatch in column %d of tables"
+
+int	ctod(), word_fetch(), strdic()
+
+begin
+	# Allocate memory for temporary strings
+
+	call smark (sp)
+	call salloc (value, SZ_FNAME, TY_CHAR)
+	call salloc (unit1, SZ_FNAME, TY_CHAR)
+	call salloc (unit2, SZ_FNAME, TY_CHAR)
+	call salloc (errmsg, SZ_LINE, TY_CHAR)
+
+	# Get each string from the list and convert to a number
+
+	ic = 1
+	icol = 0
+	while (word_fetch (factor, ic, Memc[value], SZ_FNAME) > 0) {
+	    icol = icol + 1
+
+	    jc = 1
+	    nc = ctod (Memc[value], jc, weight[icol])
+	    if (Memc[value+jc-1] != EOS) {
+		call sprintf (Memc[errmsg], SZ_LINE, badvalue)
+		call pargstr (Memc[value])
+
+		call error (1, Memc[errmsg])
+	    }
+	}
+
+	# Set remaining weights according to column units
+
+	do jcol = icol+1, ncol {
+	    # Read units from table
+
+	    call tbcigt (col1[jcol], TBL_COL_UNITS, Memc[unit1], SZ_FNAME)
+	    call tbcigt (col2[jcol], TBL_COL_UNITS, Memc[unit2], SZ_FNAME)
+
+	    # Search for units in dictionary
+
+	    call strlwr (Memc[unit1])
+	    call strlwr (Memc[unit2])
+
+	    type1 = strdic (Memc[unit1], Memc[unit1], SZ_FNAME, unitlist)
+	    type2 = strdic (Memc[unit2], Memc[unit2], SZ_FNAME, unitlist)
+
+	    # Take exit if units do not match
+
+	    if (type1 != type2) {
+		call sprintf (Memc[errmsg], SZ_LINE, badunits)
+		call pargi (jcol)
+
+		call error (1, Memc[errmsg])
+	    }
+
+	    # Read corresponding weight from unit value array
+	    # The first weight (1.0) is for missing or unknown units
+
+	    item = type1 + 1
+
+	    if (invert[item] == NO) {
+		weight[jcol] = unitval[item]
+	    } else {
+		weight[jcol] = 1.0 / unitval[item]
+	    }
+	}
+
+	call sfree (sp)
+end
diff --git a/pkg/utilities/nttools/tmatch/infomatch.x b/pkg/utilities/nttools/tmatch/infomatch.x
new file mode 100644
index 00000000..d90f89a8
--- /dev/null
+++ b/pkg/utilities/nttools/tmatch/infomatch.x
@@ -0,0 +1,219 @@
+include	<tbset.h>
+
+#* HISTORY *
+#* B.Simon	25-Aug-94	original
+
+# INFOMATCH -- Print diagnostic information for tmatch
+
+procedure infomatch (diagfile, in1, in2, nmcol1, nmcol2, maxnorm, 
+		     nclosest, closest, dist)
+
+char	diagfile[ARB]	# i: diagnostic output file
+pointer	in1		# i: first table's descriptor
+pointer	in2		# i: second table's descriptor
+char	nmcol1[ARB]	# i: name columns in first table
+char	nmcol2[ARB]	# i: name columns in second table
+double	maxnorm		# i: maximum allowed distance between matched rows
+int	nclosest	# i: length of closest array
+int	closest[ARB]	# i: array of closest matches between tables
+double	dist[ARB]	# i: distance between matched rows
+#--
+bool	first, same
+int	fd, namelen, mxcol1, mxcol2, ncol1, ncol2, idx, jdx, irow, jrow
+pointer	sp, index, name, col1, col2
+
+string	ziptitle  "\nThe following objects were not matched:\n"
+string	duptitle  "\nThe following objects matched the same object:\n"
+string	bigtitle  "\nThe following objects have the largest norms:\n"
+string	normfmt   "Norm = %0.7g\n"
+string	rowformat "%d:%d %s\n"
+
+bool	is_blank()
+int	open(), envgeti(), tbpsta()
+
+begin
+	# Open the diagnostics file
+
+	if (is_blank (diagfile))
+	    return
+
+	fd = open (diagfile, WRITE_ONLY, TEXT_FILE)
+
+	# Get maximum length of diagnostic string
+
+	iferr {
+	    namelen = envgeti ("ttyncols") - 10
+	} then {
+	    namelen = 70
+	}
+
+	# Allocate dynamic memory
+
+	call smark (sp)
+	call salloc (index, nclosest, TY_INT)
+	call salloc (name, namelen, TY_CHAR)
+
+	# Get column descriptors for name columns
+
+	mxcol1 = tbpsta (in1, TBL_NCOLS)
+	mxcol2 = tbpsta (in2, TBL_NCOLS)
+
+	call salloc (col1, mxcol1, TY_INT)
+	call salloc (col2, mxcol2, TY_INT)
+
+	if (is_blank (nmcol1)) {
+	    ncol1 = 0
+	} else {
+	    call tctexp (in1, nmcol1, mxcol1, ncol1, Memi[col1])
+	}
+
+	if (is_blank (nmcol2)) {
+	    ncol2 = 0
+	} else {
+	    call tctexp (in2, nmcol2, mxcol2, ncol2, Memi[col2])
+	}
+
+	# Sort the closest array
+
+	call setindex (Memi[index], nclosest)
+	call sortclose (nclosest, closest, Memi[index])
+
+	# Print the objects that were not matched
+
+	first = true
+	do idx = 1, nclosest {
+	    irow = Memi[index+idx-1]
+	    if (closest[irow] != 0)
+		break
+
+	    if (first) {
+		first = false
+		call fprintf (fd, ziptitle)
+	    }
+
+	    call rowname (in1, irow, ncol1, Memi[col1], Memc[name], namelen)
+	    call fprintf (fd, rowformat)
+	    call pargi (1)
+	    call pargi (irow)
+	    call pargstr (Memc[name])
+	}
+
+	# Print the objects which are matched more than once
+
+	same = false
+	first = true
+	do idx = 2, nclosest {
+	    irow = Memi[index+idx-1]
+	    jrow = Memi[index+idx-2]
+
+	    if (closest[irow] == 0)
+		next
+
+	    if (closest[irow] == closest[jrow]) {
+		same = true
+
+		if (first) {
+		    first = false
+		    call fprintf (fd, duptitle)
+		}
+
+		call rowname (in1, jrow, ncol1, Memi[col1], 
+			      Memc[name], namelen)
+
+		call fprintf (fd, rowformat)
+		call pargi (1)
+		call pargi (jrow)
+		call pargstr (Memc[name])
+
+	    } else if (same) {
+		same = false
+
+		call rowname (in1, jrow, ncol1, Memi[col1], 
+			      Memc[name], namelen)
+
+		call fprintf (fd, rowformat)
+		call pargi (1)
+		call pargi (jrow)
+		call pargstr (Memc[name])
+
+		call rowname (in2, closest[jrow], ncol2, Memi[col2],
+			      Memc[name], namelen)
+
+		call fprintf (fd, rowformat)
+		call pargi (2)
+		call pargi (closest[jrow])
+		call pargstr (Memc[name])
+
+		call fprintf (fd, "\n")
+	    }
+	}
+
+	if (same) {
+	    same = false
+	    irow = Memi[index+nclosest-1]
+
+	    call rowname (in1, irow, ncol1, Memi[col1], 
+			  Memc[name], namelen)
+
+	    call fprintf (fd, rowformat)
+	    call pargi (1)
+	    call pargi (irow)
+	    call pargstr (Memc[name])
+
+	    call rowname (in2, closest[irow], ncol2, Memi[col2],
+			  Memc[name], namelen)
+
+	    call fprintf (fd, rowformat)
+	    call pargi (2)
+	    call pargi (closest[irow])
+	    call pargstr (Memc[name])
+
+	    call fprintf (fd, "\n")
+	}
+
+	# Sort the dist array
+
+	call setindex (Memi[index], nclosest)
+	call sortdist (nclosest, dist, Memi[index])
+
+	# Print the ten objects with the largest norms
+
+	jdx = 0
+	do idx = nclosest, 1, -1 {
+	    irow = Memi[index+idx-1]
+	    if (dist[irow] == maxnorm)
+		next
+
+	    if (jdx == 0)
+		call fprintf (fd, bigtitle)
+
+	    jdx = jdx + 1
+	    if (jdx > 10)
+		break
+
+	    call fprintf (fd, normfmt)
+	    call pargd (dist[irow])
+
+	    call rowname (in1, irow, ncol1, Memi[col1], 
+			  Memc[name], namelen)
+
+	    call fprintf (fd, rowformat)
+	    call pargi (1)
+	    call pargi (irow)
+	    call pargstr (Memc[name])
+
+	    call rowname (in2, closest[irow], ncol2, Memi[col2],
+			  Memc[name], namelen)
+
+	    call fprintf (fd, rowformat)
+	    call pargi (2)
+	    call pargi (closest[irow])
+	    call pargstr (Memc[name])
+
+	    call fprintf (fd, "\n")
+
+	}
+
+	call close (fd)
+	call sfree (sp)
+end
diff --git a/pkg/utilities/nttools/tmatch/mkpkg b/pkg/utilities/nttools/tmatch/mkpkg
new file mode 100644
index 00000000..b6e6be18
--- /dev/null
+++ b/pkg/utilities/nttools/tmatch/mkpkg
@@ -0,0 +1,20 @@
+# Update the tmatch application code in the ttools package library
+# Author: B.Simon, 30-Aug-94
+
+$checkout libpkg.a ../
+$update   libpkg.a
+$checkin  libpkg.a ../
+$exit
+
+libpkg.a:
+	getmatch.x	
+	getnorm.x	<math.h>
+	getweight.x	<math.h>  <tbset.h>
+	infomatch.x	<tbset.h>
+	putmatch.x	<tbset.h>
+	rowname.x	
+	setindex.x	
+	sortclose.x	
+	sortdist.x	
+	tmatch.x	<tbset.h>
+	;
diff --git a/pkg/utilities/nttools/tmatch/putmatch.x b/pkg/utilities/nttools/tmatch/putmatch.x
new file mode 100644
index 00000000..27407107
--- /dev/null
+++ b/pkg/utilities/nttools/tmatch/putmatch.x
@@ -0,0 +1,102 @@
+include	<tbset.h>
+
+#* HISTORY *
+# B.Simon	25-Aug-94	Original
+
+# PUTMATCH -- Write matched rows in input as a single row in output table
+
+procedure putmatch (output, incol1, incol2, in1, in2, nclosest, closest)
+
+char	output[ARB]	# i: output table name
+char	incol1[ARB]	# i: list of columns to copy from first table
+char	incol2[ARB]	# i: list of columns to copy from second table
+pointer	in1		# i: first table's descriptor
+pointer	in2		# i: second table's descriptor
+int	nclosest	# i: length of closest array
+int	closest[ARB]	# i: indices of rows in second table closest to first
+#--
+int	mxcol1, mxcol2, maxcol, ncol1, ncol2, ncol, type1, type2
+pointer	colnum, datatype, lendata, lenfmt, icol, irow, jrow
+pointer	sp, colname, colunits, colfmt, oldcol, newcol,out
+
+string	nomatch "WARNING: No rows matched between tables, output \
+table is empty\n"
+
+int	tbpsta()
+pointer	tbtopn()
+
+begin
+	# Allocate memory for temporary strings
+
+	call smark (sp)
+        call salloc (colname, SZ_COLNAME, TY_CHAR)
+        call salloc (colunits, SZ_COLUNITS, TY_CHAR)
+        call salloc (colfmt, SZ_COLFMT, TY_CHAR)
+
+	# Get column descriptors from input tables
+
+	mxcol1 = tbpsta (in1, TBL_NCOLS)
+	mxcol2 = tbpsta (in2, TBL_NCOLS)
+	maxcol = mxcol1 + mxcol2
+
+	call salloc (oldcol, maxcol, TY_INT)
+	call salloc (newcol, maxcol, TY_INT)
+
+	call tctexp (in1, incol1, mxcol1, ncol1, Memi[oldcol])
+	call tctexp (in2, incol2, mxcol2, ncol2, Memi[oldcol+ncol1])
+	ncol = ncol1 + ncol2
+
+	# Create output table
+
+	out = tbtopn (output, NEW_FILE, NULL)
+
+	# Set type (text, row ordered, column ordered)
+
+	type1 = tbpsta (in1, TBL_WHTYPE)
+	type2 = tbpsta (in2, TBL_WHTYPE)
+	if (type1 == type2)
+	    call tbpset (out, TBL_WHTYPE, type1)
+
+	# Create columns in output table
+
+	do icol = 1, ncol {
+	    call tbcinf (Memi[oldcol+icol-1], colnum, Memc[colname], 
+			 Memc[colunits], Memc[colfmt], datatype, 
+			 lendata, lenfmt)
+
+	    call newcolnam (ncol, Memi[oldcol], icol, 
+			    Memc[colname], SZ_COLNAME)
+
+	    call tbcdef (out, Memi[newcol+icol-1], Memc[colname], 
+			 Memc[colunits], Memc[colfmt], datatype, lendata, 1)
+	}
+
+	# Copy header keywords from first input table
+
+	call tbtcre (out)
+	call tbhcal (in1, out)
+
+	# Copy rows from input table to output
+
+	jrow = 0
+	do irow = 1, nclosest {
+	    if (closest[irow] == 0)
+		next
+
+	    jrow = jrow + 1
+	    call tbrcsc (in1, out, Memi[oldcol], Memi[newcol], 
+			 irow, jrow, ncol1)
+	    call tbrcsc (in2, out, Memi[oldcol+ncol1], Memi[newcol+ncol1],
+			 closest[irow], jrow, ncol2)
+	}
+
+	# Write warning message if no rows matched
+
+	if (jrow == 0)
+	    call eprintf (nomatch)
+
+	# Clean up
+
+	call tbtclo (out)
+	call sfree (sp)
+end
diff --git a/pkg/utilities/nttools/tmatch/rowname.x b/pkg/utilities/nttools/tmatch/rowname.x
new file mode 100644
index 00000000..484fcd75
--- /dev/null
+++ b/pkg/utilities/nttools/tmatch/rowname.x
@@ -0,0 +1,61 @@
+#* HISTORY *
+#* B.Simon	26-Aug-94
+
+# ROWNAME -- Create name for table row by concatenating column values
+
+procedure rowname (in, irow, ncol, col, name, namelen)
+
+pointer	in		# i: table descriptor
+int	irow		# i: table row number
+int	ncol		# i: number of table columns
+pointer	col[ARB]	# i: table column pointers
+char	name[ARB]	# o: concatenated values of columns
+int	namelen		# i: maximum name length
+#--
+int	ic, jc, icol
+pointer	sp, value
+
+begin
+	# Allocate memory for column buffer
+
+	call smark (sp)
+	call salloc (value, SZ_LINE, TY_CHAR)
+
+	# Concatenate column values into name string
+
+	jc = 0
+	icol = 0
+	for (ic = 1; ic <= namelen; ic = ic + 1) {
+
+	    # A value of zero is a flag to read the next coumn
+
+	    if (jc == 0) {
+		icol = icol + 1
+		if (icol > ncol) {
+		    if (ic > 1)
+			ic = ic - 1	# remove trailing blank
+
+		    break
+		}
+
+		call tbegtt (in, col[icol], irow, Memc[value], SZ_LINE)
+	    }
+
+	    # Copy a single character from the buffer to the output string
+	    # until the buffer is exhausted. At this point we copy a blank
+	    # as a spacer and set the value of jc as a flag to read the
+	    # next column.
+
+	    if (Memc[value+jc] == EOS) {
+		name[ic] = ' '
+		jc = 0
+	    } else {
+		name[ic] = Memc[value+jc]
+		jc = jc + 1
+	    }
+	}
+
+	name[ic] = EOS
+	call sfree (sp)
+end
+
diff --git a/pkg/utilities/nttools/tmatch/setindex.x b/pkg/utilities/nttools/tmatch/setindex.x
new file mode 100644
index 00000000..9640f1cd
--- /dev/null
+++ b/pkg/utilities/nttools/tmatch/setindex.x
@@ -0,0 +1,13 @@
+# SETINDEX -- Initialize an index array
+
+procedure setindex (index, len)
+
+int	index[ARB]	# o: index rray
+int	len		# i: array length
+#--
+int	i
+
+begin
+	do i = 1, len
+	    index[i] = i
+end
diff --git a/pkg/utilities/nttools/tmatch/sortclose.x b/pkg/utilities/nttools/tmatch/sortclose.x
new file mode 100644
index 00000000..303ac250
--- /dev/null
+++ b/pkg/utilities/nttools/tmatch/sortclose.x
@@ -0,0 +1,50 @@
+#* HISTORY *
+#* B.Simon	25-Aug-94	original
+#  Phil Hodge	12-Jul-2005	in sortclose, add 'int cmpclose()'
+
+# SORTCLOSE -- Sort the closest array 
+
+procedure sortclose (nclosest, closest, index)
+
+int	nclosest	# i: length of closest and index arrays
+int	closest[ARB]	# i: indices of second table rows matching first
+int	index[ARB]	# u: indices of first table rows, in sort order on exit
+#--
+int	cmpclose()
+extern	cmpclose
+pointer	sp, close2
+
+begin
+	call smark (sp)
+	call salloc (close2, nclosest, TY_INT)
+
+	call amovi (closest, Memi[close2], nclosest)
+	call gqsort (index, nclosest, cmpclose, close2)
+
+	call sfree (sp)
+end
+
+# CMPCLOSE -- Compare two elements in the close array
+
+int procedure cmpclose (closest, ielem, jelem)
+
+pointer	closest		# address of the closest array
+int	ielem		# first element
+int	jelem		# second element
+#--
+int	order
+
+begin
+	if (Memi[closest+ielem-1] < Memi[closest+jelem-1]) {
+	    order = -1
+	} else if (Memi[closest+ielem-1] > Memi[closest+jelem-1]) {
+	    order = 1
+	} else if (ielem < jelem) {
+	    order = -1
+	} else if (ielem < jelem) {
+	    order = 1
+	} else {
+	    order = 0
+	}
+	return (order)
+end
diff --git a/pkg/utilities/nttools/tmatch/sortdist.x b/pkg/utilities/nttools/tmatch/sortdist.x
new file mode 100644
index 00000000..b0745676
--- /dev/null
+++ b/pkg/utilities/nttools/tmatch/sortdist.x
@@ -0,0 +1,50 @@
+#* HISTORY *
+#* B.Simon	25-Aug-94	original
+#  Phil Hodge	28-Sept-2005	in sortdist, add 'int cmpdist()'
+
+# SORTDIST -- Sort the dist array 
+
+procedure sortdist (ndist, dist, index)
+
+int	ndist		# i: length of dist and index arrays
+double	dist[ARB]	# i: indices of second table rows matching first
+int	index[ARB]	# u: indices of first table rows, in sort order on exit
+#--
+int	cmpdist()
+extern	cmpdist
+pointer	sp, dist2
+
+begin
+	call smark (sp)
+	call salloc (dist2, ndist, TY_DOUBLE)
+
+	call amovd (dist, Memd[dist2], ndist)
+	call gqsort (index, ndist, cmpdist, dist2)
+
+	call sfree (sp)
+end
+
+# CMPDIST -- Compare two elements in the dist array
+
+int procedure cmpdist (dist, ielem, jelem)
+
+pointer	dist		# address of the dist array
+int	ielem		# first element
+int	jelem		# second element
+#--
+int	order
+
+begin
+	if (Memd[dist+ielem-1] < Memd[dist+jelem-1]) {
+	    order = -1
+	} else if (Memd[dist+ielem-1] > Memd[dist+jelem-1]) {
+	    order = 1
+	} else if (ielem < jelem) {
+	    order = -1
+	} else if (ielem < jelem) {
+	    order = 1
+	} else {
+	    order = 0
+	}
+	return (order)
+end
diff --git a/pkg/utilities/nttools/tmatch/tmatch.x b/pkg/utilities/nttools/tmatch/tmatch.x
new file mode 100644
index 00000000..2d3ea22d
--- /dev/null
+++ b/pkg/utilities/nttools/tmatch/tmatch.x
@@ -0,0 +1,138 @@
+include	<tbset.h>
+
+#* HISTORY *
+#* B.Simon	24-Aug-1994	original
+# Phil Hodge	 8-Apr-1999	Call tbfpri.
+
+# TMATCH -- Find closest matching rows between two tables
+
+procedure tmatch ()
+
+#--
+pointer	input1		# First input table
+pointer	input2		# Second input table
+pointer	output		# Output table
+pointer	match1		# Columns from first table used to match
+pointer	match2		# Columns from second table used to match
+double	maxnorm		# Maximum value of norm for allowed match
+pointer	incol1		# Columns from first table copied to output
+pointer	incol2		# Columns from second table copied to output
+pointer	factor		# Multiplicative factors used in computing norm
+pointer	diagfile	# Diagnostic output file
+pointer	nmcol1		# Columns from first table in diagnostic output
+pointer	nmcol2		# Columns from second table in diagnostic output
+bool	sphere		# Apply spherical correction to first column?
+
+bool	fold
+int	mxcol1, mxcol2, ncol1, ncol2, nrow1, nrow2
+int	phu_copied	# set by tbfpri and ignored
+pointer	sp, in1, in2, col1, col2, index1, index2, weight, dist, closest
+
+data	fold	/ false /
+
+string	mismatch  "Both lists of match columns must have same length"
+string	nomatch   "Match columns not found in table"
+
+bool	clgetb()
+double	clgetd()
+int	tbpsta()
+pointer	tbtopn()
+
+begin
+	# Allocate memory for strings
+
+	call smark (sp)
+	call salloc (input1, SZ_FNAME, TY_CHAR)
+	call salloc (input2, SZ_FNAME, TY_CHAR)
+	call salloc (output, SZ_FNAME, TY_CHAR)
+	call salloc (match1, SZ_FNAME, TY_CHAR)
+	call salloc (match2, SZ_FNAME, TY_CHAR)
+	call salloc (incol1, SZ_FNAME, TY_CHAR)
+	call salloc (incol2, SZ_FNAME, TY_CHAR)
+	call salloc (factor, SZ_FNAME, TY_CHAR)
+	call salloc (diagfile, SZ_FNAME, TY_CHAR)
+	call salloc (nmcol1, SZ_FNAME, TY_CHAR)
+	call salloc (nmcol2, SZ_FNAME, TY_CHAR)
+
+	# Read task parameters
+
+	call clgstr ("input1", Memc[input1], SZ_FNAME)
+	call clgstr ("input2", Memc[input2], SZ_FNAME)
+	call clgstr ("output", Memc[output], SZ_FNAME)
+	call clgstr ("match1", Memc[match1], SZ_FNAME)
+	call clgstr ("match2", Memc[match2], SZ_FNAME)
+	maxnorm = clgetd ("maxnorm")
+
+	call clgstr ("incol1", Memc[incol1], SZ_FNAME)
+	call clgstr ("incol2", Memc[incol2], SZ_FNAME)
+	call clgstr ("factor", Memc[factor], SZ_FNAME)
+	call clgstr ("diagfile", Memc[diagfile], SZ_FNAME)
+	call clgstr ("nmcol1", Memc[nmcol1], SZ_FNAME)
+	call clgstr ("nmcol2", Memc[nmcol2], SZ_FNAME)
+	sphere = clgetb ("sphere")
+
+	# Open input tables and get list of match colums
+
+	in1 = tbtopn (Memc[input1], READ_ONLY, NULL)
+	in2 = tbtopn (Memc[input2], READ_ONLY, NULL)
+
+	mxcol1 = tbpsta (in1, TBL_NCOLS)
+	mxcol2 = tbpsta (in2, TBL_NCOLS)
+
+	call salloc (col1, mxcol1, TY_INT)
+	call salloc (col2, mxcol2, TY_INT)
+
+	call tctexp (in1, Memc[match1], mxcol1, ncol1, Memi[col1])
+	call tctexp (in2, Memc[match2], mxcol2, ncol2, Memi[col2])
+
+	if (ncol1 != ncol2)
+	    call error (1, mismatch)
+
+	if (ncol1 == 0)
+	    call error (1, nomatch)
+
+	if (ncol1 < 2)
+	    sphere = false
+
+	# Sort input tables
+
+	call allrows (in1, nrow1, index1)
+	call allrows (in2, nrow2, index2)
+	
+	call tbtsrt (in1, ncol1, Memi[col1], fold, nrow1, Memi[index1])
+	call tbtsrt (in2, ncol2, Memi[col2], fold, nrow2, Memi[index2])
+
+	call salloc (weight, ncol1, TY_DOUBLE)
+	call salloc (dist, nrow1, TY_DOUBLE)
+	call salloc (closest, nrow1, TY_INT)
+
+	# Compute weights from list of factors or table column units
+
+	call getweight (ncol1, Memi[col1], Memi[col2], 
+			Memc[factor], Memd[weight])
+
+	# Compute closest match between the two tables
+
+	call getmatch (in1, in2, ncol1, Memi[col1], Memi[col2], Memd[weight], 
+		       nrow1, Memi[index1], nrow2, Memi[index2], maxnorm, 
+		       sphere, Memi[closest], Memd[dist])
+
+	# Write output table
+
+	call tbfpri (Memc[input1], Memc[output], phu_copied)
+	call putmatch (Memc[output], Memc[incol1], Memc[incol2], in1, in2, 
+		       nrow1, Memi[closest])
+
+	# Write diagnostic info
+
+	call infomatch (Memc[diagfile], in1, in2, Memc[nmcol1], Memc[nmcol2], 
+			maxnorm, nrow1, Memi[closest], Memd[dist])
+
+	# Clean up
+
+	call mfree (index1, TY_INT)
+	call mfree (index2, TY_INT)
+	call tbtclo (in1)
+	call tbtclo (in2)
+	call sfree (sp)
+end