Initial commit

41 files changed, 3663 insertions, 0 deletions

diff --git a/math/surfit/doc/iscoeff.hlp b/math/surfit/doc/iscoeff.hlp
new file mode 100644
index 00000000..ea0b292c
--- /dev/null
+++ b/math/surfit/doc/iscoeff.hlp
@@ -0,0 +1,63 @@
+.help iscoeff Apr85 "Surfit Package"
+.ih
+NAME
+iscoeff - get the number and values of the surface coefficients
+.ih
+SYNOPSIS
+iscoeff (sf, coeff, ncoeff)
+
+.nf
+pointer	sf		# surface descriptor
+real	coeff[ncoeff]	# coefficient array
+int	ncoeff		# number of coefficients
+.fi
+.ih
+ARGUMENTS
+.ls sf   
+Pointer to the surface descriptor.
+.le
+.ls coeff
+Array of coefficients.
+.le
+.ls ncoeff
+The number of coefficients.
+.le
+
+.ih
+DESCRIPTION
+ISCOEFF fetches the coefficient array and the number of coefficients from
+the surface descriptor structure. A 1-D array ncoeff elements long
+is required to hold the coefficients where ncoeff is defined as follows:
+
+.nf
+	SF_LEGENDRE:		nxcoeff = xorder
+				nycoeff = yorder
+				ncoeff = xorder * yorder        xterms = yes
+				ncoeff = nycoeff + nxcoeff - 1  xterms = no
+
+	SF_CHEBYSHEV:		nxcoeff = xorder
+				nycoeff = yorder
+				ncoeff = xorder * yorder        xterms = yes
+				ncoeff = nycoeff + nxcoeff - 1  xterms = no
+
+	SF_SPLINE3:		nxcoeff = xorder + 3
+				nycoeff = yorder + 3
+				ncoeff = (xorder + 3) * (yorder + 3)
+
+	SF_SPLINE1:		nxcoeff = xorder + 1
+				nycoeff = yorder + 1
+				ncoeff = (xorder + 1) * (yorder + 1)
+
+.fi
+
+The coefficient of basis function B(i,x) * B(j,y) will be stored in
+element (i - 1) *
+nycoeff + j of the array coeff if the xterms parameter was set
+to yes by ISINIT. Otherwise the nycoeff y-term coefficients will be output
+first followed by the (nxcoeff - 1) x-term coefficients.
+
+.ih
+NOTES
+.ih
+SEE ALSO
+.endhelp
diff --git a/math/surfit/doc/iseval.hlp b/math/surfit/doc/iseval.hlp
new file mode 100644
index 00000000..53853cfe
--- /dev/null
+++ b/math/surfit/doc/iseval.hlp
@@ -0,0 +1,34 @@
+.help iseval Apr85 "Surfit Package"
+.ih
+NAME
+iseval -- evaluate the fitted surface at x and y
+.ih
+SYNOPSIS
+y = iseval (sf, x, y)
+
+.nf
+pointer	sf		# surface descriptor
+real	x		# x value, 1 <= x <= ncols
+real	y		# y value, 1 <= y <= nlines
+.fi
+.ih
+ARGUMENTS
+.ls sf      
+Pointer to the surface descriptor structure.
+.le
+.ls x, y
+X and y values at which the surface is to be evaluated.
+.le
+.ih
+DESCRIPTION
+Evaluate the surface at the specified value of x and y. ISEVAL is a real
+valued function which returns the fitted value.
+.ih
+NOTES
+ISEVAL uses the coefficient array stored in the surface descriptor structure.
+Checking for out of bounds x and y values is the responsibility of the calling
+program.
+.ih
+SEE ALSO
+isvector
+.endhelp
diff --git a/math/surfit/doc/isfree.hlp b/math/surfit/doc/isfree.hlp
new file mode 100644
index 00000000..04f52b5f
--- /dev/null
+++ b/math/surfit/doc/isfree.hlp
@@ -0,0 +1,27 @@
+.help isfree Apr85 "Surfit Package"
+.ih
+NAME
+isfree -- free the surface grid descriptor structure
+.ih
+SYNOPSIS
+isfree (sf)
+
+.nf
+pointer	sf		# surface descriptor
+.fi
+.ih
+ARGUMENTS
+.ls sf    
+Pointer to the surface descriptor structure.
+.le
+.ih
+DESCRIPTION
+Frees the surface descriptor structure and all the vectors and arrays used
+by the numerical routines.
+.ih
+NOTES
+ISFREE should be called after the surface fit is complete.
+.ih
+SEE ALSO
+isinit
+.endhelp
diff --git a/math/surfit/doc/isinit.hlp b/math/surfit/doc/isinit.hlp
new file mode 100644
index 00000000..e0f2123c
--- /dev/null
+++ b/math/surfit/doc/isinit.hlp
@@ -0,0 +1,61 @@
+.help isinit Apr85 "Surfit Package"
+.ih
+NAME
+isinit -- initialize surface grid descriptor
+.ih
+SYNOPSIS
+include <math/surfit.h>
+
+.nf
+isinit (sf, surface_type, xorder, yorder, xterms, ncols, nlines)
+.fi
+
+.nf
+pointer	sf		# surface descriptor
+int	surface_type	# surface function
+int	xorder		# order of function in x
+int	yorder		# order of function in y
+int	xterms		# include cross-terms? (YES/NO)
+int	ncols		# number of columns in the surface grid
+int	nlines		# number of lines in the surface grid
+.fi
+.ih
+ARGUMENTS
+.ls sf    
+Pointer to the surface descriptor structure.
+.le
+.ls surface_type
+Fitting function. Permitted values are SF_LEGENDRE and SF_CHEBYSHEV for
+the Legendre and Chebyshev polynomials and SF_SPLINE1 and SF_SPLINE3
+for the linear and bicubic splines.
+.le
+.ls xorder, yorder
+Order of the polynomial to be fit in x and y or the number of spline pieces to 
+be fit in x and y. The orders must be greater than or equal to 1.
+.le
+.ls xterms
+Include cross-terms? If xterms = YES coefficients are fit to terms containing
+the cross products of x and y polynomials. Xterms defaults to YES for the
+spline functions.
+.le
+.ls ncols
+The number of columns in the surface grid. The surface is assumed to lie
+on a rectangular grid such that 1 <= x <= ncols.
+.le
+.ls nlines
+The number of lines in the surface to be grid. The surface is assumed to lie
+on a rectangular grid such that 1 <= y <= nlines.
+.le
+.ih
+DESCRIPTION
+ISINIT allocates space for the surface descriptor and the arrays and vectors
+used by the numerical routines. It initializes all arrays and vectors to zero,
+calculates and stores the basis functions in x and y
+and returns the surface descriptor to the calling routine.
+.ih
+NOTES
+ISINIT must be the first SURFIT routine called.
+.ih
+SEE ALSO
+isfree
+.endhelp
diff --git a/math/surfit/doc/islaccum.hlp b/math/surfit/doc/islaccum.hlp
new file mode 100644
index 00000000..48481cf0
--- /dev/null
+++ b/math/surfit/doc/islaccum.hlp
@@ -0,0 +1,60 @@
+.help islaccum Apr85 "Surfit Package"
+.ih
+NAME
+islaccum -- accumulate a surface grid line into the fit
+.ih
+SYNOPSIS
+include <math/surfit.h>
+
+islaccum (sf, cols, lineno, z, w, npts, wtflag)
+
+.nf
+pointer	sf		# surface grid descriptor
+int	cols[npts]	# column values, 1 <= col[i] <= ncols
+int	lineno		# number of line, 1 <= lineno <= nlines 
+real	z[npts]		# surface values on lineno at cols
+real	w[npts]		# array of weights
+int	npts		# number of surface values, npts <= ncols
+int	wtflag		# type of weighting desired
+.fi
+.ih
+ARGUMENTS
+.ls sf      
+Pointer to the surface grid descriptor structure.
+.le
+.ls cols
+The column numbers of surface grid points on line lineno to be added to the
+dataset.
+.le
+.ls lineno
+The line number of the surface grid line to be added to the data set.
+.le
+.ls z      
+The data values.
+.le
+.ls w
+The array of weights for the data points.
+.le
+.ls npts
+The number of surface grid values on line number lineno to be added to the
+data set.
+.le
+.ls wtflag
+Type of weighting. The options are SF_USER and SF_UNIFORM. If wtflag
+equals SF_USER the weight for each data point is supplied by the user
+and points with zero-valued weights are not included in the fit.
+If wtflag equals SF_UNIFORM the routine sets the weights to 1.
+.le
+.ih
+DESCRIPTION
+ISLACCUM computes the contribution of each data point to the normal
+equations in x, sums that contribution into the appropriate arrays and
+vectors and stores these intermediate results for use by ISLSOLVE.
+.ih
+NOTES
+Checking for out of bounds col values and INDEF valued data points is
+the responsibility of the calling program.
+.ih
+SEE ALSO
+isinit, islfit, islrefit, islsolve, islzero, issolve, isfree
+.endhelp
diff --git a/math/surfit/doc/islfit.hlp b/math/surfit/doc/islfit.hlp
new file mode 100644
index 00000000..5adf0a59
--- /dev/null
+++ b/math/surfit/doc/islfit.hlp
@@ -0,0 +1,67 @@
+.help islfit Apr85 "Surfit Package"
+.ih
+NAME
+islfit -- fit a surface grid line
+.ih
+SYNOPSIS
+include <math/surfit.h>
+
+islfit (sf, cols, lineno, z, w, npts, wtflag, ier)
+
+.nf
+pointer	sf		# surface grid descriptor
+real	cols[npts]	# array of column numbers, 1 <= cols[i] <= ncols
+int	lineno		# number of surface grid line to be added
+real	z[npts]		# data values
+real	w[npts]		# weight array
+int	npts		# number of data points, npts <= ncols
+int	wtflag		# type of weighting
+int	ier		# error code
+.fi
+.ih
+ARGUMENTS
+.ls sf    
+Pointer to the surface grid descriptor structure.
+.le
+.ls cols
+The column numbers of surface grid points to be added to the dataset.
+.le
+.ls lineno
+The line number of the surface grid line to be added to the dataset.
+.le
+.ls z      
+Array of data values.
+.le
+.ls w
+Array of weights.
+.le
+.ls npts
+Number of data points.
+.le
+.ls wtflag
+Type of weighting. The options are SF_USER and SF_UNIFORM. If wtflag =
+SF_USER individual weights for each data point are supplied by the calling
+program and points with zero-valued weights are not included in the fit.
+If wtflag = SF_UNIFORM, all weights are assigned values of 1.
+.le
+.ls ier     
+Error code for the fit. The options are OK, SINGULAR and NO_DEG_FREEDOM.
+If ier = SINGULAR, the numerical routines will compute a solution but one
+or more of the coefficients will be zero. If ier = NO_DEG_FREEDOM there
+were too few data points to solve the matrix equations and the routine
+returns without fitting the data.
+.le
+.ih
+DESCRIPTION
+ISLFIT zeroes the appropriate arrays and vectors,
+computes the contribution of each data point to the normal equations
+in x and accumulates it into the appropriate array and vector elements.  The
+x coefficients are stored for later use by ISSOLVE.
+.ih
+NOTES
+Checking for out of bounds col values and INDEF values pixels is the
+responsibility of the user.
+.ih
+SEE ALSO
+isinit, islrefit, islaccum, islsolve, islzero, issolve, isfree
+.endhelp
diff --git a/math/surfit/doc/islrefit.hlp b/math/surfit/doc/islrefit.hlp
new file mode 100644
index 00000000..c34476d6
--- /dev/null
+++ b/math/surfit/doc/islrefit.hlp
@@ -0,0 +1,51 @@
+.help islrefit Aug85 "Gsurfit Package"
+.ih
+NAME
+islrefit -- refit surface grid line  assuming old cols, w vector
+.ih
+SYNOPSIS
+include <math/surfit.h>
+
+islrefit (sf, cols, lineno, z, w)
+
+.nf
+pointer	sf			# surface grid descriptor
+int	cols[ARB]		# columns to be fit, 1 <= cols[i] <= ncols
+int	lineno			# line number of surface grid line to be fit
+real	z[ARB]			# array of data values
+real	w[ARB]			# array of weights
+.fi
+.ih
+ARGUMENTS
+.ls sf          
+Pointer to the surface grid descriptor structure.
+.le
+.ls cols
+The column numbers of the surface grid line to be added to the dataset.
+.le
+.ls lineno
+The number of the surface grid line to be added to the dataset.
+.le
+.ls z    
+Array of data values.
+.le
+.ls w    
+Array of weights.
+.le
+.ih
+DESCRIPTION
+In some applications the cols, and w values remain unchanged from fit
+to fit and only the z values vary. In this case it is redundant to
+reaccumulate the matrix and perform the Cholesky factorization for each
+succeeding surface grid line. ISLREFIT
+zeros and reaccumulates the vector on the right hand side of the matrix
+equation and performs the forward and back substitution phase to fit for
+a new coefficient vector.
+.ih
+NOTES
+It is the responsibility of the calling program to check for out of bounds
+column values and INDEF valued pixels.
+.ih
+SEE ALSO
+isinit, islfit, islaccum, islsolve, islzero, issolve, isfree
+.endhelp
diff --git a/math/surfit/doc/islsolve.hlp b/math/surfit/doc/islsolve.hlp
new file mode 100644
index 00000000..752cbc19
--- /dev/null
+++ b/math/surfit/doc/islsolve.hlp
@@ -0,0 +1,45 @@
+.help islsolve Apr85 "Surfit Package"
+.ih
+NAME
+islsolve -- solve the equations for a surface grid line
+.ih
+SYNOPSIS
+include <math/surfit.h>
+
+islsolve (sf, lineno, ier)
+
+.nf
+pointer	sf		# surface grid descriptor
+int	lineno		# surface grid line number
+int	ier		# error code
+.fi
+.ih
+ARGUMENTS
+.ls sf     
+Pointer to the surface grid descriptor structure.
+.le
+.ls lineno
+The number of the surface grid line to be added to the dataset.
+.le
+.ls ier      
+Error code returned by the fitting routines. The options are OK, SINGULAR,
+and NO_DEG_FREEDOM. If ier = SINGULAR the matrix is singular, ISLSOLVE
+will compute a solution to the normal equations but one or more of the
+coefficients will be zero. If ier = NO_DEG_FREEDOM, too few data points
+exist for a reasonable solution to be computed and ISLSOLVE returns without
+fitting the data.
+.le
+.ih
+DESCRIPTION
+ISLSOLVE computes the Cholesky factorization of the data matrix and
+solves for the coefficients
+of the x fitting function by forward and back substitution.
+An error code is
+returned by ISLSOLVE if it is unable to solve the normal equations as
+formulated.
+.ih
+NOTES
+.ih
+SEE ALSO
+isinit, islfit, islrefit, islaccum, islzero, issolve, isfree
+.endhelp
diff --git a/math/surfit/doc/islzero.hlp b/math/surfit/doc/islzero.hlp
new file mode 100644
index 00000000..cbdb9515
--- /dev/null
+++ b/math/surfit/doc/islzero.hlp
@@ -0,0 +1,32 @@
+.help islzero Apr85 "Surfit Package"
+.ih
+NAME
+islzero -- set up for a new surface grid line fit
+.ih
+SYNOPSIS
+islzero (sf, lineno)
+
+.nf
+pointer	sf		# surface grid descriptor
+int	lineno		# surface grid line number
+.fi
+.ih
+ARGUMENTS
+.ls sf     
+Pointer to the surface grid descriptor structure.
+.le
+.ls lineno
+Line number of the surface grid line fit to be zeroed.
+.le
+.ih
+DESCRIPTION
+ISLZERO zeros the appropriate arrays and vectors for the surface grid line
+number line number. The line can then be refit by calls to ISLACCUM and
+ISLSOLVE. Alternatively a single call to ISLFIT can perform the combined
+functions of ISLZERO, ISLACCUM and ISLSOLVE.
+.ih
+NOTES
+.ih
+SEE ALSO
+isinit, islfit, islrefit, islaccum, islsolve, issolve
+.endhelp
diff --git a/math/surfit/doc/isreplace.hlp b/math/surfit/doc/isreplace.hlp
new file mode 100644
index 00000000..d248e5f8
--- /dev/null
+++ b/math/surfit/doc/isreplace.hlp
@@ -0,0 +1,33 @@
+.help isreplace Apr1985 "Surfit Package"
+.ih
+NAME
+isreplace -- restore a surface fit saved by issave
+.ih
+SYNOPSIS
+isreplace (sf, fit)
+
+.nf
+pointer	sf		# surface grid descriptor
+real	fit[ARB]	# fit array
+.fi
+
+.ih
+ARGUMENTS
+
+.ls sf
+The pointer to the surface grid descriptor.
+.le
+.ls fit
+The array containing the fit parameters and coefficients.
+.le
+
+.ih
+DESCRIPTION
+
+ISREPLACE restores the fit saved by ISSAVE for use by ISEVAL or ISVECTOR.
+.ih
+NOTES
+.ih
+SEE ALSO
+issave
+.endhelp
diff --git a/math/surfit/doc/isresolve.hlp b/math/surfit/doc/isresolve.hlp
new file mode 100644
index 00000000..afed8ed8
--- /dev/null
+++ b/math/surfit/doc/isresolve.hlp
@@ -0,0 +1,45 @@
+.help isresolve Aug85 "Surfit Package"
+.ih
+NAME
+isresolve -- resolve the surface, same lines array
+.ih
+SYNOPSIS
+include <math/surfit.h>
+
+isresolve (sf, lines, ier)
+
+.nf
+pointer	sf		# surface grid descriptor
+int	lines[nlines]	# surface grid line numbers, 1 <= lines[i] <= nlines
+int	ier		# error code
+.fi
+.ih
+ARGUMENTS
+.ls sf     
+Pointer to the surface grid descriptor structure.
+.le
+.ls nlines
+The number of surface grid lines to be fit.
+.le
+.ls ier      
+Error code returned by the fitting routines. The options are OK, SINGULAR,
+and NO_DEG_FREEDOM. If ier = SINGULAR the matrix is singular, ISSOLVE
+will compute a solution to the normal equations but one or more of the
+coefficients will be zero. If ier = NO_DEG_FREEDOM, too few data points
+exist for a reasonable solution to be computed. ISSOLVE returns without
+fitting the data.
+.le
+.ih
+DESCRIPTION
+In some applications it is necessary to refit the same surface several
+times with the same lines array. In this case it is inefficient to reaccumulate
+the matrices and calculate the Cholesky factorization for each fit.
+ISERESOLVE reaccumulates only the right side of the equation.
+.ih
+NOTES
+It is the responsibility of the calling program to check for out of bounds
+lines values.
+.ih
+SEE ALSO
+issolve, isinit, slzero, islfit, islrefit, islaccum, islresolve, isfree
+.endhelp
diff --git a/math/surfit/doc/issave.hlp b/math/surfit/doc/issave.hlp
new file mode 100644
index 00000000..6ec97cd6
--- /dev/null
+++ b/math/surfit/doc/issave.hlp
@@ -0,0 +1,55 @@
+.help issave Apr85 "Surfit Package"
+.ih
+NAME
+issave - save the surface fit
+.ih
+SYNOPSIS
+issave (sf, fit)
+
+.nf
+pointer	sf		# surface descriptor
+real	fit[ARB]	# fit array
+.fi
+.ih
+ARGUMENTS
+.ls sf   
+Pointer to the surface grid descriptor.
+.le
+.ls fit
+Array for storing the fit. The fit array must be at least SAVE_COEFF + ncoeff
+elements long.
+.le
+
+.ih
+DESCRIPTION
+ISSAVE stores the surface parameters and coefficient array
+in the 1-D reall array fit. Fit must be at least ncoeff + SAVE_COEFF
+elements long where ncoeff is defined as follows.
+
+.nf
+	SF_LEGENDRE:		nxcoeff = xorder
+				nycoeff = yorder
+				ncoeff = xorder * yorder        xterms = yes
+				ncoeff = nycoeff + nxcoeff - 1  xterms = no
+
+	SF_CHEBYSHEV:		nxcoeff = xorder
+				nycoeff = yorder
+				ncoeff = xorder * yorder        xterms = yes
+				ncoeff = nycoeff + nxcoeff - 1  xterms = no
+
+	SF_SPLINE3:		nxcoeff = xorder + 3
+				nycoeff = yorder + 3
+				ncoeff = (xorder + 3) * (yorder + 3)
+
+	SF_SPLINE1:		nxcoeff = xorder + 1
+				nycoeff = yorder + 1
+				ncoeff = (xorder + 1) * (yorder + 1)
+
+.fi
+
+.ih
+NOTES
+.ih
+SEE ALSO
+isreplace
+.endhelp
diff --git a/math/surfit/doc/issolve.hlp b/math/surfit/doc/issolve.hlp
new file mode 100644
index 00000000..312d69a3
--- /dev/null
+++ b/math/surfit/doc/issolve.hlp
@@ -0,0 +1,49 @@
+.help issolve Aug85 "Surfit Package"
+.ih
+NAME
+issolve -- solve the surface
+.ih
+SYNOPSIS
+include <math/surfit.h>
+
+issolve (sf, lines, nlines, ier)
+
+.nf
+pointer	sf		# surface grid descriptor
+int	lines[nlines]	# surface grid line numbers, 1 <= lines[i] <= nlines
+int	nlines		# number of lines
+int	ier		# error code
+.fi
+.ih
+ARGUMENTS
+.ls sf     
+Pointer to the surface grid descriptor structure.
+.le
+.ls lines
+The line number of surface grid lines to be included in the fit.
+.le
+.ls nlines
+The number of surface grid lines to be fit.
+.le
+.ls ier      
+Error code returned by the fitting routines. The options are OK, SINGULAR,
+and NO_DEG_FREEDOM. If ier = SINGULAR the matrix is singular, ISSOLVE
+will compute a solution to the normal equations but one or more of the
+coefficients will be zero. If ier = NO_DEG_FREEDOM, too few data points
+exist for a reasonable solution to be computed. ISSOLVE returns without
+fitting the data.
+.le
+.ih
+DESCRIPTION
+ISSOLVE solves for the surface coefficients.
+An error code is
+returned by ISSOLVE if it is unable to solve the normal equations as
+formulated.
+.ih
+NOTES
+It is the responsibility of the calling program to check for out of bounds
+lines values.
+.ih
+SEE ALSO
+isinit, slzero, islfit, islrefit, islaccum, islresolve, isfree
+.endhelp
diff --git a/math/surfit/doc/isvector.hlp b/math/surfit/doc/isvector.hlp
new file mode 100644
index 00000000..fabbdcc7
--- /dev/null
+++ b/math/surfit/doc/isvector.hlp
@@ -0,0 +1,41 @@
+.help isvector Aug85 "Surfit Package"
+.ih
+NAME
+isvector -- evaluate the fitted surface at a set of points
+.ih
+SYNOPSIS
+isvector (sf, x, y, zfit, npts)
+
+.nf
+pointer	sf		# surface grid descriptor
+real	x[npts]		# x array, 1 <= x[i] <= ncols
+real	y[npts]		# y array, 1 <= y[i] <= nlines
+real	zfit[npts]	# data values
+int	npts		# number of data points
+.fi
+.ih
+ARGUMENTS
+.ls sf    
+Pointer to the surface grid descriptor structure.
+.le
+.ls x, y
+Arrays of x and y values.
+.le
+.ls zfit
+Array of fitted values.
+.le
+.ls npts
+The number of points to be fit.
+.le
+.ih
+DESCRIPTION
+Fit the surface to an array of x and y values. ISVECTOR uses the coefficients
+stored in the surface descriptor structure.
+.ih
+NOTES
+Checking for out of bounds x and y values is the responsibility of the
+calling program.
+.ih
+SEE ALSO
+iseval
+.endhelp
diff --git a/math/surfit/doc/iszero.hlp b/math/surfit/doc/iszero.hlp
new file mode 100644
index 00000000..4c338807
--- /dev/null
+++ b/math/surfit/doc/iszero.hlp
@@ -0,0 +1,26 @@
+.help iszero Aug85 "Surfit Package"
+.ih
+NAME
+iszero -- set up for a new surface fit
+.ih
+SYNOPSIS
+iszero (sf)
+
+.nf
+pointer	sf		# surface grid descriptor
+.fi
+.ih
+ARGUMENTS
+.ls sf     
+Pointer to the surface grid descriptor structure.
+.le
+.ih
+DESCRIPTION
+ISZERO zeros all matrices and vectors used by the fitting program. It
+should be used only if one wishes to refit the surface from scratch.
+.ih
+NOTES
+.ih
+SEE ALSO
+isinit, islzero, islaccum, islfit, islrefit, islsolve, issolve, isfree
+.endhelp
diff --git a/math/surfit/doc/surfit.hd b/math/surfit/doc/surfit.hd
new file mode 100644
index 00000000..18e7ec1a
--- /dev/null
+++ b/math/surfit/doc/surfit.hd
@@ -0,0 +1,19 @@
+# HELP Directory for SURFIT package
+
+$surfit = "math$surfit/"
+
+iscoeff		hlp = iscoeff.hlp,		src = surfit$iscoeff.x
+iseval		hlp = iseval.hlp,		src = surfit$iseval.x
+isfree		hlp = isfree.hlp,		src = surfit$isfree.x
+isinit		hlp = isinit.hlp,		src = surfit$isinit.x
+islaccum	hlp = islaccum.hlp,		src = surfit$islaccum.x
+islfit		hlp = islfit.hlp,		src = surfit$islfit.x
+islrefit	hlp = islrefit.hlp,		src = surfit$islrefit.x
+islsolve	hlp = islsolve.hlp,		src = surfit$islsolve.x
+islzero		hlp = islzero.hlp,		src = surfit$islzero.x
+isreplace	hlp = isreplace.hlp,		src = surfit$isreplace.x
+isresolve	hlp = isresolve.hlp,		src = surfit$isresolve.x
+issave		hlp = issave.hlp,		src = surfit$issave.x
+issolve		hlp = issolve.hlp,		src = surfit$issolve.x
+isvector	hlp = isvector.hlp,		src = surfit$isvector.x
+iszero		hlp = iszero.hlp,		src = surfit$iszero.x
diff --git a/math/surfit/doc/surfit.hlp b/math/surfit/doc/surfit.hlp
new file mode 100644
index 00000000..a7c347cc
--- /dev/null
+++ b/math/surfit/doc/surfit.hlp
@@ -0,0 +1,157 @@
+.help surfit Aug84 "Math Package"
+
+.ih
+NAME
+surfit -- set of routines for fitting gridded surface data
+
+
+.ih
+SYNOPSIS
+
+The surface is defined in the following way. The x and y values are
+assumed to be integers and lie in the range 0 <= x <= ncols + 1 and
+0 <= y <= nlines + 1. The package prefix is for image surface fit.
+Package routines with a prefix followed by l operate on individual
+image lines only.
+
+.nf
+        isinit	(sf, surf_type, xorder, yorder, xterms, ncols, nlines)
+	iszero	(sf)
+       islzero	(sf, lineno)
+      islaccum	(sf, cols, lineno, z, w, ncols, wtflag)
+      islsolve	(sf, lineno, ier)
+        islfit	(sf, cols, lineno, z, w, ncols, wtflag, ier)
+      islrefit	(sf, cols, lineno, z, w, ier)
+       issolve	(sf, lines, nlines, ier)
+     isresolve	(sf, lines, ier)
+    z = iseval	(sf, x, y)
+      isvector	(sf, x, y, zfit, npts)
+        issave	(sf, surface)
+     isrestore	(sf, surface)
+        isfree	(sf)
+.fi
+
+.ih
+DESCRIPTION
+
+The SURFIT package provides a set of routines for fitting surfaces to functions
+linear in their coefficients using the tensor product method and least squares
+techniques. The basic numerical
+technique employed is the solution of the normal equations by the Cholesky
+method.
+
+.ih
+NOTES
+
+The following series of steps illustrates the use of the package.
+
+.ls 4
+.ls (1)
+Include an include <math/surfit.h> statement in the calling program to make the
+SURFIT package definitions available to the user program.
+.le
+.ls (2)
+Call ISINIT to initialize the surface fitting parameters.
+.le
+.ls (3)
+Call ISLACCUM to select a weighting function and accumulate data points
+for each line into the appropriate arrays and vectors. Call ISLSOLVE
+to compute the coefficients in x for each line. The coefficients are
+stored inside SURFIT. Repeat this procedure for each line. ISLACCUM
+and SFLSOLVE can be combined by a call to SFLFIT. If the x values
+and weights remain the same from line to line ISLREFIT can be called.
+.le
+.ls (4)
+Call ISSOLVE to solve for the surface coefficients.
+.le
+.ls (5)
+Call ISEVAL or ISVECTOR to evaluate the fitted surface at the x and y
+value(s) of interest.
+.le
+.ls (6)
+Call ISFREE to release the space allocated for the fit.
+.le
+.le
+
+.ih
+EXAMPLES
+
+.nf
+Example 1: Fit a 2nd order Lengendre polynomial in x and y to an image
+and output the fitted image
+
+include <imhdr.h>
+include <math/surfit.h>
+
+ 	
+	old = immap (old_image, READ_ONLY, 0)
+	new = immap (new_image, NEW_COPY, 0)
+	
+	ncols = IM_LEN(old, 1)
+	nlines = IM_LEN(old, 2)
+	
+	# initialize surface fit
+	call isinit (sf, LEGENDRE, 3, 3, YES, ncols, nlines)
+        
+	# allocate space for lines, columns and weight arrays
+	call malloc (cols, ncols, TY_INT)
+	call malloc (lines, nlines, TY_INT)
+	call malloc (weight, ncols, TY_REAL)
+
+	# initialize lines and columns arrays
+	do i = 1, ncols
+	    Memi[cols - 1 + i] = i
+	do i = 1, nlines
+	    Memi[lines - 1 + i] = i
+
+	# fit the surface in x line by line
+	call amovkl (long(1), v, IM_MAXDIM)
+	do i = 1, nlines {
+	    if (imgnlr (old, inbuf, v) == EOF)
+		call error (0, "Error reading image.")
+	    if (i == 1) {
+		call islfit (sf, Memi[cols], i, Memr[inbuf], Memr[weight],
+				ncols, SF_WTSUNIFORM, ier)
+		if (ier != OK)
+		    ...
+	    } else
+		call islrefit (sf, Memi[cols], i, Memr[inbuf], Memr[weight],
+				ier)
+	}
+
+	# solve for surface coefficients
+	call issolve (sf, Memi[lines], nlines, ier)
+
+	# free space used in fitting arrays
+	call mfree (cols, TY_INT)
+	call mfree (lines, TY_INT)
+	call mfree (weight, TY_REAL)
+
+	# allocate space for x and y arrays
+	call malloc (x, ncols, TY_REAL)
+	call malloc (y, ncols, TY_REAL)
+
+	# intialize z array
+	do i = 1, ncols
+	    Memr[x - 1 + i] = real (i)
+
+	# create fitted image
+	call amovkl (long(10, v, IM_MAXDIM)
+	do i = 1, nlines {
+	    if (impnlr (new, outbuf, v) == EOF)
+		call error (0, "Error writing image.")
+	    call amovkr (real (i), Memr[y], ncols)
+	    call isvector (sf, Memr[x], Memr[y], Memr[outbuf], ncols)
+	}
+
+	# close files and cleanup
+	call mfree (x, TY_REAL)
+	call mfree (y, TY_REAL
+
+	call isfree (sf)
+
+	call imunmap (old)
+	call imunmap (new)
+
+.fi
+.endhelp
diff --git a/math/surfit/doc/surfit.men b/math/surfit/doc/surfit.men
new file mode 100644
index 00000000..4a3b3258
--- /dev/null
+++ b/math/surfit/doc/surfit.men
@@ -0,0 +1,15 @@
+         isinit - initialize surface descriptor
+        iscoeff - get coefficients of surface
+         iseval - evaluate the surface at a point
+         isfree - free surface descriptor
+       islaccum - accumulate surface at a given line number
+         islfit - fit the surface at a given line number
+       islrefit - refit the surface at a given line number, same columns
+       islsolve - solve the surface at a given line number
+        islzero - zero the fit for a given line number
+        issolve - solve the surface
+      isreplace - restore the surface fit
+      isresolve - resolve surface, same lines
+         issave - save the surface fit
+       isvector - fit surface at a set of (x,y) values
+         iszero - zero the surface fit
diff --git a/math/surfit/doc/surfit.spc b/math/surfit/doc/surfit.spc
new file mode 100644
index 00000000..07f2b934
--- /dev/null
+++ b/math/surfit/doc/surfit.spc
@@ -0,0 +1,500 @@
+.help surfit Aug84 "Math Package"
+.CE
+Specifications for the Surface Fitting Package
+.CE
+Lindsey Davis
+.CE
+August 1984
+.CE
+First Draft
+
+.sh
+1. Introduction
+
+The SURFIT package provides a set of routines for fitting surfaces to functions
+linear in their coefficients using least squares techniques. The basic numerical
+technique employed is the solution of the normal equations by the Cholesky
+method.
+
+.sh
+2. Requirements
+
+.ls (1)
+The package shall take as input a surface or section of
+a surface.
+The data are assumed to be on a rectangular grid and to have the following
+ranges, 0 <= x <= ncols + 1 and 0 <= y <= nlines + 1.
+The package routines assume that data values equal to INDEF have
+been removed from the data set or replaced with appropriate interpolated
+values prior to entering the package routines. It is not necessary that
+all the data be present to perform the fit.
+.le
+.ls (2)
+The package shall perform the following operations:
+.ls o
+Determine the coefficients of the fitting function by solving the normal
+equations. The surface fitting function is selected at run time from the
+following list: (1) SF_LEGENDRE, Legendre polynomials in x and y,
+(2) SF_CHEBYSHEV, Chebyshev polynomials in x and y, (3) SF_SPLINE3, bicubic
+spline with break points evenly spaced in x and y. The calling sequence
+must be invariant to the form of the fitting function.
+.le
+.ls o
+Set an error code if the numerical routines are unable to fit the specified
+function.
+.le
+.ls o
+Optionally output the values of the coefficients. The coefficients will be
+stored internal to the SURFIT package. However in some applications it is
+the coefficients which are of primary interest. A package routine shall
+exist to extract the coefficients from the curve descriptor structure.
+.le
+.ls o
+Evaluate the surface at arbitrary value(s) of x and y. The evaluating
+routines shall use the calculated coefficients and the user supplied
+x values(s).
+.le
+.ls o
+Be capable of storing the fitted surface parameters and coefficients in a
+user supplied array and restoring the saved fit for later use by the
+evaluating routines.
+.le
+.ls o
+Input surface parameters and surface coefficients derived external to
+the surfit package into the surfit format for use by the evaluate
+routines.
+.le
+.le
+.ls (3)
+The program shall perform a weighted fit to the surface using a user
+supplied weight array and weight flag. The weighting options will be
+SF_WTSUSER and SF_WTSUNIFORM. In SF_WTSUNIFORM mode the package routines
+set all the weights to 1. In SF_WTSUSER mode the package routines apply
+user supplied weights to the individual data points.
+.le
+.ls (4)
+The input data set, output coefficients, error and fitted z arrays are
+single precision real quantities. All package arithmetic shall be performed
+in single precision.
+.le
+
+.sh
+3. Specifications
+
+.sh
+3.1. List of Routines
+
+
+The surface is defined in the following way. The x and y values are
+assumed to be integers and lie in the range 0 <= x <= ncols + 1 and
+0 <= y <= nlines + 1. The package prefix is is for image surface fit.
+Package routines with a prefix followed by l operate on individual
+image lines only.
+
+.nf
+	z = f (col, line)
+.fi
+
+.nf
+        isinit	(sf, surf_type, xorder, yorder, xterms, ncols, nlines)
+	iszero	(sf)
+       islzero	(sf, lineno)
+      islaccum	(sf, cols, lineno, z, w, ncols, wtflag)
+      islsolve	(sf, lineno, ier)
+        islfit	(sf, cols, lineno, z, w, ncols, wtflag, ier)
+      islrefit	(sf, cols, lineno, z, w, ier)
+       issolve	(sf, lines, nlines, ier)
+     isresolve	(sf, lines, ier)
+    z = iseval	(sf, x, y)
+      isvector	(sf, x, y, zfit, npts)
+        issave	(sf, surface)
+     isrestore	(sf, surface)
+        isfree	(sf)
+.fi
+
+.sh
+3.2. Algorithms
+
+.sh
+3.2.1. Polynomial Basis Functions
+
+The approximating function is assumed to be of the following form,
+
+.nf
+    f(x,y) = sum (a[i,j] * F(i,x) * F(j,y))    i = 1...nxcoeff  j = 1...nycoeff
+.fi
+
+where the F(i,x) are the polynomial basis functions containing terms of order
+x**(i-1) and the a(i,j) are the coefficients. In order to avoid a very
+ill-conditioned linear system for moderate or large i, the Legendre and
+Chebyshev polynomials were chosen for the basis functions. The Chebyshev
+and Legendre polynomials are orthogonal over -1. <= x,y <= 1. The data x and
+y values are normalized to this regime using the number of lines and columns
+supplied by the user. For each data point the basis
+functions are calculated using the following recursion relations. The
+cross terms are optional.
+
+Legendre series:
+.nf
+
+    F(1,x) = 1.
+    F(2,x) = x
+    F(i,x) = [(2*i-1)*x*F(i-1,x)-(i-2)*F(i-2,x)]/(i-1)
+    F(1,y) = 1.
+    F(2,y) = y
+    F(j,y) = [(2*j-1)*y*F(j-1,y)-(j-2)*F(j-2,y)]/(j-1)
+
+.fi
+Chebyshev series:
+.nf
+
+    F(1,x) = 1.
+    F(2,x) = x
+    F(i,x) = 2.*x*F(i-1,x)-F(i-2,x)
+    F(1,y) = 1.
+    F(2,y) = y
+    F(j,y) = 2.*y*F(j-1,y)-F(j-2,y)
+.fi
+
+
+.sh
+3.2.2. Bicubic Cardinal B-spline
+
+The approximating function is of the form
+
+.nf
+    f(x,y) =  sum (x(i,j) * F(i,x) * F(j,y))   i=1...nxcoeff j=1...nycoeff
+.fi
+
+where the basis functions, F(i,x), are the cubic cardinal B-splines
+(Prenter 1975). The user supplies the number of columns and lines and the
+number of polynomial pieces in x and y, nxpieces and nypieces, to be fit
+to the data set. The number of bicubic spline coefficients, ncoeff, will be
+
+.nf
+    nxcoeff = (nxpieces + 3)
+    nycoeff = (nypieces + 3)
+    ncoeff = nxcoeff * nycoeff
+.fi
+
+The cardinal B-spline is stored in a lookup table. For each x and y the
+appropriate break point is selected and the four non-zero B-splines are
+calculated by nearest neighbour interpolation in the lookup table.
+
+.sh
+3.2.3. Method of Solution
+
+.sh
+3.2.3.1. Full Surface Fit
+
+The normal equations are accumulated in the following form.
+
+.nf
+    c * x = b
+    c[i,j] = (B(i,x,y), B(j,x,y))
+    b[j] = (B(j,x,y), S(x,y))
+.fi
+
+B(i,x,y) is the ith basis function at x and y, S(x,y) is the surface to
+be approximated and the inner product of two functions G and H is
+given by 
+
+.nf
+    (G,H) = sum (w[i] * G(x[i],y[i]) * H(x[i],y[i]))  i = 1...npts
+.fi
+
+Since the matrix c is symmetric and positive semi-definite it may
+be solved by the Cholesky method.
+The coefficient matrix c  can be written as
+
+.nf
+    c = l * d * l-transpose
+.fi
+
+where l is a unit lower triangular matrix and d is the diagonal of c
+(de Boor 1978). Near zero pivots are handled in the following way.
+At the nth elimination step the current value of the nth diagonal
+element is compared with the original nth diagonal element. If the
+diagonal element has been reduced by one computer word length, the
+entire nth row is declared linearly dependent on the previous n-1
+rows and x(n) = 0.
+
+The triangular system
+
+.nf
+    l * w = b
+.fi
+
+is solved for w (forward substitution), the vector d  ** (-1) * w is
+computed and the triangular system
+
+.nf
+    l-transpose * x = d ** (-1) * w
+.fi
+
+solved for the coefficients, x (backward substitution).
+
+
+.sh
+3.2.3.2. Line by Line Fit
+
+Each line of the surface can be represented by the following equation.
+
+.nf
+    S(x,yo) = a[i](yo) * B[i](x)            i = 1...nxcoeff for each yo
+.fi
+
+The normal equations for each image line are formed
+
+.nf
+    c * a = b
+    c[i,j] = (B(i,x), B(j,x))
+    b[j] = (B(j,x), S(x,yo))
+.fi
+
+and solved for a as described in the previous section.
+After fitting the entire image matrix a has nxcoeff columns and
+nlines rows.
+
+For each column i in a the normal equations
+are formed and solved for the c[i,j]
+
+.nf
+    c * x  = b
+    c[j,k] = (B(j,y), B(k,x))
+    b[i,j] = (a[i](y), B(j,y))
+.fi
+
+.sh
+3.2.4. Number of Operations
+
+It is worth while to calculate the number of operations reqired to compute
+the surface fit by the full surface method versus the line by line method.
+The number of operations required to accumulate the normal equations
+and solve the set is the following
+
+.nf
+    nop = npts * order ** 2 / 2 + order ** 3 /6
+.fi
+
+where the first term is the number of operations required to accumulate the
+matrix and the second is the number required to solve the system.
+
+.nf
+Example 1: full surface, 3 by 3 polynomial no cross terms, 512 by 512 image
+
+    norder = 5
+    npts = 512 * 512
+    nops(accum) = 3,276,800
+    nops(solve) = 21
+
+Example 2: full surface, 30 by 30 spline, 512 by 512 image
+
+    norder = 900
+    npts = 512 * 512
+    nops(accum) = 1.062 * 10 ** 11
+    nops(solve) = 1.216 * 10 ** 8
+
+Example 3: line by line method, 3 by 3 polynomial, 512 by 512 image
+
+    step 1 solve for a coefficients
+    norder = 3
+    npts = 512
+    nlines = 512
+    nops(accum) = 1,179,648
+    nops(solve) = 2304
+
+    step 2 solve for c coefficients
+    norder = 3
+    npts = 512
+    nlines = 3
+    nops(accum) = 6912
+    nops(solve) = 14
+
+Example 4: line by line method, 30 by 30 spline, 512 by 512 image
+
+    step 1 solve coefficients
+    norder = 30
+    npts = 512
+    nlines = 512
+    nops(accum) = 117,964,800
+    nops(solve) = 2,304,000
+
+    step 2 solve for c coefficients
+    norder = 30
+    npts = 512
+    nlines = 30
+    nops(accum) = 6,912,000
+    nops(solve) = 135,000
+.fi
+
+The figures for the line by line method are worst case numbers. If the
+x values remain the same from line to line then the coefficient matrix
+only has to be accumulated and inverted twice.
+For the bicubic spline function the number of operations is significantly
+reduced by taking advantage of the banded nature of the matrices.
+
+.sh
+4. Usage
+
+.sh
+4.1. User Notes
+
+The following series of steps illustrates the use of the package.
+
+.ls 4
+.ls (1)
+Include an include <surfit.h> statement in the calling program to make the
+SURFIT package definitions available to the user program.
+.le
+.ls (2)
+Call SFINIT to initialize the surface fitting parameters.
+.le
+.ls (3)
+Call SFLACCUM to select a weighting function and accumulate data points
+for each line into the appropriate arrays and vectors. Call SFLSOLVE
+to compute the coefficients in x for each line. The coefficents are
+stored inside SURFIT. Repeat this procedure for each line. SFLACCUM
+and SFLSOLVE can be combined by a call to SFLFIT. If the x values
+and weights remain the same from line to line SFLREFIT can be called.
+.le
+.ls (4)
+Call SFSOLVE to solve for the surface coefficients.
+.le
+.ls (5)
+Call SFEVAL or SFVECTOR to evaluate the fitted surface at the x and y
+value(s) of interest.
+.le
+.ls (6)
+Call SFFREE to release the space allocated for the fit.
+.le
+.le
+
+.sh
+4.2. Examples
+
+.nf
+Example 1: Fit a 2nd order Lengendre polynomial in x and y to an image
+and output the fitted image
+
+include <imhdr.h>
+include <math/surfit.h>
+
+ 	
+	old = immap (old_image, READ_ONLY, 0)
+	new = immap (new_image, NEW_COPY, 0)
+	
+	ncols = IM_LEN(old, 1)
+	nlines = IM_LEN(old, 2)
+	
+	# initialize surface fit
+	call isinit (sf, LEGENDRE, 3, 3, YES, ncols, nlines)
+        
+	# allocate space for lines, columns and weight arrays
+	call malloc (cols, ncols, TY_INT)
+	call malloc (lines, nlines, TY_INT)
+	call malloc (weight, ncols, TY_REAL)
+
+	# initialize lines and columns arrays
+	do i = 1, ncols
+	    Memi[cols - 1 + i] = i
+	do i = 1, nlines
+	    Memi[lines - 1 + i] = i
+
+	# fit the surface in x line by line
+	call amovkl (long(1), v, IM_MAXDIM)
+	do i = 1, nlines {
+	    if (imgnlr (old, inbuf, v) == EOF)
+		call error (0, "Error reading image.")
+	    if (i == 1) {
+		call islfit (sf, Memi[cols], i, Memr[inbuf], Memr[weight],
+				ncols, SF_WTSUNIFORM, ier)
+		if (ier != OK)
+		    ...
+	    } else
+		call sflrefit (sf, Memi[cols], i, Memr[inbuf], Memr[weight],
+				ier)
+	}
+
+	# solve for surface coefficients
+	call issolve (sf, Memi[lines], nlines, ier)
+
+	# free space used in fitting arrays
+	call mfree (cols, TY_INT)
+	call mfree (lines, TY_INT)
+	call mfree (weight, TY_REAL)
+
+	# allocate space for x and y arrays
+	call malloc (x, ncols, TY_REAL)
+	call malloc (y, ncols, TY_REAL)
+
+	# intialize z array
+	do i = 1, ncols
+	    Memr[x - 1 + i] = real (i)
+
+	# create fitted image
+	call amovkl (long(10, v, IM_MAXDIM)
+	do i = 1, nlines {
+	    if (impnlr (new, outbuf, v) == EOF)
+		call error (0, "Error writing image.")
+	    call amovkr (real (i), Memr[y], ncols)
+	    call isvector (sf, Memr[x], Memr[y], Memr[outbuf], ncols)
+	}
+
+	# close files and cleanup
+	call mfree (x, TY_REAL)
+	call mfree (y, TY_REAL
+
+	call isfree (sf)
+
+	call imunmap (old)
+	call imunmap (new)
+
+.fi
+.sh
+5. Detailed Design
+
+.sh
+5.1. Surface Descriptor Structure
+
+To be written when specifications are finalised.
+
+.sh
+5.2. Storage Requirements
+
+The minimum storage requirements for the full surface fit method
+assuming that points are to be rejected from the fit without revaluating
+the matrix are the following.
+
+.nf
+    (nxcoeff*nycoeff) ** 2    -- coefficient array plus triangularized matrix
+    nxcoeff*nycoeff           -- right side
+    nxcoeff*nycoeff           -- solution vector
+.fi
+
+For a 30 by 30 spline roughly 811,800 storage units are required.
+For a 3 by 3 polynomial the requirements are 675 storage units. The
+requirements are roughly half when rejection is disabled.
+
+The minimum storage requirements for the line by line method are
+
+.nf
+    nx*nx*2   -- The x matrix and its factorization
+    nx*nlines -- The x coefficient matrix
+    ny*ny*2 -- The y matrix and its factorization
+    nx*ny -- The solution
+.fi
+
+
+.sh
+6. References
+
+.ls (1)
+Carl de Boor, "A Practical Guide to Splines", 1978, Springer-Verlag New York
+Inc.
+.le
+.ls (2)
+P.M. Prenter, "Splines and Variational Methods", 1975, John Wiley and Sons
+Inc.
+.le
+.endhelp
diff --git a/math/surfit/iscoeff.x b/math/surfit/iscoeff.x
new file mode 100644
index 00000000..f656e7e6
--- /dev/null
+++ b/math/surfit/iscoeff.x
@@ -0,0 +1,37 @@
+# Copyright(c) 1986 Association of Universities for Research in Astronomy Inc.
+
+include "surfitdef.h"
+
+# ISCOEFF -- Procedure to fetch the number and magnitude of the coefficients.
+# If the surface type is a bicubic spline or polynomials with cross terms
+# then the  number of coefficients is SF_NXCOEFF(sf) * SF_NYCOEFF(sf) and
+# the coefficient of B(i,x) * B(j,y) will be stored in element
+# (i - 1) * SF_NYCOEFF(sf) + j of the array coeff. Otherwise the number
+# of coefficients will be (SF_NXCOEFF(sf) + SF_NYCOEFF(sf) - 1), and
+# the SF_NYCOEFF(sf) y coefficients will be output first followed by
+# the (SF_NXCOEFF(sf) - 1) x coefficients.
+
+procedure iscoeff (sf, coeff, ncoeff)
+
+pointer	sf		# pointer to the surface fitting descriptor
+real	coeff[ARB]	# the coefficients of the fit
+int	ncoeff		# the number of coefficients
+
+int	i
+pointer	cptr
+
+begin
+	# calculate the number of coefficients
+	if (SF_XTERMS(sf) == NO) {
+	    ncoeff = SF_NXCOEFF(sf) + SF_NYCOEFF(sf) - 1
+	    call amovr (COEFF(SF_COEFF(sf)), coeff, SF_NYCOEFF(sf))
+	    cptr = SF_COEFF(sf) + SF_NYCOEFF(sf)
+	    do i = SF_NYCOEFF(sf) + 1, ncoeff {
+		coeff[i] = COEFF(cptr)
+		cptr = cptr + SF_NYCOEFF(sf)
+	    }
+	} else {
+	    ncoeff = SF_NXCOEFF(sf) * SF_NYCOEFF(sf)
+	    call amovr (COEFF(SF_COEFF(sf)), coeff, ncoeff)
+	}
+end
diff --git a/math/surfit/iseval.x b/math/surfit/iseval.x
new file mode 100644
index 00000000..08ef5f89
--- /dev/null
+++ b/math/surfit/iseval.x
@@ -0,0 +1,92 @@
+# Copyright(c) 1986 Association of Universities for Research in Astronomy Inc.
+
+include <math/surfit.h>
+include "surfitdef.h"
+
+# ISEVAL -- Procedure to evaluate the fitted surface at a single point.
+# The SF_NXCOEFF(sf) by SF_NYCOEFF(sf) coefficients are stored in the
+# SF_NYCOEFF(sf) by SF_NXCOEFF(sf) matrix COEFF. The j-th element of the ith
+# row of COEFF contains the coefficient of the i-th basis function in x and
+# the j-th basis function in y.
+
+real procedure iseval (sf, x, y)
+
+pointer	sf		# pointer to surface descriptor structure
+real	x		# x value
+real	y		# y value
+
+real	sum, accum
+int	i, k, leftx, lefty, yorder
+pointer	sp, xb, xzb, yb, yzb, czptr
+
+begin
+	# allocate space for the basis functions
+	call smark (sp)
+	call salloc (xb, SF_XORDER(sf), MEM_TYPE)
+	xzb = xb - 1
+	call salloc (yb, SF_YORDER(sf), MEM_TYPE)
+	yzb = yb - 1
+
+	# calculate the basis functions
+	switch (SF_TYPE(sf)) {
+	case SF_CHEBYSHEV:
+	    leftx = 0
+	    lefty = 0
+	    czptr = SF_COEFF(sf) - 1
+	    call sf_b1cheb (x, SF_XORDER(sf), SF_XMAXMIN(sf), SF_XRANGE(sf),
+		XBS(xb))
+	    call sf_b1cheb (y, SF_YORDER(sf), SF_YMAXMIN(sf), SF_YRANGE(sf),
+		YBS(yb))
+
+	case SF_LEGENDRE:
+	    leftx = 0
+	    lefty = 0
+	    czptr = SF_COEFF(sf) - 1
+	    call sf_b1leg (x, SF_XORDER(sf), SF_XMAXMIN(sf), SF_XRANGE(sf),
+		XBS(xb))
+	    call sf_b1leg (y, SF_YORDER(sf), SF_YMAXMIN(sf), SF_YRANGE(sf),
+		YBS(yb))
+
+	case SF_SPLINE3:
+	    call sf_b1spline3 (x, SF_NXPIECES(sf), -SF_XMIN(sf),
+	        SF_XSPACING(sf), XBS(xb), leftx)
+	    call sf_b1spline3 (y, SF_NYPIECES(sf), -SF_YMIN(sf),
+	        SF_YSPACING(sf), YBS(yb), lefty)
+	    czptr = SF_COEFF(sf) - 1 + lefty + leftx * SF_NYCOEFF(sf)
+
+	case SF_SPLINE1:
+	    call sf_b1spline1 (x, SF_NXPIECES(sf), -SF_XMIN(sf),
+	        SF_XSPACING(sf), XBS(xb), leftx)
+	    call sf_b1spline1 (y, SF_NYPIECES(sf), -SF_YMIN(sf),
+	        SF_YSPACING(sf), YBS(yb), lefty)
+	    czptr = SF_COEFF(sf) - 1 + lefty + leftx * SF_NYCOEFF(sf)
+	}
+
+	# initialize accumulator
+	# basis functions
+	sum = 0.
+
+	# loop over y basis functions
+	yorder = SF_YORDER(sf)
+	do i = 1, SF_XORDER(sf) {
+
+	    # loop over the x basis functions
+	    accum = 0.
+	    do k = 1, yorder {
+		accum = accum + COEFF(czptr+k) * YBS(yzb+k) 
+	    }
+	    accum = accum * XBS(xzb+i)
+	    sum = sum + accum
+
+	    # elements of COEFF where neither k = 1 or i = 1
+	    # are not calculated if SF_XTERMS(sf) = NO
+	    if (SF_XTERMS(sf) == NO)
+		yorder = 1
+
+	    czptr = czptr + SF_NYCOEFF(sf)
+	}
+
+	call sfree (sp)
+
+	return (sum)
+end
diff --git a/math/surfit/isfree.x b/math/surfit/isfree.x
new file mode 100644
index 00000000..2eb25891
--- /dev/null
+++ b/math/surfit/isfree.x
@@ -0,0 +1,45 @@
+# Copyright(c) 1986 Association of Universities for Research in Astronomy Inc.
+
+include "surfitdef.h"
+
+# ISFREE -- Procedure to free the surface descriptor
+
+procedure isfree (sf)
+
+pointer	sf		# pointer to the surface descriptor
+errchk	mfree
+
+begin
+	# free arrays in memory
+
+	# first the basis functions
+	if (SF_XBASIS(sf) != NULL)
+	    call mfree (SF_XBASIS(sf), MEM_TYPE)
+	if (SF_XLEFT(sf) != NULL)
+	    call mfree (SF_XLEFT(sf), TY_INT)
+	if (SF_YBASIS(sf) != NULL)
+	    call mfree (SF_YBASIS(sf), MEM_TYPE)
+	if (SF_YLEFT(sf) != NULL)
+	    call mfree (SF_YLEFT(sf), TY_INT)
+
+	# next the x and y matrices
+	if (SF_XMATRIX(sf) != NULL)
+	    call mfree (SF_XMATRIX(sf), MEM_TYPE)
+	if (SF_YMATRIX(sf) != NULL)
+	    call mfree (SF_YMATRIX(sf), MEM_TYPE)
+
+	# last the coefficient matrices
+	if (SF_XCOEFF(sf) != NULL)
+	    call mfree (SF_XCOEFF(sf), MEM_TYPE)
+	if (SF_COEFF(sf) != NULL)
+	    call mfree (SF_COEFF(sf), MEM_TYPE)
+
+	if (SF_WZ(sf) != NULL)
+	    call mfree (SF_WZ(sf), MEM_TYPE)
+	if (SF_TLEFT(sf) != NULL)
+	    call mfree (SF_TLEFT(sf), TY_INT)
+
+	# free surface descriptor
+	if (sf != NULL)
+	    call mfree (sf, TY_STRUCT)
+end
diff --git a/math/surfit/isinit.x b/math/surfit/isinit.x
new file mode 100644
index 00000000..ade35b47
--- /dev/null
+++ b/math/surfit/isinit.x
@@ -0,0 +1,167 @@
+# Copyright(c) 1986 Association of Universities for Research in Astronomy Inc.
+
+include <math/surfit.h>
+include "surfitdef.h"
+
+# ISINIT --  Procedure to set up the surface  descriptor.
+
+procedure isinit (sf, surf_type, xorder, yorder, xterms, ncols, nlines)
+
+pointer	sf		# pointer to surface descriptor structure
+int	surf_type	# type of surface to be fitted
+int	xorder		# x order of surface to be fit, or in the case of the
+			# spline the number of polynomial pieces in x to be fit
+int	yorder		# y order of surface to be fit, or in the case of the
+			# spline the number of polynomial pieces in y to be fit
+int	xterms		# cross terms for polynomials?
+int	ncols		# number of columns in the surface
+int	nlines		# number of lines in the surface
+
+int	i
+pointer	x, y
+pointer	sp
+errchk	malloc, calloc
+
+begin
+	# allocate space for the surface descriptor
+	call malloc (sf, LEN_SFSTRUCT, TY_STRUCT)
+
+	if (xorder < 1 || yorder < 1)
+	    call error (0, "SFLINIT: Illegal order.")
+
+	if (ncols < 1)
+	    call error (0, "SFLINIT: x data range is 0.")
+	if (nlines < 1)
+	    call error (0, "SFLINIT: y data range is 0.")
+
+	switch (surf_type) {
+
+	case SF_CHEBYSHEV, SF_LEGENDRE:
+	    SF_NXCOEFF(sf) = xorder
+	    SF_XORDER(sf) = xorder
+	    SF_XRANGE(sf) = 2. / real (ncols + 1)
+	    SF_XMAXMIN(sf) = - real (ncols + 1) / 2.
+	    SF_XMIN(sf) = 0.
+	    SF_XMAX(sf) = real (ncols + 1)
+	    SF_NYCOEFF(sf) = yorder
+	    SF_YORDER(sf) = yorder
+	    SF_YRANGE(sf) = 2. / real (nlines + 1)
+	    SF_YMAXMIN(sf) = - real (nlines + 1) / 2.
+	    SF_YMIN(sf) = 0.
+	    SF_YMAX(sf) = real (nlines + 1)
+	    SF_XTERMS(sf) = xterms
+
+	case SF_SPLINE3:
+	    SF_NXCOEFF(sf) = (xorder + SPLINE3_ORDER - 1)
+	    SF_XORDER(sf) = SPLINE3_ORDER
+	    SF_NXPIECES(sf) = xorder - 1
+	    SF_XSPACING(sf) = xorder / real (ncols + 1)
+	    SF_NYCOEFF(sf) = (yorder + SPLINE3_ORDER - 1)
+	    SF_YORDER(sf) = SPLINE3_ORDER
+	    SF_NYPIECES(sf) = yorder - 1
+	    SF_YSPACING(sf) = yorder / real (nlines + 1)
+	    SF_XMIN(sf) = 0.
+	    SF_XMAX(sf) = real (ncols + 1)
+	    SF_YMIN(sf) = 0.
+	    SF_YMAX(sf) = real (nlines + 1)
+	    SF_XTERMS(sf) = YES
+
+	case SF_SPLINE1:
+	    SF_NXCOEFF(sf) = (xorder + SPLINE1_ORDER - 1)
+	    SF_XORDER(sf) = SPLINE1_ORDER
+	    SF_NXPIECES(sf) = xorder - 1
+	    SF_XSPACING(sf) = xorder / real (ncols + 1)
+	    SF_NYCOEFF(sf) = (yorder + SPLINE1_ORDER - 1)
+	    SF_YORDER(sf) = SPLINE1_ORDER
+	    SF_NYPIECES(sf) = yorder - 1
+	    SF_YSPACING(sf) = yorder / real (nlines + 1)
+	    SF_XMIN(sf) = 0.
+	    SF_XMAX(sf) = real (ncols + 1)
+	    SF_YMIN(sf) = 0.
+	    SF_YMAX(sf) = real (nlines + 1)
+	    SF_XTERMS(sf) = YES
+
+	default:
+	    call error (0, "SFINIT: Unknown surface type.")
+	}
+
+	SF_TYPE(sf) = surf_type
+	SF_NLINES(sf) = nlines
+	SF_NCOLS(sf) = ncols
+
+	# allocate space for the matrix and vectors
+	call calloc (SF_XBASIS(sf), SF_XORDER(sf) * SF_NCOLS(sf),
+		    MEM_TYPE)
+	call calloc (SF_YBASIS(sf), SF_YORDER(sf) * SF_NLINES(sf),
+		    MEM_TYPE)
+	call calloc (SF_XMATRIX(sf), SF_XORDER(sf) * SF_NXCOEFF(sf), MEM_TYPE)
+	call calloc (SF_XCOEFF(sf), SF_NLINES(sf) * SF_NXCOEFF(sf), MEM_TYPE)
+	call calloc (SF_YMATRIX(sf), SF_YORDER(sf) * SF_NYCOEFF(sf), MEM_TYPE)
+	call calloc (SF_COEFF(sf), SF_NXCOEFF(sf) * SF_NYCOEFF(sf), MEM_TYPE)
+
+	# allocate temporary space
+	call smark (sp)
+	call salloc (x, SF_NCOLS(sf), MEM_TYPE)
+	call salloc (y, SF_NLINES(sf), MEM_TYPE)
+
+	# calculate all possible x basis functions and store
+	do i = 1, SF_NCOLS(sf)
+	    Memr[x+i-1] = i
+
+	switch (SF_TYPE(sf)) {
+	case SF_LEGENDRE:
+	    SF_XLEFT(sf) = NULL
+	    call sf_bleg (Memr[x], SF_NCOLS(sf), SF_XORDER(sf), SF_XMAXMIN(sf),
+		SF_XRANGE(sf), XBASIS(SF_XBASIS(sf)))
+
+	case SF_CHEBYSHEV:
+	    SF_XLEFT(sf) = NULL
+	    call sf_bcheb (Memr[x], SF_NCOLS(sf), SF_XORDER(sf), SF_XMAXMIN(sf),
+		SF_XRANGE(sf), XBASIS(SF_XBASIS(sf)))
+
+	case SF_SPLINE3:
+	    call calloc (SF_XLEFT(sf), SF_NCOLS(sf), TY_INT)
+	    call sf_bspline3 (Memr[x], SF_NCOLS(sf), SF_NXPIECES(sf),
+	        -SF_XMIN(sf), SF_XSPACING(sf), XBASIS(SF_XBASIS(sf)),
+	        XLEFT(SF_XLEFT(sf)))
+
+	case SF_SPLINE1:
+	    call calloc (SF_XLEFT(sf), SF_NCOLS(sf), TY_INT)
+	    call sf_bspline1 (Memr[x], SF_NCOLS(sf), SF_NXPIECES(sf),
+	        -SF_XMIN(sf), SF_XSPACING(sf), XBASIS(SF_XBASIS(sf)),
+	        XLEFT(SF_XLEFT(sf)))
+	}
+
+	# calculate all possible y basis functions and store
+	do i = 1, SF_NLINES(sf)
+	    Memr[y+i-1] = i
+
+	switch (SF_TYPE(sf)) {
+	case SF_LEGENDRE:
+	    SF_YLEFT(sf) = NULL
+	    call sf_bleg (Memr[y], SF_NLINES(sf), SF_YORDER(sf),
+	        SF_YMAXMIN(sf), SF_YRANGE(sf), YBASIS(SF_YBASIS(sf)))
+
+	case SF_CHEBYSHEV:
+	    SF_YLEFT(sf) = NULL
+	    call sf_bcheb (Memr[y], SF_NLINES(sf), SF_YORDER(sf),
+	        SF_YMAXMIN(sf), SF_YRANGE(sf), YBASIS(SF_YBASIS(sf)))
+
+	case SF_SPLINE3:
+	    call calloc (SF_YLEFT(sf), SF_NLINES(sf), TY_INT)
+	    call sf_bspline3 (Memr[y], SF_NLINES(sf), SF_NYPIECES(sf),
+	        -SF_YMIN(sf), SF_YSPACING(sf), YBASIS(SF_YBASIS(sf)),
+	        YLEFT(SF_YLEFT(sf)))
+
+	case SF_SPLINE1:
+	    call calloc (SF_YLEFT(sf), SF_NLINES(sf), TY_INT)
+	    call sf_bspline1 (Memr[y], SF_NLINES(sf), SF_NYPIECES(sf),
+	        -SF_YMIN(sf), SF_YSPACING(sf), YBASIS(SF_YBASIS(sf)),
+	        YLEFT(SF_YLEFT(sf)))
+	}
+
+	SF_WZ(sf) = NULL
+	SF_TLEFT(sf) = NULL
+
+	call sfree (sp)
+end
diff --git a/math/surfit/islaccum.x b/math/surfit/islaccum.x
new file mode 100644
index 00000000..cc69323a
--- /dev/null
+++ b/math/surfit/islaccum.x
@@ -0,0 +1,117 @@
+# Copyright(c) 1986 Association of Universities for Research in Astronomy Inc.
+
+include <math/surfit.h>
+include "surfitdef.h"
+
+# ISLACCUM -- Procedure to add points on a line to the data set.
+# The inner products of the non-zero x basis functions are stored
+# in the SF_XORDER(sf) by SF_NXCOEFF(sf) matrix XMATRIX. The
+# main diagonal is stored in the first row of XMATRIX. Successive
+# non-zero diagonals are stored in the succeeding rows. This method
+# of storage is particularly efficient for the large symmetric
+# banded matrices produced during spline fits. The inner products
+# of the data ordinates and the non-zero x basis functions are
+# caculated and stored in the lineno-th row of the SF_NXCOEFF(sf)
+# by SF_NLINES(sf) matrix XCOEFF.
+
+procedure islaccum (sf, cols, lineno, z, w, ncols, wtflag)
+
+pointer	sf		# pointer to surface descriptor
+int	cols[ncols]	# column values
+int	lineno		# lineno of data being accumulated
+real	z[ncols]	# surface values on lineno at cols
+real	w[ncols]	# weight of the data points
+int	ncols		# number of data points
+int	wtflag		# type of weighting desired
+
+int	i, ii, j, k
+pointer	xbzptr, xbptr
+pointer	xlzptr
+pointer	xmzptr, xmindex
+pointer	xczptr, xcindex
+pointer	bw, rows, left
+pointer sp
+
+begin
+	# count the number of points
+	SF_NXPTS(sf) = SF_NXPTS(sf) + ncols
+
+	# calculate the weights, default is uniform weighting
+	switch (wtflag) {
+	case SF_UNIFORM:
+	    call amovkr (1.0, w, ncols)
+	case SF_USER:
+	    # do not alter weights
+	default:
+	    call amovkr (1.0, w, ncols)
+	}
+
+	# set up temporary storage
+	call smark (sp)
+	call salloc (bw, ncols, TY_REAL)
+	call salloc (left, ncols, TY_INT)
+	call salloc (rows, ncols, TY_INT)
+
+	# set up the pointers
+	xbzptr = SF_XBASIS(sf) - 1
+	xmzptr = SF_XMATRIX(sf)
+	xczptr = SF_XCOEFF(sf) + (lineno - 1) * SF_NXCOEFF(sf) - 1
+
+	# accumulate the line
+	switch (SF_TYPE(sf)) {
+	case SF_LEGENDRE, SF_CHEBYSHEV:
+
+	    do i = 1, SF_XORDER(sf) {
+		do j = 1, ncols
+		    Memr[bw+j-1] = w[j] * XBASIS(xbzptr+cols[j]) 
+		xcindex = xczptr + i
+		do j = 1, ncols
+		    XCOEFF(xcindex) = XCOEFF(xcindex) + Memr[bw+j-1] * z[j]
+		xbptr  = xbzptr
+		ii = 0
+		do k = i, SF_XORDER(sf) {
+		    xmindex = xmzptr + ii
+		    do j = 1, ncols
+			XMATRIX(xmindex) = XMATRIX(xmindex) + Memr[bw+j-1] *
+			    XBASIS(xbptr+cols[j])
+		    ii = ii + 1
+		    xbptr = xbptr + SF_NCOLS(sf)
+		}
+		xbzptr = xbzptr + SF_NCOLS(sf)
+		xmzptr = xmzptr + SF_XORDER(sf)
+	    }
+
+	case SF_SPLINE3, SF_SPLINE1:
+
+	    xlzptr = SF_XLEFT(sf) - 1
+	    do j = 1, ncols
+	        Memi[left+j-1] = XLEFT(xlzptr+cols[j])
+	    call amulki (Memi[left], SF_XORDER(sf), Memi[rows], ncols)
+	    call aaddki (Memi[rows], SF_XMATRIX(sf), Memi[rows], ncols)
+	    call aaddki (Memi[left], xczptr, Memi[left], ncols) 
+
+	    do i = 1, SF_XORDER(sf) {
+		do j = 1, ncols {
+		    Memr[bw+j-1] = w[j] * XBASIS(xbzptr+cols[j])
+		    xcindex = Memi[left+j-1] + i
+		    XCOEFF(xcindex) = XCOEFF(xcindex) + Memr[bw+j-1] * z[j]
+		}
+		xbptr = xbzptr
+		ii = 0
+		do k = i, SF_XORDER(sf) {
+		    do j = 1, ncols {
+			xmindex = Memi[rows+j-1] + ii
+			XMATRIX(xmindex) = XMATRIX(xmindex) + Memr[bw+j-1] *
+			    XBASIS(xbptr+cols[j])
+		    }
+		    ii = ii + 1
+		    xbptr = xbptr + SF_NCOLS(sf)
+		}
+		xbzptr = xbzptr + SF_NCOLS(sf)
+		call aaddki (Memi[rows], SF_XORDER(sf), Memi[rows], ncols)
+	    }
+	}
+
+	# release space
+	call sfree (sp)
+end
diff --git a/math/surfit/islfit.x b/math/surfit/islfit.x
new file mode 100644
index 00000000..7b60c361
--- /dev/null
+++ b/math/surfit/islfit.x
@@ -0,0 +1,150 @@
+# Copyright(c) 1986 Association of Universities for Research in Astronomy Inc.
+
+include <math/surfit.h>
+include "surfitdef.h"
+
+# ISLFIT -- Procedure to fit a single line of a surface. The inner products
+# of the x basis functions are calculated and accumulated into the
+# SF_XORDER(sf) by SF_NXCOEFF(sf) matrix XMATRIX. The main diagonal is stored
+# in the first row of XMATRIX followed by the non-zero lower diagonals. This
+# method of storage is very efficient for the large symmetric banded matrices
+# generated by spline fits. The Cholesky factorization of XMATRIX is calculated
+# and stored in XMATRIX destroying the original data. The inner products
+# of the x basis functions and the data ordinates are stored in the lineno-th
+# row of the SF_NXCOEFF(sf) by SF_NLINES(sf) matrix XCOEFF. The x coefficients
+# for each line are calculated by forward and back substitution and
+# stored in the lineno-th line of XCOEFF destroying the original data.
+
+procedure islfit (sf, cols, lineno, z, w, ncols, wtflag, ier)
+
+pointer	sf		# pointer to the surface descriptor
+int	cols[ncols]	# array of columns
+int	lineno    	# lineno
+real	z[ncols]	# the surface values
+real	w[ncols]	# array of weights
+int	ncols		# the number of columns
+int	wtflag		# type of weighting
+int	ier		# error codes
+
+int	i, ii, j, k
+pointer	xbzptr, xbptr
+pointer	xmzptr, xmindex
+pointer	xczptr, xcindex
+pointer	xlzptr
+pointer	left, rows, bw
+pointer	sp
+real	adotr()
+
+begin
+	# index the pointers
+	xbzptr = SF_XBASIS(sf) - 1
+	xmzptr = SF_XMATRIX(sf)
+	xczptr = SF_XCOEFF(sf) + (lineno - 1) * SF_NXCOEFF(sf) - 1
+
+	# zero the accumulators
+	call aclrr (XMATRIX(SF_XMATRIX(sf)), SF_NXCOEFF(sf) * SF_XORDER(sf))
+	call aclrr (XCOEFF(xczptr + 1), SF_NXCOEFF(sf))
+
+	# free space used by previous islrefit calls
+	if (SF_WZ(sf) != NULL)
+	    call mfree (SF_WZ(sf), MEM_TYPE)
+	if (SF_TLEFT(sf) != NULL)
+	    call mfree (SF_TLEFT(sf), TY_INT)
+
+	# reset number of points
+	SF_NXPTS(sf) = ncols
+
+	# calculate the weights, default is uniform weighting
+	switch (wtflag) {
+	case SF_UNIFORM:
+	    call amovkr (1.0, w, ncols)
+	case SF_USER:
+	    # do not alter weights
+	default:
+	    call amovkr (1.0, w, ncols)
+	}
+
+	# allocate temporary space
+	call smark (sp)
+	call salloc (bw, ncols, TY_REAL)
+
+	switch (SF_TYPE(sf)) {
+	case SF_LEGENDRE, SF_CHEBYSHEV:
+
+	    do i = 1, SF_XORDER(sf) {
+		do j = 1, ncols
+		    Memr[bw+j-1] = w[j] * XBASIS(xbzptr+cols[j])
+
+		xcindex = xczptr + i
+		XCOEFF(xcindex) = XCOEFF(xcindex) + adotr (Memr[bw], z, ncols)
+
+		xbptr = xbzptr
+		ii = 0
+		do k = i, SF_XORDER(sf) {
+		    xmindex = xmzptr + ii
+		    do j = 1, ncols
+		        XMATRIX(xmindex) = XMATRIX(xmindex) + Memr[bw+j-1] *
+			    XBASIS(xbptr+cols[j])
+		    ii = ii + 1
+		    xbptr = xbptr + SF_NCOLS(sf)
+		}
+
+		xbzptr = xbzptr + SF_NCOLS(sf)
+		xmzptr = xmzptr + SF_XORDER(sf)
+	    }
+
+	case SF_SPLINE3, SF_SPLINE1:
+	    xlzptr = SF_XLEFT(sf) - 1
+
+	    call salloc (left, ncols, TY_INT)
+	    call salloc (rows, ncols, TY_INT)
+
+	    do j = 1, ncols
+		Memi[left+j-1] = XLEFT(xlzptr+cols[j])
+	    call amulki (Memi[left], SF_XORDER(sf), Memi[rows], ncols)
+	    call aaddki (Memi[rows], SF_XMATRIX(sf), Memi[rows], ncols)
+	    call aaddki (Memi[left], xczptr, Memi[left], ncols)
+
+	    do i = 1, SF_XORDER(sf) {
+		do j = 1, ncols {
+		    Memr[bw+j-1] = w[j] * XBASIS(xbzptr+cols[j])
+		    xcindex = Memi[left+j-1] + i
+		    XCOEFF(xcindex) = XCOEFF(xcindex) + Memr[bw+j-1] * z[j]
+		}
+
+		xbptr = xbzptr
+		ii = 0
+		do k = i, SF_XORDER(sf) {
+		    do j = 1, ncols {
+			xmindex = Memi[rows+j-1] + ii
+			XMATRIX(xmindex) = XMATRIX(xmindex) + Memr[bw+j-1] *
+			    XBASIS(xbptr+cols[j])
+		    }
+		    ii = ii + 1
+		    xbptr = xbptr + SF_NCOLS(sf)
+		}
+
+		xbzptr = xbzptr + SF_NCOLS(sf)
+		call aaddki (Memi[rows], SF_XORDER(sf), Memi[rows], ncols)
+	    }
+	}
+
+	# release space
+	call sfree (sp)
+
+	# return if not enough data points
+	ier = OK
+	if ((SF_NXPTS(sf) - SF_NXCOEFF(sf)) < 0) {
+	    ier = NO_DEG_FREEDOM
+	    return
+	}
+
+	# calculate the Cholesky factorization of XMATRIX
+	call sfchofac (XMATRIX(SF_XMATRIX(sf)), SF_XORDER(sf), SF_NXCOEFF(sf),
+	    XMATRIX(SF_XMATRIX(sf)), ier)
+
+	# calculate the x coefficients for lineno-th image line and place in the
+	# lineno-th row of XCOEFF
+	call sfchoslv (XMATRIX(SF_XMATRIX(sf)), SF_XORDER(sf), SF_NXCOEFF(sf),
+	    XCOEFF(xczptr + 1), XCOEFF(xczptr + 1))
+end
diff --git a/math/surfit/islrefit.x b/math/surfit/islrefit.x
new file mode 100644
index 00000000..b0430e5b
--- /dev/null
+++ b/math/surfit/islrefit.x
@@ -0,0 +1,74 @@
+# Copyright(c) 1986 Association of Universities for Research in Astronomy Inc.
+
+include <math/surfit.h>
+include "surfitdef.h"
+
+# ISLREFIT -- Procedure to refit the data assuming that the cols and w
+# arrays do not change. SIFLREFIT assumes that the Cholesky factorization
+# of XMATRIX is stored in XMATRIX. The inner products of the x basis
+# functions and the data ordinates are accumulated into the lineno-th
+# row of the SF_NXCOEFF(sf) by SF_NLINES(sf) matrix XCOEFF. The coefficients
+# for line number lineno are calculated by forward and back
+# substitution and placed in the lineno-th row of XCOEFF replacing the
+# original data.
+
+procedure islrefit (sf, cols, lineno, z, w)
+
+pointer	sf		# pointer to surface descriptor
+int	cols[ARB]	# columns to be fit 
+int	lineno		# line number
+real	z[ARB]		# surface values
+real	w[ARB]		# weight values
+
+int	i, j
+pointer	xbzptr, xczptr, xcindex, xlzptr
+
+begin
+	# set pointers
+	xbzptr = SF_XBASIS(sf) - 1
+	xczptr = SF_XCOEFF(sf) + (lineno - 1) * SF_NXCOEFF(sf) - 1
+	xlzptr = SF_XLEFT(sf) - 1
+
+	# reset lineno-th row of the x coefficient matrix
+	call aclrr (XCOEFF(xczptr+1), SF_NXCOEFF(sf))
+
+	if (SF_WZ(sf) == NULL)
+	    call malloc (SF_WZ(sf), SF_NXPTS(sf), MEM_TYPE)
+
+	# calculate new right sides
+	call amulr (w, z, Memr[SF_WZ(sf)], SF_NXPTS(sf))
+
+	switch (SF_TYPE(sf)) {
+	case SF_LEGENDRE, SF_CHEBYSHEV:
+
+	    do i = 1, SF_XORDER(sf) {
+		xcindex = xczptr + i
+		do j = 1, SF_NXPTS(sf)
+		    XCOEFF(xcindex) = XCOEFF(xcindex) + Memr[SF_WZ(sf)+j-1] *
+			XBASIS(xbzptr+cols[j])
+		xbzptr = xbzptr + SF_NCOLS(sf)
+	    }
+
+	case SF_SPLINE3, SF_SPLINE1:
+
+	    if (SF_TLEFT(sf) == NULL)
+	        call malloc (SF_TLEFT(sf), SF_NXPTS(sf), TY_INT)
+
+	    do i = 1, SF_NXPTS(sf)
+		Memi[SF_TLEFT(sf)+i-1] = XLEFT(xlzptr+cols[i]) + xczptr
+
+	    do i = 1, SF_XORDER(sf) {
+		do j = 1, SF_NXPTS(sf) {
+		    xcindex = Memi[SF_TLEFT(sf)+j-1] + i
+		    XCOEFF(xcindex) = XCOEFF(xcindex) + Memr[SF_WZ(sf)+j-1] *
+			XBASIS(xbzptr+cols[j])
+		}
+		xbzptr = xbzptr + SF_NCOLS(sf)
+	    }
+
+	}
+
+	# solve for the new x coefficients for line number lineno
+	call sfchoslv (XMATRIX(SF_XMATRIX(sf)), SF_XORDER(sf), SF_NXCOEFF(sf),
+	    XCOEFF(xczptr+1), XCOEFF(xczptr+1))
+end
diff --git a/math/surfit/islsolve.x b/math/surfit/islsolve.x
new file mode 100644
index 00000000..11d0d313
--- /dev/null
+++ b/math/surfit/islsolve.x
@@ -0,0 +1,48 @@
+# Copyright(c) 1986 Association of Universities for Research in Astronomy Inc.
+
+include <math/surfit.h>
+include "surfitdef.h"
+
+# ISLSOLVE -- Procedure to  solve for the x coefficients of image line
+# number lineno. The inner products of the x basis functions are assumed
+# to be in the SF_XORDER(sf) by SF_NXCOEFF(sf) array XMATRIX,
+# while the inner products of the basis functions and
+# the data ordinated for line number lineno are assumed to be in the
+# lineno-th row of the SF_NXCOEFF(sf) by SF_NLINES(sf) matrix XCOEFF.
+# The Cholesky factorization of XMATRIX is calculated and placed
+# in XMATRIX overwriting the original data. The x coefficients for
+# line number lineno are calculated and placed in the lineno-th row
+# of XCOEFF replacing the original data.
+
+procedure islsolve (sf, lineno, ier)
+
+pointer	sf 		# pointer to the surface descriptor structure
+int	lineno		# line being fitted in x
+int	ier		# ier = 0, everything OK
+			# ier = 1, matrix is singular
+			# ier = 2, no degree of freedom
+
+pointer	xcptr
+
+begin
+	# return if there are insuffucient points to solve the matrix
+	ier = OK
+	if ((SF_NXPTS(sf) - SF_NXCOEFF(sf)) < 0 ) {
+	    ier = NO_DEG_FREEDOM
+	    return
+	}
+
+	# calculate the Cholesky factorization of the x matrix and store
+	# separately for possible use by SFLREFIT
+
+	call sfchofac (XMATRIX(SF_XMATRIX(sf)), SF_XORDER(sf), SF_NXCOEFF(sf),
+	    XMATRIX(SF_XMATRIX(sf)), ier)
+
+	# solve for the x coefficients for line lineno assuming the
+	# data are in row lineno of xcoeff, the solution is placed
+	# on top of the data
+
+	xcptr = SF_XCOEFF(sf) + (lineno - 1) * SF_NXCOEFF(sf)
+	call sfchoslv (XMATRIX(SF_XMATRIX(sf)), SF_XORDER(sf), SF_NXCOEFF(sf),
+	    XCOEFF(xcptr), XCOEFF(xcptr))
+end
diff --git a/math/surfit/islzero.x b/math/surfit/islzero.x
new file mode 100644
index 00000000..03034f01
--- /dev/null
+++ b/math/surfit/islzero.x
@@ -0,0 +1,25 @@
+# Copyright(c) 1986 Association of Universities for Research in Astronomy Inc.
+
+include "surfitdef.h"
+
+# ISLZERO -- Procedure to zero the accumulators for line number lineno.
+
+procedure islzero (sf, lineno)
+
+pointer	sf	# pointer to the surface descriptor
+int	lineno	# line number
+pointer	xcptr
+
+begin
+	SF_NXPTS(sf) = 0
+
+	call aclrr (XMATRIX(SF_XMATRIX(sf)),
+		     SF_XORDER(sf) * SF_NXCOEFF(sf))
+	xcptr = SF_XCOEFF(sf) + (lineno - 1) * SF_NXCOEFF(sf)
+	call aclrr (XCOEFF(xcptr), SF_NXCOEFF(sf))
+
+	if (SF_WZ(sf) != NULL)
+	    call mfree (SF_WZ(sf), MEM_TYPE)
+	if (SF_TLEFT(sf) != NULL)
+	    call mfree (SF_TLEFT(sf), TY_INT)
+end
diff --git a/math/surfit/isreplace.x b/math/surfit/isreplace.x
new file mode 100644
index 00000000..1ba69479
--- /dev/null
+++ b/math/surfit/isreplace.x
@@ -0,0 +1,114 @@
+# Copyright(c) 1986 Association of Universities for Research in Astronomy Inc.
+
+include <math/surfit.h>
+include "surfitdef.h"
+
+# ISREPLACE -- Procedure to restore the surface fit stored by SIFSAVE
+# to the surface descriptor for use by the evaluating routines. The
+# surface parameters, surface type, xorder (or number of polynomial
+# pieces in x), yorder (or number of polynomial pieces in y), xterms,
+# number of columns and number of lines, are stored in the first
+# six elements of the real array fit, followed by the SF_NYCOEFF(sf) *
+# SF_NYCOEFF(sf) surface coefficients. The coefficient of B(i,x) * B(j,y)
+# is stored in element number 6 + (i - 1) * SF_NYCOEFF(sf) + j of the
+# array fit.
+
+procedure isreplace (sf, fit)
+
+pointer	sf		# surface descriptor
+real	fit[ARB]	# array containing the surface parameters and
+			# coefficients
+
+int	surface_type, xorder, yorder, ncols, nlines
+
+begin
+	# allocate space for the surface descriptor
+	call calloc (sf, LEN_SFSTRUCT, TY_STRUCT)
+
+	xorder = nint (SF_SAVEXORDER(fit))
+	if (xorder < 1)
+	    call error (0, "SFRESTORE: Illegal x order.")
+	yorder = nint (SF_SAVEYORDER(fit))
+	if (yorder < 1)
+	    call error (0, "SFRESTORE: Illegal y order.")
+
+	ncols = nint (SF_SAVENCOLS(fit))
+	if (ncols < 1)
+	    call error (0, "SFRESTORE: Illegal x range.")
+	nlines = nint (SF_SAVENLINES(fit))
+	if (nlines < 1)
+	    call error (0, "SFRESTORE: Illegal y range.")
+
+	# set surface type dependent surface descriptor parameters
+	surface_type = nint (SF_SAVETYPE(fit))
+
+	switch (surface_type) {
+	case SF_LEGENDRE, SF_CHEBYSHEV:
+	    SF_NXCOEFF(sf) = xorder
+	    SF_XORDER(sf) = xorder
+	    SF_XRANGE(sf) = 2. / real (ncols + 1)
+	    SF_XMAXMIN(sf) =  - real (ncols + 1) / 2.
+	    SF_XMIN(sf) = 0.
+	    SF_XMAX(sf) = real (ncols + 1)
+	    SF_NYCOEFF(sf) = yorder
+	    SF_YORDER(sf) = yorder
+	    SF_YRANGE(sf) = 2. / real (nlines + 1)
+	    SF_YMAXMIN(sf) =  - real (nlines + 1) / 2.
+	    SF_YMIN(sf) = 0.
+	    SF_YMAX(sf) = real (nlines + 1)
+	    SF_XTERMS(sf) = SF_SAVEXTERMS(fit)
+
+	case SF_SPLINE3:
+	    SF_NXCOEFF(sf) = (xorder + SPLINE3_ORDER - 1)
+	    SF_XORDER(sf) = SPLINE3_ORDER
+	    SF_NXPIECES(sf) = xorder - 1
+	    SF_XSPACING(sf) = xorder / real (ncols + 1)
+	    SF_XMIN(sf) = 0.
+	    SF_XMAX(sf) = real (ncols + 1)
+	    SF_NYCOEFF(sf) = (yorder + SPLINE3_ORDER - 1)
+	    SF_YORDER(sf) = SPLINE3_ORDER
+	    SF_NYPIECES(sf) = yorder - 1
+	    SF_YSPACING(sf) = yorder / real (nlines + 1)
+	    SF_YMIN(sf) = 0.
+	    SF_YMAX(sf) = real (nlines + 1)
+	    SF_XTERMS(sf) = YES
+
+	case SF_SPLINE1:
+	    SF_NXCOEFF(sf) = (xorder + SPLINE1_ORDER - 1)
+	    SF_XORDER(sf) = SPLINE1_ORDER
+	    SF_NXPIECES(sf) = xorder - 1
+	    SF_XSPACING(sf) = xorder / real (ncols + 1)
+	    SF_XMIN(sf) = 0.
+	    SF_XMAX(sf) = real (ncols + 1)
+	    SF_NYCOEFF(sf) = (yorder + SPLINE1_ORDER - 1)
+	    SF_YORDER(sf) = SPLINE1_ORDER
+	    SF_NYPIECES(sf) = yorder - 1
+	    SF_YSPACING(sf) = yorder / real (nlines + 1)
+	    SF_YMIN(sf) = 0.
+	    SF_YMAX(sf) = real (nlines + 1)
+	    SF_XTERMS(sf) = YES
+
+	default:
+	    call error (0, "SFRESTORE: Unknown surface type.")
+	}
+
+	# set remaining curve parameters
+	SF_TYPE(sf) = surface_type
+	SF_NLINES(sf) = nlines
+	SF_NCOLS(sf) = ncols
+
+	# allocate space for the coefficient array
+	SF_XBASIS(sf) = NULL
+	SF_YBASIS(sf) = NULL
+	SF_XMATRIX(sf) = NULL
+	SF_YMATRIX(sf) = NULL
+	SF_XCOEFF(sf) = NULL
+	SF_WZ(sf) = NULL
+	SF_TLEFT(sf) = NULL
+
+	call calloc (SF_COEFF(sf), SF_NXCOEFF(sf) * SF_NYCOEFF(sf), MEM_TYPE)
+
+	# restore coefficient array
+	call amovr (fit[SF_SAVECOEFF+1], COEFF(SF_COEFF(sf)), SF_NYCOEFF(sf) *
+	    SF_NXCOEFF(sf))
+end
diff --git a/math/surfit/isresolve.x b/math/surfit/isresolve.x
new file mode 100644
index 00000000..7af6a8da
--- /dev/null
+++ b/math/surfit/isresolve.x
@@ -0,0 +1,127 @@
+# Copyright(c) 1986 Association of Universities for Research in Astronomy Inc.
+
+include <math/surfit.h>
+include "surfitdef.h"
+
+# ISRESOLVE -- Procedure to solve the x coefficient matrix for the surface
+# coefficients assuming that the lines array is unchanged since the last
+# call to SIFSOLVE. The Cholesky factorization of YMATRIX is assumed to be
+# stored in YMATRIX. The inner product of the y basis functions and i-th
+# column of XCOEFF containing the i-th x coefficients for each line are
+# calculated and stored in the i-th row of the SF_NYCOEFF(sf) by
+# SF_NXCOEFF(sf) array COEFF. Each of the SF_NXCOEFF(sf) rows of COEFF
+# is solved to determine the SF_NYCOEFF(sf) by SF_NXCOEFF(sf) surface
+# coefficients. After a call to SIFSOLVE the coefficient for the i-th
+# y basis function and the j-th x coefficient is found in the j-th row
+# and i-th column of COEFF.
+
+procedure isresolve (sf, lines, ier)
+
+pointer	sf		# pointer to the surface descriptor structure
+int	lines[ARB]	# line numbers included in the fit
+int	ier		# error code
+
+int	i, j, k, nxcoeff
+pointer	ybzptr
+pointer	ylzptr
+pointer	xczptr, xcptr, xcindex
+pointer	czptr, cptr
+pointer	left, tleft
+pointer	sp
+
+begin
+	# define pointers
+	ybzptr = SF_YBASIS(sf) - 1
+	xczptr = SF_XCOEFF(sf) - SF_NXCOEFF(sf) - 1
+	czptr = SF_COEFF(sf) - 1
+
+	# zero out coefficient matrix
+	call aclrr (COEFF(SF_COEFF(sf)), SF_NXCOEFF(sf) * SF_NYCOEFF(sf))
+
+	switch (SF_TYPE(sf)) {
+	case SF_LEGENDRE, SF_CHEBYSHEV:
+
+	    # loop over the y basis functions
+	    nxcoeff = SF_NXCOEFF(sf)
+	    do i = 1, SF_YORDER(sf) {
+		cptr = czptr + i
+		do k = 1, nxcoeff {
+		    xcptr = xczptr + k
+		    do j = 1, SF_NYPTS(sf) {
+			xcindex = xcptr + lines[j] * SF_NXCOEFF(sf)
+		        COEFF(cptr) = COEFF(cptr) +
+		    		  YBASIS(ybzptr+lines[j]) * XCOEFF(xcindex)
+		    }
+		    cptr = cptr + SF_NYCOEFF(sf)
+		}
+
+		ybzptr = ybzptr + SF_NYPTS(sf)
+
+		# if SF_XTERMS(sf) = NO do not accumulate elements
+		# of COEFF where i != 1 and k != 1
+		if (SF_XTERMS(sf) == NO)
+		    nxcoeff = 1
+	    }
+
+	case SF_SPLINE3, SF_SPLINE1:
+	    call smark (sp)
+	    call salloc (left, SF_NYPTS(sf), TY_INT)
+	    call salloc (tleft, SF_NYPTS(sf), TY_INT)
+
+	    ylzptr = SF_YLEFT(sf) - 1
+	    do j = 1, SF_NYPTS(sf)
+		Memi[left+j-1] = YLEFT(ylzptr+lines[j])
+	    call aaddki (Memi[left], czptr, Memi[left], SF_NYPTS(sf))
+
+	    nxcoeff = SF_NXCOEFF(sf)
+	    do i = 1, SF_YORDER(sf) {
+		call aaddki (Memi[left], i, Memi[tleft], SF_NYPTS(sf))
+		do k = 1, nxcoeff {
+		    xcptr = xczptr + k
+		    do j = 1, SF_NYPTS(sf) {
+			cptr = Memi[tleft+j-1]
+			xcindex = xcptr + lines[j] * SF_NXCOEFF(sf)
+			COEFF(cptr) = COEFF(cptr) + YBASIS(ybzptr+lines[j]) *
+			    XCOEFF(xcindex)
+		    }
+		    call aaddki (Memi[tleft], SF_NYCOEFF(sf), Memi[tleft],
+			SF_NYPTS(sf))
+		}
+
+		ybzptr = ybzptr + SF_NYPTS(sf)
+	    }
+
+	    call sfree (sp)
+	}
+
+	# return if not enough points
+	ier = OK
+	if ((SF_NYPTS(sf) - SF_NYCOEFF(sf)) < 0) {
+	    ier = NO_DEG_FREEDOM
+	    return
+	}
+
+	if (SF_XTERMS(sf) == YES) {
+
+	    # solve for the nxcoeff right sides
+	    cptr = SF_COEFF(sf)
+	    do i = 1, SF_NXCOEFF(sf) {
+	        call sfchoslv (YMATRIX(SF_YMATRIX(sf)), SF_YORDER(sf),
+	    		       SF_NYCOEFF(sf), COEFF(cptr), COEFF(cptr))
+	        cptr = cptr + SF_NYCOEFF(sf)
+	    }
+
+	} else {
+
+	    # solve for the y coefficients
+	    cptr = SF_NYCOEFF(sf)
+	    call sfchoslv (YMATRIX(SF_YMATRIX(sf)), SF_YORDER(sf),
+			      SF_NYCOEFF(sf), COEFF(cptr), COEFF(cptr))
+
+	    # solve for the x coefficients
+	    do i = 2, SF_NXCOEFF(sf) {
+		cptr = czptr + SF_NYCOEFF(sf)
+		COEFF(cptr) = COEFF(cptr) / SF_NYPTS(sf)
+	    }
+	}
+end
diff --git a/math/surfit/issave.x b/math/surfit/issave.x
new file mode 100644
index 00000000..5a3ecab8
--- /dev/null
+++ b/math/surfit/issave.x
@@ -0,0 +1,44 @@
+# Copyright(c) 1986 Association of Universities for Research in Astronomy Inc.
+
+include	<math/surfit.h>
+include	"surfitdef.h"
+
+# ISSAVE -- Procedure to save the surface fit for later use by the
+# evaluate routines. After a call to SIFSAVE the first six elements
+# of fit contain the surface type, xorder (or number of polynomial pieces
+# in x), yorder (or the number of polynomial pieces in y), xterms, ncols
+# and nlines. The remaining spaces are filled by the SF_NYCOEFF(sf) *
+# SF_NXCOEFF(sf) surface coefficients. The coefficient of B(i,x) * B(j,y)
+# is located in element number 6 + (i - 1) * SF_NYCOEFF(sf) + j of the
+# array fit where i <= SF_NXCOEFF(sf) and j <= SF_NYCOEFF(sf).
+
+procedure issave (sf, fit)
+
+pointer	sf		# pointer to the surface descriptor
+real	fit[ARB]	# array for storing fit
+
+begin
+	# get the surface parameters
+
+	# order is surface type dependent
+	switch (SF_TYPE(sf)) {
+	case SF_LEGENDRE, SF_CHEBYSHEV:
+	    SF_SAVEXORDER(fit) = SF_XORDER(sf)
+	    SF_SAVEYORDER(fit) = SF_YORDER(sf)
+	case SF_SPLINE3, SF_SPLINE1:
+	    SF_SAVEXORDER(fit) = SF_NXPIECES(sf) + 1
+	    SF_SAVEYORDER(fit) = SF_NYPIECES(sf) + 1
+	default:
+	    call error (0, "SIFSAVE: Unknown surface type.")
+	}
+
+	# save remaining parameters
+	SF_SAVETYPE(fit) = SF_TYPE(sf)
+	SF_SAVENLINES(fit) = SF_NLINES(sf)
+	SF_SAVENCOLS(fit) = SF_NCOLS(sf)
+	SF_SAVEXTERMS(fit) = SF_XTERMS(sf)
+
+	# save the coefficients
+	call amovr (COEFF(SF_COEFF(sf)), fit[SF_SAVECOEFF+1], SF_NXCOEFF(sf) *
+	    SF_NYCOEFF(sf))
+end
diff --git a/math/surfit/issolve.x b/math/surfit/issolve.x
new file mode 100644
index 00000000..7790ef60
--- /dev/null
+++ b/math/surfit/issolve.x
@@ -0,0 +1,169 @@
+# Copyright(c) 1986 Association of Universities for Research in Astronomy Inc.
+
+include <math/surfit.h>
+include "surfitdef.h"
+
+# ISSOLVE -- Procedure to solve the x coefficient matrix for the surface
+# coefficients. The inner products of the y basis functions are accumulated
+# and stored in the SF_YORDER(sf) by SF_NYCOEFF(sf) array YMATRIX. The
+# main diagonal of YMATRIX is stored in the first row of YMATRIX followed
+# by the remaining non-zero lower diagonals. The Cholesky factorization
+# of YMATRIX is calculated and stored on top of YMATRIX destroying the
+# original data. The inner products of the y basis functions and 
+# and the i-th column of the SF_NXCOEFF(sf) by SF_NLINES(sf) matrix XCOEFF
+# containing the i-th x coefficients for each line are calculated and
+# placed in the i-th row of the SF_NYCOEFF(sf) by SF_NXCOEFF(sf) matrix
+# COEFF. Each of the SF_NXCOEFF(sf) rows of COEFF is solved to determine
+# the SF_NXCOEFF(sf) by SF_NYCOEFF(sf) surface coefficients. After a
+# call to SIFSOLVE the coefficient of the i-th x basis function and the
+# j-th y basis function will be found in the j-th column and i-th row
+# of COEFF.
+
+procedure issolve (sf, lines, nlines, ier)
+
+pointer	sf		# pointer to the curve descriptor structure
+int	lines[ARB]	# line numbers included in the fit
+int	nlines		# number of lines fit
+int	ier		# error code
+
+int	i, ii, j, k, nxcoeff
+pointer	ybzptr, ybptr
+pointer	ylzptr
+pointer	ymzptr, ymindex
+pointer	xczptr, xcptr, xcindex
+pointer	czptr, cptr
+pointer	left, tleft, rows
+pointer	sp
+
+begin
+	# define pointers
+	ybzptr = SF_YBASIS(sf) - 1
+	ymzptr = SF_YMATRIX(sf)
+	xczptr = SF_XCOEFF(sf) - SF_NXCOEFF(sf) - 1
+	czptr = SF_COEFF(sf) - 1
+
+	# zero out coefficient matrix and the y coefficient matrix
+	call aclrr (YMATRIX(SF_YMATRIX(sf)), SF_YORDER(sf) * SF_NYCOEFF(sf))
+	call aclrr (COEFF(SF_COEFF(sf)), SF_NXCOEFF(sf) * SF_NYCOEFF(sf))
+
+	# increment the number of points
+	SF_NYPTS(sf) =  nlines
+
+	switch (SF_TYPE(sf)) {
+	case SF_LEGENDRE, SF_CHEBYSHEV:
+
+	    # accumulate the y value in the y matrix 
+	    nxcoeff = SF_NXCOEFF(sf)
+	    do i = 1, SF_YORDER(sf) {
+		cptr = czptr + i
+		do k = 1, nxcoeff {
+		    xcptr = xczptr + k
+		    do j = 1, nlines {
+		        xcindex = xcptr + lines[j] * SF_NXCOEFF(sf)
+		        COEFF(cptr) = COEFF(cptr) +
+			    YBASIS(ybzptr+lines[j]) * XCOEFF(xcindex)
+		    }
+		    cptr = cptr + SF_NYCOEFF(sf)
+		}
+		ii = 0
+		ybptr = ybzptr
+		do k = i, SF_YORDER(sf) {
+		    ymindex = ymzptr + ii
+		    do j = 1, nlines  
+		        YMATRIX(ymindex) = YMATRIX(ymindex) +
+				     YBASIS(ybzptr+lines[j]) *
+				     YBASIS(ybptr+lines[j])
+		    ii = ii + 1
+		    ybptr = ybptr + SF_NLINES(sf)
+		}
+
+		if (SF_XTERMS(sf) == NO)
+		    nxcoeff = 1
+
+		ybzptr = ybzptr + SF_NLINES(sf)
+		ymzptr = ymzptr + SF_YORDER(sf)
+	    }
+
+	case SF_SPLINE3, SF_SPLINE1:
+
+	    call smark (sp)
+	    call salloc (left, nlines, TY_INT)
+	    call salloc (tleft, nlines, TY_INT)
+	    call salloc (rows, nlines, TY_INT)
+
+	    ylzptr = SF_YLEFT(sf) - 1
+	    do j = 1, nlines
+		Memi[left+j-1] = YLEFT(ylzptr+lines[j])
+	    call amulki (Memi[left], SF_YORDER(sf), Memi[rows], nlines)
+	    call aaddki (Memi[rows], SF_YMATRIX(sf), Memi[rows], nlines)
+	    call aaddki (Memi[left], czptr, Memi[left], nlines)
+
+	    # accumulate the y value in the y matrix 
+	    nxcoeff = SF_NXCOEFF(sf)
+	    do i = 1, SF_YORDER(sf) {
+		call aaddki (Memi[left], i, Memi[tleft], nlines)
+		do k = 1, nxcoeff {
+		    xcptr = xczptr + k
+		    do j = 1, nlines {
+			cptr = Memi[tleft+j-1]
+			xcindex = xcptr + lines[j] * SF_NXCOEFF(sf)
+		        COEFF(cptr) = COEFF(cptr) + YBASIS(ybzptr+lines[j]) *
+			    XCOEFF(xcindex)
+		    }
+		    call aaddki (Memi[tleft], SF_NYCOEFF(sf), Memi[tleft],
+		        nlines)
+		}
+		ii = 0
+		ybptr = ybzptr
+		do k = i, SF_YORDER(sf) {
+		    do j = 1, nlines {  
+		        ymindex = Memi[rows+j-1] + ii
+		        YMATRIX(ymindex) = YMATRIX(ymindex) +
+				     YBASIS(ybzptr+lines[j]) *
+				     YBASIS(ybptr+lines[j])
+		    }
+		    ii = ii + 1
+		    ybptr = ybptr + SF_NLINES(sf)
+		}
+
+		ybzptr = ybzptr + SF_NLINES(sf)
+		call aaddki (Memi[rows], SF_YORDER(sf), Memi[rows], nlines)
+	    }
+
+	    call sfree (sp)
+
+	}
+
+	# return if not enough points
+	ier = OK
+	if ((SF_NYPTS(sf) - SF_NYCOEFF(sf)) < 0) {
+	    ier = NO_DEG_FREEDOM
+	    return
+	}
+
+	# calculate the Cholesky factorization of the y matrix
+	call sfchofac (YMATRIX(SF_YMATRIX(sf)), SF_YORDER(sf), SF_NYCOEFF(sf),
+		       YMATRIX(SF_YMATRIX(sf)), ier)
+
+	if (SF_XTERMS(sf) == YES) {
+
+	    # solve for the nxcoeff right sides
+	    cptr = SF_COEFF(sf)
+	    do i = 1, SF_NXCOEFF(sf) {
+	        call sfchoslv (YMATRIX(SF_YMATRIX(sf)), SF_YORDER(sf),
+	    		       SF_NYCOEFF(sf), COEFF(cptr), COEFF(cptr))
+		cptr = cptr + SF_NYCOEFF(sf)
+	    }
+
+	} else {
+
+	    cptr = SF_COEFF(sf)
+	    call sfchoslv (YMATRIX(SF_YMATRIX(sf)), SF_YORDER(sf),
+			      SF_NYCOEFF(sf), COEFF(cptr), COEFF(cptr))
+
+	    do i = 2, SF_NXCOEFF(sf) {
+		cptr = cptr + SF_NYCOEFF(sf)
+		COEFF(cptr) = COEFF(cptr) / SF_NYPTS(sf)
+	    }
+	}
+end
diff --git a/math/surfit/isvector.x b/math/surfit/isvector.x
new file mode 100644
index 00000000..023d3f4d
--- /dev/null
+++ b/math/surfit/isvector.x
@@ -0,0 +1,76 @@
+# Copyright(c) 1986 Association of Universities for Research in Astronomy Inc.
+
+include <math/surfit.h>
+include "surfitdef.h"
+
+# ISVECTOR -- Procedure to evaluate the fitted surface at an array of points.
+# The SF_NXCOEFF(sf) by SF_NYCOEFF(sf) coefficients are stored in the
+# SF_NYCOEFF(sf) by SF_NXCOEFF(sf) matrix COEFF. The j-th element of the ith
+# row of COEFF contains the coefficient of the i-th basis function in x and
+# the j-th basis function in y.
+
+procedure isvector (sf, x, y, zfit, npts)
+
+pointer	sf		# pointer to surface descriptor structure
+real	x[ARB]		# x value
+real	y[ARB]		# y value
+real	zfit[ARB]	# fits surface values
+int	npts		# number of data points
+
+int	i
+pointer	xcoeff, cptr, sp
+
+begin
+	# evaluate the surface along the vector
+	switch (SF_TYPE(sf)) {
+	case SF_CHEBYSHEV:
+	    if (SF_XORDER(sf) == 1) {
+		call cv_evcheb (COEFF(SF_COEFF(sf)), y, zfit, npts,
+		    SF_YORDER(sf), SF_YMAXMIN(sf), SF_YRANGE(sf))
+	    } else if (SF_YORDER(sf) == 1) {
+		call smark (sp)
+		call salloc (xcoeff, SF_NXCOEFF(sf), MEM_TYPE)
+		cptr = SF_COEFF(sf)
+		do i = 1, SF_NXCOEFF(sf) {
+		    Memr[xcoeff+i-1] = COEFF(cptr)
+		    cptr = cptr + SF_NYCOEFF(sf)
+		}
+		call cv_evcheb (Memr[xcoeff], x, zfit, npts,
+		    SF_XORDER(sf), SF_XMAXMIN(sf), SF_XRANGE(sf))
+		call sfree (sp)
+	    } else
+	        call sf_evcheb (COEFF(SF_COEFF(sf)), x, y, zfit, npts,
+	            SF_XTERMS(sf), SF_XORDER(sf), SF_YORDER(sf), SF_XMAXMIN(sf),
+		    SF_XRANGE(sf), SF_YMAXMIN(sf), SF_YRANGE(sf))
+
+	case SF_LEGENDRE:
+	    if (SF_XORDER(sf) == 1) {
+		call cv_evleg (COEFF(SF_COEFF(sf)), y, zfit, npts,
+		    SF_YORDER(sf), SF_YMAXMIN(sf), SF_YRANGE(sf))
+	    } else if (SF_YORDER(sf) == 1) {
+		call smark (sp)
+		call salloc (xcoeff, SF_NXCOEFF(sf), MEM_TYPE)
+		cptr = SF_COEFF(sf)
+		do i = 1, SF_NXCOEFF(sf) {
+		    Memr[xcoeff+i-1] = COEFF(cptr)
+		    cptr = cptr + SF_NYCOEFF(sf)
+		}
+		call cv_evcheb (Memr[xcoeff], x, zfit, npts,
+		    SF_XORDER(sf), SF_XMAXMIN(sf), SF_XRANGE(sf))
+		call sfree (sp)
+	    } else
+	        call sf_evleg (COEFF(SF_COEFF(sf)), x, y, zfit, npts,
+		    SF_XTERMS(sf), SF_XORDER(sf), SF_YORDER(sf), SF_XMAXMIN(sf),
+		    SF_XRANGE(sf), SF_YMAXMIN(sf), SF_YRANGE(sf))
+
+	case SF_SPLINE3:
+	    call sf_evspline3 (COEFF(SF_COEFF(sf)), x, y, zfit, npts,
+	        SF_NXPIECES(sf), SF_NYPIECES(sf), -SF_XMIN(sf), SF_XSPACING(sf),
+		-SF_YMIN(sf), SF_YSPACING(sf))
+
+	case SF_SPLINE1:
+	    call sf_evspline1 (COEFF(SF_COEFF(sf)), x, y, zfit, npts,
+	        SF_NXPIECES(sf), SF_NYPIECES(sf), -SF_XMIN(sf), SF_XSPACING(sf),
+		-SF_YMIN(sf), SF_YSPACING(sf))
+	}
+end
diff --git a/math/surfit/iszero.x b/math/surfit/iszero.x
new file mode 100644
index 00000000..d453d1fd
--- /dev/null
+++ b/math/surfit/iszero.x
@@ -0,0 +1,26 @@
+# Copyright(c) 1986 Association of Universities for Research in Astronomy Inc.
+
+include "surfitdef.h"
+
+# ISZERO -- Procedure to zero the accumulators for line number lineno.
+
+procedure iszero (sf)
+
+pointer	sf	# pointer to the surface descriptor
+
+begin
+	SF_NXPTS(sf) = 0
+	SF_NYPTS(sf) = 0
+
+	call aclrr (XMATRIX(SF_XMATRIX(sf)),
+		     SF_XORDER(sf) * SF_NXCOEFF(sf))
+	call aclrr (YMATRIX(SF_YMATRIX(sf)),
+		    SF_YORDER(sf) * SF_NYCOEFF(sf))
+	call aclrr (XCOEFF(SF_XCOEFF(sf)), SF_NXCOEFF(sf) * SF_NLINES(sf))
+	call aclrr (COEFF(SF_COEFF(sf)), SF_NXCOEFF(sf) * SF_NYCOEFF(sf))
+
+	if (SF_WZ(sf) != NULL)
+	    call mfree (SF_WZ(sf), MEM_TYPE)
+	if (SF_TLEFT(sf) != NULL)
+	    call mfree (SF_TLEFT(sf), TY_INT)
+end
diff --git a/math/surfit/mkpkg b/math/surfit/mkpkg
new file mode 100644
index 00000000..58948e29
--- /dev/null
+++ b/math/surfit/mkpkg
@@ -0,0 +1,29 @@
+# Surface fitting tools library.
+
+$checkout libsurfit.a lib$
+$update   libsurfit.a
+$checkin  libsurfit.a lib$
+$exit
+
+libsurfit.a:
+	iscoeff.x	surfitdef.h
+	iseval.x	surfitdef.h <math/surfit.h>
+	isfree.x	surfitdef.h
+	isinit.x	surfitdef.h <math/surfit.h>
+	islaccum.x	surfitdef.h <math/surfit.h>
+	islfit.x	surfitdef.h <math/surfit.h>
+	islrefit.x	surfitdef.h <math/surfit.h>
+	islsolve.x	surfitdef.h <math/surfit.h>
+	islzero.x	surfitdef.h
+	isreplace.x	surfitdef.h <math/surfit.h>
+	isresolve.x	surfitdef.h <math/surfit.h>
+	issave.x	surfitdef.h <math/surfit.h>
+	issolve.x	surfitdef.h <math/surfit.h>
+	isvector.x	surfitdef.h <math/surfit.h>
+	iszero.x	surfitdef.h
+	sf_b1eval.x	
+	sf_beval.x	
+	sf_f1deval.x	
+	sf_feval.x	
+	sfchomat.x	surfitdef.h <mach.h> <math/surfit.h>
+	;
diff --git a/math/surfit/sf_b1eval.x b/math/surfit/sf_b1eval.x
new file mode 100644
index 00000000..d07006fc
--- /dev/null
+++ b/math/surfit/sf_b1eval.x
@@ -0,0 +1,108 @@
+# Copyright(c) 1986 Association of Universities for Research in Astronomy Inc.
+
+# SF_B1LEG -- Procedure to evaluate all the non-zero Legendrefunctions for
+# a single point and given order.
+
+procedure sf_b1leg (x, order, k1, k2, basis)
+
+real	x		# array of data points
+int	order		# order of polynomial, order = 1, constant
+real	k1, k2		# normalizing constants
+real	basis[ARB]	# basis functions
+
+int	i
+real	ri, xnorm
+
+begin
+	basis[1] = 1.
+	if (order == 1)
+	    return
+
+	xnorm = (x + k1) * k2 
+	basis[2] = xnorm
+	if (order == 2)
+	    return
+
+	do i = 3, order {
+	    ri = i
+	    basis[i] = ((2. * ri - 3.) * xnorm * basis[i-1] -
+		       (ri - 2.) * basis[i-2]) / (ri - 1.)	
+	}
+end
+
+
+# SF_B1CHEB -- Procedure to evaluate all the non zero Chebyshev function
+# for a given x and order.
+
+procedure sf_b1cheb (x, order, k1, k2, basis)
+
+real	x		# number of data points
+int	order		# order of polynomial, 1 is a constant
+real	k1, k2		# normalizing constants
+real	basis[ARB]	# array of basis functions
+
+int	i
+real	xnorm
+
+begin
+	basis[1] = 1.
+	if (order == 1)
+	    return
+
+	xnorm = (x + k1) * k2
+	basis[2] = xnorm
+	if (order == 2)
+	    return
+
+	do i = 3, order
+	    basis[i] = 2. * xnorm * basis[i-1] - basis[i-2]
+end
+
+
+# SF_B1SPLINE1 -- Evaluate all the non-zero spline1 functions for a
+# single point.
+
+procedure sf_b1spline1 (x, npieces, k1, k2, basis, left)
+
+real	x		# set of data points
+int	npieces		# number of polynomial pieces minus 1
+real	k1, k2		# normalizing constants
+real	basis[ARB]	# basis functions
+int	left		# index of the appropriate spline functions
+
+real	xnorm
+
+begin
+	xnorm = (x + k1) * k2
+	left = min (int (xnorm), npieces)
+
+	basis[2] = xnorm - left
+	basis[1] = 1. - basis[2]
+end
+
+
+# SF_B1SPLINE3 --  Procedure to evaluate all the non-zero basis functions
+# for a cubic spline.
+
+procedure sf_b1spline3 (x, npieces, k1, k2, basis, left)
+
+real	x		# array of data points
+int	npieces		# number of polynomial pieces
+real	k1, k2		# normalizing constants
+real	basis[ARB]	# array of basis functions
+int	left		# array of indices for first non-zero spline
+
+real	sx, tx
+
+begin
+	sx = (x + k1) * k2
+	left = min (int (sx), npieces)
+
+	sx = sx - left
+	tx = 1. - sx
+
+	basis[1] = tx * tx * tx
+	basis[2] = 1. + tx * (3. + tx * (3. - 3. * tx))
+	basis[3] = 1. + sx * (3. + sx * (3. - 3. * sx))
+	basis[4] = sx * sx * sx
+end
diff --git a/math/surfit/sf_beval.x b/math/surfit/sf_beval.x
new file mode 100644
index 00000000..301967f9
--- /dev/null
+++ b/math/surfit/sf_beval.x
@@ -0,0 +1,143 @@
+# Copyright(c) 1986 Association of Universities for Research in Astronomy Inc.
+
+# SF_BCHEB -- Procedure to evaluate all the non-zero Chebyshev functions for
+# a set of points and given order.
+
+procedure sf_bcheb (x, npts, order, k1, k2, basis)
+
+real	x[npts]		# array of data points
+int	npts		# number of points
+int	order		# order of polynomial, order = 1, constant
+real	k1, k2		# normalizing constants
+real	basis[ARB]	# basis functions
+
+int	k, bptr
+
+begin
+	bptr = 1
+	do k = 1, order {
+
+	    if (k == 1)
+	        call amovkr (1., basis, npts)
+	    else if (k == 2)
+		call altar (x, basis[bptr], npts, k1, k2)
+	    else {
+		call amulr (basis[1+npts], basis[bptr-npts], basis[bptr],
+				npts)
+		call amulkr (basis[bptr], 2., basis[bptr], npts)
+		call asubr (basis[bptr], basis[bptr-2*npts], basis[bptr], npts)
+	    }
+		
+	    bptr = bptr + npts
+	}
+end
+
+
+# SF_BLEG -- Procedure to evaluate all the non zero Legendre function
+# for a given order and set of points.
+
+procedure sf_bleg (x, npts, order, k1, k2, basis)
+
+real	x[npts]		# number of data points
+int	npts		# number of points
+int	order		# order of polynomial, 1 is a constant
+real	k1, k2		# normalizing constants
+real	basis[ARB]	# array of basis functions
+
+int	k, bptr
+real	ri, ri1, ri2
+
+begin
+	bptr = 1
+
+	do k = 1, order {
+	    if (k == 1)
+		call amovkr (1., basis, npts)
+	    else if (k == 2)
+		call altar (x, basis[bptr], npts, k1, k2)
+	    else {
+		ri = k
+		ri1 = (2. * ri - 3.) / (ri - 1.)
+		ri2 = - (ri - 2.) / (ri - 1.)
+		call amulr (basis[1+npts], basis[bptr-npts], basis[bptr],
+				npts)
+		call awsur (basis[bptr], basis[bptr-2*npts],
+			basis[bptr], npts, ri1, ri2)
+	    }
+			
+	    bptr = bptr + npts
+	}
+end
+
+
+# SF_BSPLINE1 -- Evaluate all the non-zero spline1 functions for a set
+# of points.
+
+procedure sf_bspline1 (x, npts, npieces, k1, k2, basis, left)
+
+real	x[npts]		# set of data points
+int	npts		# number of points
+int	npieces		# number of polynomial pieces minus 1
+real	k1, k2		# normalizing constants
+real	basis[ARB]	# basis functions
+int	left[ARB]	# indices of the appropriate spline functions
+
+int	k
+
+begin
+	call altar (x, basis[1+npts], npts, k1, k2)
+	call achtri (basis[1+npts], left, npts)
+	call aminki (left, npieces, left, npts)
+
+	do k = 1, npts {
+	    basis[npts+k] = basis[npts+k] - left[k]
+	    basis[k] = 1. - basis[npts+k]
+	}
+end
+
+
+# SF_BSPLINE3 --  Procedure to evaluate all the non-zero basis functions
+# for a cubic spline.
+
+procedure sf_bspline3 (x, npts, npieces, k1, k2, basis, left)
+
+real	x[npts]		# array of data points
+int	npts		# number of data points
+int	npieces		# number of polynomial pieces minus 1
+real	k1, k2		# normalizing constants
+real	basis[ARB]	# array of basis functions
+int	left[ARB]	# array of indices for first non-zero spline
+
+int	i
+pointer	sp, sx, tx
+
+begin
+	# allocate space
+	call smark (sp)
+	call salloc (sx, npts, TY_REAL)
+	call salloc (tx, npts, TY_REAL)
+
+	# calculate the index of the first non-zero coeff
+	call altar (x, Memr[sx], npts, k1, k2)
+	call achtri (Memr[sx], left, npts)
+	call aminki (left, npieces, left, npts)
+
+	# normalize x to 0 to 1
+	do i = 1, npts {
+	    Memr[sx+i-1] = Memr[sx+i-1] - left[i]
+	    Memr[tx+i-1] = 1. - Memr[sx+i-1]
+	}
+
+	# calculate the basis function
+	call apowkr (Memr[tx], 3, basis, npts)
+	do i = 1, npts {
+	    basis[npts+i] = 1. + Memr[tx+i-1] * (3. + Memr[tx+i-1] * (3. -
+		3. * Memr[tx+i-1]))
+	    basis[2*npts+i] = 1. + Memr[sx+i-1] * (3. + Memr[sx+i-1] * (3. -
+		3. * Memr[sx+i-1]))
+	}
+	call apowkr (Memr[sx], 3, basis[1+3*npts], npts)
+
+	# release space
+	call sfree (sp)
+end
diff --git a/math/surfit/sf_f1deval.x b/math/surfit/sf_f1deval.x
new file mode 100644
index 00000000..01cb98d0
--- /dev/null
+++ b/math/surfit/sf_f1deval.x
@@ -0,0 +1,233 @@
+# Copyright(c) 1986 Association of Universities for Research in Astronomy Inc.
+
+# CV_EVCHEB -- Procedure to evaluate a Chebyshev polynomial assuming that
+# the coefficients have been calculated. 
+
+procedure cv_evcheb (coeff, x, yfit, npts, order, k1, k2)
+
+real	coeff[ARB]		# EV array of coefficients
+real	x[npts]			# x values of points to be evaluated
+real	yfit[npts]		# the fitted points
+int	npts			# number of points to be evaluated
+int	order			# order of the polynomial, 1 = constant
+real	k1, k2			# normalizing constants
+
+int	i
+pointer	sx, pn, pnm1, pnm2
+pointer sp
+real	c1, c2
+
+begin
+	# fit a constant
+	call amovkr (coeff[1], yfit, npts)
+	if (order == 1)
+	    return
+
+	# fit a linear function
+	c1 = k2 * coeff[2]
+	c2 = c1 * k1 + coeff[1]
+	call altmr (x, yfit, npts, c1, c2)
+	if (order == 2)
+	    return
+
+	# allocate temporary space
+	call smark (sp)
+	call salloc (sx, npts, TY_REAL)
+	call salloc (pn, npts, TY_REAL)
+	call salloc (pnm1, npts, TY_REAL)
+	call salloc (pnm2, npts, TY_REAL)
+
+	# a higher order polynomial
+	call amovkr (1., Memr[pnm2], npts)
+	call altar (x, Memr[sx], npts, k1, k2)
+	call amovr (Memr[sx], Memr[pnm1], npts)
+	call amulkr (Memr[sx], 2., Memr[sx], npts)
+	do i = 3, order {
+	    call amulr (Memr[sx], Memr[pnm1], Memr[pn], npts)
+	    call asubr (Memr[pn], Memr[pnm2], Memr[pn], npts)
+	    if (i < order) {
+	        call amovr (Memr[pnm1], Memr[pnm2], npts)
+	        call amovr (Memr[pn], Memr[pnm1], npts)
+	    }
+	    call amulkr (Memr[pn], coeff[i], Memr[pn], npts)
+	    call aaddr (yfit, Memr[pn], yfit, npts)
+	}
+
+	# free temporary space
+	call sfree (sp)
+
+end
+
+
+# CV_EVLEG -- Procedure to evaluate a Legendre polynomial assuming that
+# the coefficients have been calculated. 
+
+procedure cv_evleg (coeff, x, yfit, npts, order, k1, k2)
+
+real	coeff[ARB]		# EV array of coefficients
+real	x[npts]			# x values of points to be evaluated
+real	yfit[npts]		# the fitted points
+int	npts			# number of data points
+int	order			# order of the polynomial, 1 = constant
+real	k1, k2			# normalizing constants
+
+int	i
+pointer	sx, pn, pnm1, pnm2
+pointer	sp
+real	ri, ri1, ri2
+
+begin
+	# fit a constant
+	call amovkr (coeff[1], yfit, npts)
+	if (order == 1)
+	    return
+
+	# fit a linear function
+	ri1 = k2 * coeff[2]
+	ri2 = ri1 * k1 + coeff[1]
+	call altmr (x, yfit, npts, ri1, ri2)
+	if (order == 2)
+	    return
+
+	# allocate temporary space
+	call smark (sp)
+	call salloc (sx, npts, TY_REAL)
+	call salloc (pn, npts, TY_REAL)
+	call salloc (pnm1, npts, TY_REAL)
+	call salloc (pnm2, npts, TY_REAL)
+
+	# a higher order polynomial
+	call amovkr (1., Memr[pnm2], npts)
+	call altar (x, Memr[sx], npts, k1, k2)
+	call amovr (Memr[sx], Memr[pnm1], npts)
+	do i = 3, order {
+	    ri = i
+	    ri1 = (2. * ri - 3.) / (ri - 1.)
+	    ri2 = - (ri - 2.) / (ri - 1.)
+	    call amulr (Memr[sx], Memr[pnm1], Memr[pn], npts)
+	    call awsur (Memr[pn], Memr[pnm2], Memr[pn], npts, ri1, ri2)
+	    if (i < order) {
+	        call amovr (Memr[pnm1], Memr[pnm2], npts)
+	        call amovr (Memr[pn], Memr[pnm1], npts)
+	    }
+	    call amulkr (Memr[pn], coeff[i], Memr[pn], npts)
+	    call aaddr (yfit, Memr[pn], yfit, npts)
+	}
+
+	# free temporary space
+	call sfree (sp)
+
+end
+
+
+# CV_EVSPLINE1 -- Procedure to evaluate a piecewise linear spline function
+# assuming that the coefficients have been calculated.
+
+procedure cv_evspline1 (coeff, x, yfit, npts, npieces, k1, k2)
+
+real	coeff[ARB]		# array of coefficients
+real	x[npts]			# array of x values
+real	yfit[npts]		# array of fitted values
+int	npts			# number of data points
+int	npieces			# number of fitted points minus 1
+real	k1, k2			# normalizing constants
+
+int	j
+pointer sx, tx, index
+pointer	sp
+
+begin
+	# allocate the required space
+	call smark (sp)
+	call salloc (sx, npts, TY_REAL)
+	call salloc (tx, npts, TY_REAL)
+	call salloc (index, npts, TY_INT)
+
+	# calculate the index of the first non-zero coefficient
+	# for each point
+	call altar (x, Memr[sx], npts, k1, k2)
+	call achtri (Memr[sx], Memi[index], npts)
+	call aminki (Memi[index], npieces, Memi[index], npts)
+
+	# transform sx to range 0 to 1
+	do j = 1, npts {
+	    Memr[sx+j-1] = Memr[sx+j-1] - Memi[index+j-1]
+	    Memr[tx+j-1] = 1. - Memr[sx+j-1]
+	}
+
+    	# calculate yfit using the two non-zero basis function
+	call aclrr (yfit, npts)
+	do j = 1, npts
+	    yfit[j] = Memr[tx+j-1] * coeff[1+Memi[index+j-1]] +
+		      Memr[sx+j-1] * coeff[2+Memi[index+j-1]]
+
+        # free space
+	call sfree (sp)
+
+end
+
+
+# CV_EVSPLINE3 -- Procedure to evaluate the cubic spline assuming that
+# the coefficients of the fit are known.
+
+procedure cv_evspline3 (coeff, x, yfit, npts, npieces, k1, k2)
+
+real	coeff[ARB]	# array of coeffcients
+real	x[npts]		# array of x values
+real	yfit[npts]	# array of fitted values
+int	npts		# number of data points
+int	npieces		# number of polynomial pieces
+real	k1, k2		# normalizing constants
+
+int	i, j
+pointer	sx, tx, temp, index, sp
+
+begin
+
+	# allocate the required space
+	call smark (sp)
+        call salloc (sx, npts, TY_REAL)
+	call salloc (tx, npts, TY_REAL)
+	call salloc (temp, npts, TY_REAL)
+	call salloc (index, npts, TY_INT)
+
+	# calculate to which coefficients the x values contribute to
+        call altar (x, Memr[sx], npts, k1, k2)
+        call achtri (Memr[sx], Memi[index], npts)
+        call aminki (Memi[index], npieces, Memi[index], npts)
+
+        # transform sx to range 0 to 1
+	do j = 1, npts {
+	    Memr[sx+j-1] = Memr[sx+j-1] - Memi[index+j-1]
+	    Memr[tx+j-1] = 1. - Memr[sx+j-1]
+	}
+
+        # calculate yfit using the four non-zero basis function
+	call aclrr (yfit, npts)
+        do i = 1, 4 {
+
+	    switch (i) {
+	    case 1:
+		call apowkr (Memr[tx], 3, Memr[temp], npts)
+	    case 2:
+		do j = 1, npts {
+		    Memr[temp+j-1] = 1. + Memr[tx+j-1] * (3. + Memr[tx+j-1] *
+			(3. - 3. * Memr[tx+j-1]))
+		}
+	    case 3:
+		do j = 1, npts {
+		    Memr[temp+j-1] = 1. + Memr[sx+j-1] * (3. + Memr[sx+j-1] *
+			(3. - 3. * Memr[sx+j-1]))
+		}
+	    case 4:
+		call apowkr (Memr[sx], 3, Memr[temp], npts)
+	    }
+
+	    do j = 1, npts
+		Memr[temp+j-1] = Memr[temp+j-1] * coeff[i+Memi[index+j-1]]
+	    call aaddr (yfit, Memr[temp], yfit, npts)
+	}
+
+	# free space
+	call sfree (sp)
+end
diff --git a/math/surfit/sf_feval.x b/math/surfit/sf_feval.x
new file mode 100644
index 00000000..58b2f765
--- /dev/null
+++ b/math/surfit/sf_feval.x
@@ -0,0 +1,280 @@
+# Copyright(c) 1986 Association of Universities for Research in Astronomy Inc.
+
+# SF_EVCHEB -- Procedure to evaluate a Chebyshev polynomial assuming that
+# the coefficients have been calculated. 
+
+procedure sf_evcheb (coeff, x, y, zfit, npts, xterms, xorder, yorder, k1x, k2x,
+	k1y, k2y)
+
+real	coeff[ARB]		# 1D array of coefficients
+real	x[npts]			# x values of points to be evaluated
+real	y[npts]
+real	zfit[npts]		# the fitted points
+int	npts			# number of points to be evaluated
+int	xterms			# cross terms ?
+int	xorder,yorder		# order of the polynomials in x and y
+real	k1x, k2x		# normalizing constants
+real	k1y, k2y
+
+int	i, k, j
+int	ytorder, cptr
+pointer	sp
+pointer	xb, yb, accum
+pointer ybzptr, ybptr, xbzptr
+
+begin
+	# fit a constant
+	if (xorder == 1 && yorder == 1) {
+	    call amovkr (coeff[1], zfit, npts)
+	    return
+	}
+
+	# allocate temporary space for the basis functions
+	call smark (sp)
+	call salloc (xb, xorder * npts, TY_REAL)
+	call salloc (yb, yorder * npts, TY_REAL)
+	call salloc (accum, npts, TY_REAL)
+
+	# calculate basis functions
+	call sf_bcheb (x, npts, xorder, k1x, k2x, Memr[xb])
+	call sf_bcheb (y, npts, yorder, k1y, k2y, Memr[yb])
+
+	# clear the accumulator
+	call aclrr (zfit, npts)
+
+	# accumulate the values
+	cptr = 0
+	ybzptr = yb - 1
+	xbzptr = xb - 1
+	ytorder = yorder
+
+	do i = 1, xorder {
+	    call aclrr (Memr[accum], npts)
+	    ybptr = ybzptr
+	    do k = 1, ytorder {
+		do j = 1, npts
+		    Memr[accum+j-1] = Memr[accum+j-1] + coeff[cptr+k] *
+		        Memr[ybptr+j]
+		ybptr = ybptr + npts
+	    }
+	    do j = 1, npts
+		zfit[j] = zfit[j] + Memr[accum+j-1] * Memr[xbzptr+j]
+
+	    if (xterms == NO)
+		ytorder = 1
+
+	    cptr = cptr + yorder
+	    xbzptr = xbzptr + npts
+	}
+
+	# free temporary space
+	call sfree (sp)
+end
+
+
+# SF_EVLEG -- Procedure to evaluate a Chebyshev polynomial assuming that
+# the coefficients have been calculated. 
+
+procedure sf_evleg (coeff, x, y, zfit, npts, xterms, xorder, yorder, k1x, k2x,
+	k1y, k2y)
+
+real	coeff[ARB]		# 1D array of coefficients
+real	x[npts]			# x values of points to be evaluated
+real	y[npts]
+real	zfit[npts]		# the fitted points
+int	npts			# number of points to be evaluated
+int	xterms			# cross terms ?
+int	xorder,yorder		# order of the polynomials in x and y
+real	k1x, k2x		# normalizing constants
+real	k1y, k2y
+
+int	i, k, j
+int	ytorder, cptr
+pointer	sp
+pointer	xb, yb, accum
+pointer ybzptr, ybptr, xbzptr
+
+begin
+	# fit a constant
+	if (xorder == 1 && yorder == 1) {
+	    call amovkr (coeff[1], zfit, npts)
+	    return
+	}
+
+	# allocate temporary space for the basis functions
+	call smark (sp)
+	call salloc (xb, xorder * npts, TY_REAL)
+	call salloc (yb, yorder * npts, TY_REAL)
+	call salloc (accum, npts, TY_REAL)
+
+	# calculate basis functions
+	call sf_bleg (x, npts, xorder, k1x, k2x, Memr[xb])
+	call sf_bleg (y, npts, yorder, k1y, k2y, Memr[yb])
+
+	# clear the accumulator
+	call aclrr (zfit, npts)
+
+	# accumulate the values
+	cptr = 0
+	ybzptr = yb - 1
+	xbzptr = xb - 1
+	ytorder = yorder
+
+	do i = 1, xorder {
+	    call aclrr (Memr[accum], npts)
+	    ybptr = ybzptr
+	    do k = 1, ytorder {
+		do j = 1, npts
+		    Memr[accum+j-1] = Memr[accum+j-1] + coeff[cptr+k] *
+		        Memr[ybptr+j]
+		ybptr = ybptr + npts
+	    }
+	    do j = 1, npts
+		zfit[j] = zfit[j] + Memr[accum+j-1] * Memr[xbzptr+j]
+
+	    if (xterms == NO)
+		ytorder = 1
+
+	    cptr = cptr + yorder
+	    xbzptr = xbzptr + npts
+	}
+
+	# free temporary space
+	call sfree (sp)
+end
+
+
+# SF_EVSPLINE3 -- Procedure to evaluate a piecewise linear spline function
+# assuming that the coefficients have been calculated.
+
+procedure sf_evspline3 (coeff, x, y, zfit, npts, nxpieces, nypieces, k1x, k2x,
+	k1y, k2y)
+
+real	coeff[ARB]		# array of coefficients
+real	x[npts]			# array of x values
+real	y[npts]			# array of y values
+real	zfit[npts]		# array of fitted values
+int	npts			# number of data points
+int	nxpieces, nypieces	# number of fitted points minus 1
+real	k1x, k2x		# normalizing constants
+real	k1y, k2y
+
+int	i, j, k, cindex
+pointer	xb, xbzptr, yb, ybzptr, ybptr
+pointer	accum, leftx, lefty
+pointer	sp
+
+begin
+	# allocate temporary space for the basis functions
+	call smark (sp)
+	call salloc (xb, 4 * npts, TY_REAL)
+	call salloc (yb, 4 * npts, TY_REAL)
+	call salloc (accum, npts, TY_REAL)
+	call salloc (leftx, npts, TY_INT)
+	call salloc (lefty, npts, TY_INT)
+
+	# calculate basis functions
+	call sf_bspline3 (x, npts, nxpieces, k1x, k2x, Memr[xb], Memi[leftx])
+	call sf_bspline3 (y, npts, nypieces, k1y, k2y, Memr[yb], Memi[lefty])
+
+	# set up the indexing
+	call amulki (Memi[leftx], (nypieces+4), Memi[leftx], npts)
+	call aaddi (Memi[leftx], Memi[lefty], Memi[lefty], npts)
+
+	# clear the accumulator
+	call aclrr (zfit, npts)
+
+	# accumulate the values
+
+	ybzptr = yb - 1
+	xbzptr = xb - 1
+
+	do i = 1, 4 {
+	    call aclrr (Memr[accum], npts)
+	    ybptr = ybzptr
+	    do k = 1, 4 {
+		do j = 1, npts {
+		    cindex = k + Memi[lefty+j-1]
+		    Memr[accum+j-1] = Memr[accum+j-1] + coeff[cindex] *
+		        Memr[ybptr+j]
+		}
+		ybptr = ybptr + npts
+	    }
+	    do j = 1, npts
+		zfit[j] = zfit[j] + Memr[accum+j-1] * Memr[xbzptr+j]
+
+	    xbzptr = xbzptr + npts
+	    call aaddki (Memi[lefty], (nypieces+4), Memi[lefty], npts)
+	}
+
+	# free temporary space
+	call sfree (sp)
+end
+
+
+# SF_EVSPLINE1 -- Procedure to evaluate a piecewise linear spline function
+# assuming that the coefficients have been calculated.
+
+procedure sf_evspline1 (coeff, x, y, zfit, npts, nxpieces, nypieces, k1x, k2x,
+	k1y, k2y)
+
+real	coeff[ARB]		# array of coefficients
+real	x[npts]			# array of x values
+real	y[npts]			# array of y values
+real	zfit[npts]		# array of fitted values
+int	npts			# number of data points
+int	nxpieces, nypieces	# number of fitted points minus 1
+real	k1x, k2x		# normalizing constants
+real	k1y, k2y
+
+int	i, j, k, cindex
+pointer	xb, xbzptr, yb, ybzptr, ybptr
+pointer	accum, leftx, lefty
+pointer	sp
+
+begin
+	# allocate temporary space for the basis functions
+	call smark (sp)
+	call salloc (xb, 2 * npts, TY_REAL)
+	call salloc (yb, 2 * npts, TY_REAL)
+	call salloc (accum, npts, TY_REAL)
+	call salloc (leftx, npts, TY_INT)
+	call salloc (lefty, npts, TY_INT)
+
+	# calculate basis functions
+	call sf_bspline1 (x, npts, nxpieces, k1x, k2x, Memr[xb], Memi[leftx])
+	call sf_bspline1 (y, npts, nypieces, k1y, k2y, Memr[yb], Memi[lefty])
+
+	# set up the indexing
+	call amulki (Memi[leftx], (nypieces+2), Memi[leftx], npts)
+	call aaddi (Memi[leftx], Memi[lefty], Memi[lefty], npts)
+
+	# clear the accumulator
+	call aclrr (zfit, npts)
+
+	# accumulate the values
+
+	ybzptr = yb - 1
+	xbzptr = xb - 1
+
+	do i = 1, 2 {
+	    call aclrr (Memr[accum], npts)
+	    ybptr = ybzptr
+	    do k = 1, 2 {
+		do j = 1, npts {
+		    cindex = k + Memi[lefty+j-1]
+		    Memr[accum+j-1] = Memr[accum+j-1] + coeff[cindex] *
+		        Memr[ybptr+j]
+		}
+		ybptr = ybptr + npts
+	    }
+	    do j = 1, npts
+		zfit[j] = zfit[j] + Memr[accum+j-1] * Memr[xbzptr+j]
+
+	    xbzptr = xbzptr + npts
+	    call aaddki (Memi[lefty], (nypieces+2), Memi[lefty], npts)
+	}
+
+	# free temporary space
+	call sfree (sp)
+end
diff --git a/math/surfit/sfchomat.x b/math/surfit/sfchomat.x
new file mode 100644
index 00000000..419fe777
--- /dev/null
+++ b/math/surfit/sfchomat.x
@@ -0,0 +1,105 @@
+# Copyright(c) 1986 Association of Universities for Research in Astronomy Inc.
+
+include <math/surfit.h>
+include "surfitdef.h"
+include <mach.h>
+
+# SFCHOFAC -- Routine to calculate the Cholesky factorization of a
+# symmetric, positive semi-definite banded matrix. This routines was
+# adapted from the bchfac.f routine described in "A Practical Guide
+# to Splines", Carl de Boor (1978).
+
+procedure sfchofac (matrix, nbands, nrows, matfac, ier)
+
+VAR_TYPE matrix[nbands, nrows]	# data matrix
+int	nbands			# number of bands
+int	nrows			# number of rows
+VAR_TYPE matfac[nbands, nrows]	# Cholesky factorization
+int	ier			# error code
+
+int	i, n, j, imax, jmax
+VAR_TYPE ratio
+
+begin
+	if (nrows == 1) {
+	    if (matrix[1,1] > 0.)
+	        matfac[1,1] = 1. / matrix[1,1]
+	    return
+	}
+		
+	# copy matrix into matfac
+	do n = 1, nrows {
+	    do j = 1, nbands
+		matfac[j,n] = matrix[j,n]
+	}
+
+	do n = 1, nrows {
+	    # test to see if matrix is singular
+	    if (((matfac[1,n] + matrix[1,n]) - matrix[1,n]) <= DELTA) {
+		do j = 1, nbands
+		    matfac[j,n] = 0.
+		ier = SINGULAR
+		next
+	    }
+
+	    matfac[1,n] = 1. / matfac[1,n]
+	    imax = min (nbands - 1, nrows - n)
+	    if (imax < 1)
+		next
+
+	    jmax = imax
+	    do i = 1, imax {
+		ratio = matfac[i+1,n] * matfac[1,n]
+		do j = 1, jmax
+		    matfac[j,n+i] = matfac[j,n+i] - matfac[j+i,n] * ratio
+		jmax = jmax - 1
+		matfac[i+1,n] = ratio
+	    }
+	}
+end
+
+
+# SFCHOSLV -- Solve the matrix whose Cholesky factorization was calculated in
+# SFCHOFAC for the coefficients. This routine was adapted from bchslv.f
+# described in "A Practical Guide to Splines", by Carl de Boor (1978).
+
+procedure sfchoslv (matfac, nbands, nrows, vector, coeff)
+
+VAR_TYPE matfac[nbands,nrows] 		# Cholesky factorization
+int	nbands				# number of bands
+int	nrows				# number of rows
+VAR_TYPE vector[nrows]			# right side of matrix equation
+VAR_TYPE coeff[nrows]			# coefficients
+
+int	i, n, j, jmax, nbndm1
+
+begin
+	if (nrows == 1) {
+	    coeff[1] = vector[1] * matfac[1,1]
+	    return
+	}
+
+	# copy vector to coefficients
+	do i = 1, nrows
+	    coeff[i] = vector[i]
+
+	# forward substitution
+	nbndm1 = nbands - 1
+	do n = 1, nrows {
+	    jmax = min (nbndm1, nrows - n)
+	    if (jmax >= 1) {
+	        do j = 1, jmax
+		    coeff[j+n] = coeff[j+n] - matfac[j+1,n] * coeff[n]
+	    }
+	}
+
+	# back substitution
+	for (n = nrows; n >= 1; n = n - 1) {
+	    coeff[n] = coeff[n] * matfac[1,n]
+	    jmax = min (nbndm1, nrows - n)
+	    if (jmax >= 1) {
+		do j = 1, jmax
+		    coeff[n] = coeff[n] - matfac[j+1,n] * coeff[j+n]
+	    }
+	}
+end
diff --git a/math/surfit/surfitdef.h b/math/surfit/surfitdef.h
new file mode 100644
index 00000000..846109ea
--- /dev/null
+++ b/math/surfit/surfitdef.h
@@ -0,0 +1,74 @@
+# set up the curve fitting structure
+
+define	LEN_SFSTRUCT	 35
+
+define	SF_TYPE		Memi[$1]	# Type of curve to be fitted
+define	SF_NXCOEFF	Memi[$1+1]	# Number of coefficients
+define	SF_XORDER	Memi[$1+2]	# Order of the fit in x
+define	SF_NXPIECES	Memi[$1+3]	# Number of x polynomial pieces - 1
+define	SF_NCOLS	Memi[$1+4]	# Maximum x value
+define	SF_NXPTS	Memi[$1+5]	# Number of points in x
+define	SF_NYCOEFF	Memi[$1+6]	# Number of y coefficients
+define	SF_YORDER	Memi[$1+7]	# Order of the fit in y
+define	SF_NYPIECES	Memi[$1+8]	# Number of y polynomial pieces - 1
+define	SF_NLINES	Memi[$1+9]	# Minimum x value
+define	SF_NYPTS	Memi[$1+10]	# Number of y points
+define	SF_XTERMS	Memi[$1+11]	# cross terms?
+
+define	SF_XBASIS	Memi[$1+12]	# Pointer to the x basis functions
+define	SF_XLEFT	Memi[$1+13]	# Indices to x basis functions, spline
+define	SF_YBASIS	Memi[$1+14]	# Pointer to the y basis functions
+define	SF_YLEFT	Memi[$1+15]	# Indices to y basis functions, spline
+define	SF_XMATRIX	Memi[$1+16]	# Pointer to x data matrix
+define	SF_YMATRIX	Memi[$1+17]	# Pointer to y data matrix
+define	SF_XCOEFF	Memi[$1+18]	# X coefficient matrix
+define	SF_COEFF	Memi[$1+19]	# Pointer to coefficient vector
+
+define	SF_XMIN		Memr[P2R($1+20)] # Min x value
+define	SF_XMAX		Memr[P2R($1+21)] # Max x value
+define	SF_XRANGE	Memr[P2R($1+22)] # 2. / (xmax - xmin), polynomials
+define	SF_XMAXMIN	Memr[P2R($1+23)] # - (xmax + xmin) / 2., polynomials
+define	SF_XSPACING	Memr[P2R($1+24)] # order / (xmax - xmin), splines
+define	SF_YMIN		Memr[P2R($1+25)] # Min y value
+define	SF_YMAX		Memr[P2R($1+26)] # Max y value
+define	SF_YRANGE	Memr[P2R($1+27)] # 2. / (ymax - ymin), polynomials
+define	SF_YMAXMIN	Memr[P2R($1+28)] # - (ymax + ymin) / 2., polynomials
+define	SF_YSPACING	Memr[P2R($1+29)] # order / (ymax - ymin), splines
+
+define	SF_WZ		Memi[$1+30]
+define	SF_TLEFT	Memi[$1+31]
+
+# matrix and vector element definitions
+
+define	XBASIS		Memr[P2P($1)]	#
+define	XBS		Memr[P2P($1)]	#
+define	YBASIS		Memr[P2P($1)]	#
+define	YBS		Memr[P2P($1)]	#
+define	XMATRIX		Memr[P2P($1)]	#
+define	XCHOFAC		Memr[P2P($1)]	# 
+define	YMATRIX		Memr[P2P($1)]	#
+define	XCOEFF		Memr[P2P($1)]	#
+define	COEFF		Memr[P2P($1)]	#
+define	XLEFT		Memi[$1]	#
+define	YLEFT		Memi[$1]	#
+
+# structure definitions for the save restore functions
+
+define	SF_SAVETYPE	$1[1]
+define	SF_SAVEXORDER	$1[2]
+define	SF_SAVEYORDER	$1[3]
+define	SF_SAVEXTERMS	$1[4]
+define	SF_SAVENCOLS	$1[5]
+define	SF_SAVENLINES	$1[6]
+define	SF_SAVECOEFF	6
+
+# miscellaneous
+
+define	SPLINE3_ORDER	4
+define	SPLINE1_ORDER	2
+
+# data type
+
+define	MEM_TYPE	TY_REAL
+define	VAR_TYPE	real
+define	DELTA		EPSILON