Initial commit

26 files changed, 1177 insertions, 0 deletions

diff --git a/math/gsurfit/doc/gsaccum.hlp b/math/gsurfit/doc/gsaccum.hlp
new file mode 100644
index 00000000..afa63f70
--- /dev/null
+++ b/math/gsurfit/doc/gsaccum.hlp
@@ -0,0 +1,51 @@
+.help gsaccum Aug85 "Gsurfit Package"
+.ih
+NAME
+gsaccum -- accumulate a single data point into the fit
+.ih
+SYNOPSIS
+include <math/gsurfit.h>
+
+gsaccum (sf, x, y, weight, wtflag)
+
+.nf
+pointer	sf		# surface descriptor
+real	x		# x value, xmin <= x <= xmax
+real	y		# y value, ymin <= y <= ymax
+real	z		# z value
+real	weight		# weight
+int	wtflag		# type of weighting
+.fi
+.ih
+ARGUMENTS
+.ls sf      
+Pointer to the surface descriptor structure.
+.le
+.ls x, y
+The x and y values.
+.le
+.ls z      
+Data value.
+.le
+.ls weight
+The weight assigned to the data point.
+.le
+.ls wtflag
+Type of weighting. The options are WTS_USER and WTS_UNIFORM. If wtflag
+equals WTS_USER the weight for each data point is supplied by the user.
+If wtflag equals WTS_UNIFORM the routine sets the weight to 1.
+.le
+.ih
+DESCRIPTION
+GSACCUM calculates the non-zero basis functions for the given x and
+y values, computes the contribution of each data point to the normal
+equations and sums that contribution into the appropriate arrays and
+vectors.
+.ih
+NOTES
+Checking for out of bounds x and y values and INDEF valued data points is
+the responsibility of the calling program.
+.ih
+SEE ALSO
+gsacpts, gsfit, gsrefit
+.endhelp
diff --git a/math/gsurfit/doc/gsacpts.hlp b/math/gsurfit/doc/gsacpts.hlp
new file mode 100644
index 00000000..1e253c61
--- /dev/null
+++ b/math/gsurfit/doc/gsacpts.hlp
@@ -0,0 +1,56 @@
+.help gsacpts Aug85 "Gsurfit Package"
+.ih
+NAME
+gsacpts -- accumulate an array of data points into the fit
+.ih
+SYNOPSIS
+include <math/gsurfit.h>
+
+gsacpts (sf, x, y, z, weight, npts, wtflag)
+
+.nf
+pointer	sf		# surface descriptor
+real	x[npts]		# x values, xmin <= x <= xmax
+real	y[npts]		# y values, ymin <= y <= ymax
+real	z[npts]		# z values
+real	weight[npts]	# array of weights
+int	npts		# the number of data points
+int	wtflag		# type of weighting
+.fi
+.ih
+ARGUMENTS
+.ls sf      
+Pointer to the surface descriptor structure.
+.le
+.ls x, y
+Array of x and y values.
+.le
+.ls z      
+Array of data values.
+.le
+.ls weight
+The weights assigned to the data points.
+.le
+.ls npts
+The number of data points.
+.le
+.ls wtflag
+Type of weighting. The options are WTS_USER and WTS_UNIFORM. If wtflag
+equals WTS_USER the weight for each data point is supplied by the user.
+If wtflag equals WTS_UNIFORM the routine sets the weight to 1.
+The weight definitions are contained in the package header file gsurfit.h.
+.le
+.ih
+DESCRIPTION
+GSACCUM calculates the non-zero basis functions for the given x and
+y values, computes the contribution of each data point to the normal
+equations and sums that contribution into the appropriate arrays and
+vectors.
+.ih
+NOTES
+Checking for out of bounds x and y values and INDEF valued data points is
+the responsibility of the calling program.
+.ih
+SEE ALSO
+gsaccum, gsfit, gsrefit
+.endhelp
diff --git a/math/gsurfit/doc/gsadd.hlp b/math/gsurfit/doc/gsadd.hlp
new file mode 100644
index 00000000..84d388fd
--- /dev/null
+++ b/math/gsurfit/doc/gsadd.hlp
@@ -0,0 +1,35 @@
+.help gsadd Aug85 "Gsurfit Package"
+.ih
+NAME
+gsadd -- add two surface fits together
+.ih
+SYNOPSIS
+gsadd (sf1, sf2, sf3)
+
+.nf
+pointer	sf1		# first surface descriptor
+pointer	sf2		# second surface descriptor
+pointer	sf3		# resultant surface descriptor
+.fi
+.ih
+ARGUMENTS
+.ls sf1    
+Pointer to the first surface descriptor.
+.le
+.ls sf2    
+Pointer to the second surface descriptor.
+.le
+.ls sf3    
+Pointer to the resultant surface descriptor.
+.le
+.ih
+DESCRIPTION
+The coefficients of the two surfaces are added together. GSADD checks
+that the curve_types are the same and that the fits are normalized over
+the same range of data.
+.ih
+NOTES
+.ih
+SEE ALSO
+gscopy, gssub
+.endhelp
diff --git a/math/gsurfit/doc/gscoeff.hlp b/math/gsurfit/doc/gscoeff.hlp
new file mode 100644
index 00000000..e4b792db
--- /dev/null
+++ b/math/gsurfit/doc/gscoeff.hlp
@@ -0,0 +1,39 @@
+.help gscoeff Aug85 "Gsurfit Package"
+.ih
+NAME
+gscoeff - get the number and values of the coefficients
+.ih
+SYNOPSIS
+gscoeff (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. Ncoeff may be obtained by a call
+to gsstati.
+.le
+
+.nf
+	ncoeff = gsstati (sf, GSNCOEFF)
+.fi
+.ih
+DESCRIPTION
+GSCOEFF fetches the coefficient array and the number of coefficients from
+the surface descriptor structure.
+.ih
+NOTES
+.ih
+SEE ALSO
+gserrors
+.endhelp
diff --git a/math/gsurfit/doc/gscopy.hlp b/math/gsurfit/doc/gscopy.hlp
new file mode 100644
index 00000000..46a0935f
--- /dev/null
+++ b/math/gsurfit/doc/gscopy.hlp
@@ -0,0 +1,32 @@
+.help gscopy Aug85 "Gsurfit Package"
+.ih
+NAME
+gscopy -- copy a surface fit
+.ih
+SYNOPSIS
+gscopy (sf1, sf2)
+
+.nf
+pointer	sf1		# old surface descriptor
+pointer	sf2		# new surface descriptor
+.fi
+.ih
+ARGUMENTS
+.ls sf1    
+Pointer to the old surface descriptor structure.
+.le
+.ls sf2   
+Pointer to the new surface descriptor structure.
+.le
+.ih
+DESCRIPTION
+The surface fit and parameters are copied for later use by
+GSEVAL or GSVECTOR.
+.ih
+NOTES
+The matrices and vectors used by the numerical fitting routines are not
+stored.
+.ih
+SEE ALSO
+gsadd, gssub
+.endhelp
diff --git a/math/gsurfit/doc/gsder.hlp b/math/gsurfit/doc/gsder.hlp
new file mode 100644
index 00000000..e1af66e9
--- /dev/null
+++ b/math/gsurfit/doc/gsder.hlp
@@ -0,0 +1,48 @@
+.help gsder Aug85 "Gsurfit Package"
+.ih
+NAME
+gsder -- evaluate the derivatives of the fitted surface
+.ih
+SYNOPSIS
+gsder (sf, x, y, zfit, npts, nxder, nyder)
+
+.nf
+pointer	sf		# surface descriptor
+real	x[npts]		# x array, xmin <= x[i] <= xmax
+real	y[npts]		# y array, ymin <= x[i] <= ymax
+real	zfit[npts]	# data values
+int	npts		# number of data points
+int	nxder		# order of x derivative, 0 = function
+int	nyder		# order of y derivative, 0 = function
+.fi
+.ih
+ARGUMENTS
+.ls sf    
+Pointer to the surface descriptor structure.
+.le
+.ls x, y
+Array of x and y values.
+.le
+.ls zfit
+Array of fitted values.
+.le
+.ls npts
+The number of points to be fit.
+.le
+.ls nxder, nyder
+The order of derivative to be fit. GSDER is the same as GSVECTOR if nxder = 0
+and nyder = 0. If nxder = 1 and nyder = 0 GSDER calculates the first
+derivatives of the surface with respect to x.
+.le
+.ih
+DESCRIPTION
+Evaluate the derivatives of a surface at a set of data points.
+GSDER 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
+gseval, gsvector
+.endhelp
diff --git a/math/gsurfit/doc/gserrors.hlp b/math/gsurfit/doc/gserrors.hlp
new file mode 100644
index 00000000..fed9a82e
--- /dev/null
+++ b/math/gsurfit/doc/gserrors.hlp
@@ -0,0 +1,61 @@
+.help gserrors Aug85 "Gsurfit Package"
+.ih
+NAME
+.nf
+gserrors -- calculate errors of the coefficients and the chi-square
+	    of the fit
+.fi
+.ih
+SYNOPSIS
+gserrors (sf, y, weight, yfit, chi_square, errors)
+
+.nf
+pointer	sf		# surface descriptor
+real	y[ARB]		# array of data values
+real	weight[ARB]	# array of weights
+real	yfit[ARB]	# array of fitted values
+real	chi_square	# chi_square of fit
+real	errors[ARB]	# array of errors
+.fi
+.ih
+ARGUMENTS
+.ls sf
+Pointer to the surface descriptor structure.
+.le
+.ls y
+Array of data values.
+.le
+.ls weight
+Array of weights.
+.le
+.ls yfit
+Array of fitted values.
+.le
+.ls chi_square
+The reduced chi-square of the fit.
+.le
+.ls errors
+The array of errors of the coefficients. The number of coefficients
+can be obtained by a call to gsstati.
+.le
+
+.nf
+	    nerrors = gsstati (sf, GSNCOEFF)
+.fi
+.ih
+DESCRIPTION
+GSCOEFF calculates the reduced chi-square of the fit and the standard
+deviation of the coefficients.
+The chi-square of the fit is the square root of the sum of the
+weighted squares of the residuals divided by the number of degrees
+of freedom. If the weights are equal, then the reduced chi-square is
+the variance of the fit. The error of the j-th coefficient is the
+square root of the j-th diagonal element of the inverse of the data
+matrix. If the weights are equal to one, then the errors are scaled
+by the square root of the variance of the data.
+.ih
+NOTES
+.ih
+SEE ALSO
+gscoeff
+.endhelp
diff --git a/math/gsurfit/doc/gseval.hlp b/math/gsurfit/doc/gseval.hlp
new file mode 100644
index 00000000..b9cd08bb
--- /dev/null
+++ b/math/gsurfit/doc/gseval.hlp
@@ -0,0 +1,34 @@
+.help gseval Aug85 "Gsurfit Package"
+.ih
+NAME
+gseval -- evaluate the fitted surface at x and y
+.ih
+SYNOPSIS
+y = gseval (sf, x, y)
+
+.nf
+pointer	sf		# surface descriptor
+real	x		# x value, xmin <= x <= xmax
+real	y		# y value, ymin <= y <= ymax
+.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. GSEVAL is a real
+valued function which returns the fitted value.
+.ih
+NOTES
+GSEVAL 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
+gsvector, gsder
+.endhelp
diff --git a/math/gsurfit/doc/gsfit.hlp b/math/gsurfit/doc/gsfit.hlp
new file mode 100644
index 00000000..4abdc546
--- /dev/null
+++ b/math/gsurfit/doc/gsfit.hlp
@@ -0,0 +1,64 @@
+.help gsfit Aug85 "Gsurfit Package"
+.ih
+NAME
+gsfit -- fit a surface to a set of data values
+.ih
+SYNOPSIS
+include <math/gsurfit.h>
+
+gsfit (sf, x, y, z, weight, npts, wtflag, ier)
+
+.nf
+pointer	sf		# surface descriptor
+real	x[npts]		# x array, xmin <= x[i] <= xmax
+real	y[npts]		# y array, ymin <= y[i] <= ymax
+real	z[npts]		# data values
+real	weight[npts]	# weight array
+int	npts		# number of data points
+int	wtflag		# type of weighting
+int	ier		# error coded
+.fi
+.ih
+ARGUMENTS
+.ls sf    
+Pointer to the surface descriptor structure.
+.le
+.ls x, y
+X and y value arrays.
+.le
+.ls z      
+Array of data values.
+.le
+.ls weight
+Array of weights.
+.le
+.ls npts
+Number of data points
+.le
+.ls wtflag
+Type of weighting. The options are WTS_USER and WTS_UNIFORM. If wtflag =
+WTS_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 = WTS_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
+GSFIT zeroes the matrix and vectors, calculates the non-zero basis functions,
+computes the contribution of each data point to the normal equations
+and accumulates it into the appropriate array and vector elements. The
+Cholesky factorization of the coefficient array is computed and the coefficients
+of the fitting function are calculated.
+.ih
+NOTES
+Checking for out of bounds x and y values is the responsibility of the user.
+.ih
+SEE ALSO
+gsrefit, gsaccum, gsacpts, gssolve, gszero
+.endhelp
diff --git a/math/gsurfit/doc/gsfree.hlp b/math/gsurfit/doc/gsfree.hlp
new file mode 100644
index 00000000..a576b2e1
--- /dev/null
+++ b/math/gsurfit/doc/gsfree.hlp
@@ -0,0 +1,26 @@
+.help gsfree Aug85 "Gsurfit Package"
+.ih
+NAME
+gsfree -- free the surface descriptor structure
+.ih
+SYNOPSIS
+gsfree (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.
+.ih
+NOTES
+GSFREE should be called after each surface fit.
+.ih
+SEE ALSO
+gsinit
+.endhelp
diff --git a/math/gsurfit/doc/gsgcoeff.hlp b/math/gsurfit/doc/gsgcoeff.hlp
new file mode 100644
index 00000000..78fdc707
--- /dev/null
+++ b/math/gsurfit/doc/gsgcoeff.hlp
@@ -0,0 +1,31 @@
+.help gsgcoeff Aug85 "Gsurfit Package"
+.ih
+NAME
+gsgcoeff -- Procedure to fetch a coefficient
+.ih
+SYNOPSIS
+rval = gsgcoeff (sf, xorder, yorder)
+
+.nf
+    pointer	sf	# surface descriptor	
+    int		xorder	# x order of desired coefficient
+    int		yorder	# y order of desired coefficient
+.fi
+.ih
+ARGUMENTS
+.ls sf    
+Pointer to the surface descriptor.
+.le
+.ls xorder, yorder
+The x and y order of the desired coefficient.
+.le
+.ih
+DESCRIPTION
+GSGCOEFF fetches the coefficient of x ** (xorder - 1) * y ** (yorder - 1).
+INDEF is returned if xorder and yorder are out of range.
+.ih
+NOTES
+.ih
+SEE ALSO
+gsscoeff
+.endhelp
diff --git a/math/gsurfit/doc/gsinit.hlp b/math/gsurfit/doc/gsinit.hlp
new file mode 100644
index 00000000..3a647a8c
--- /dev/null
+++ b/math/gsurfit/doc/gsinit.hlp
@@ -0,0 +1,64 @@
+.help gsinit Aug85 "Gsurfit Package"
+.ih
+NAME
+gsinit -- initialize surface descriptor
+.ih
+SYNOPSIS
+include <math/gsurfit.h>
+
+.nf
+gsinit (sf, surface_type, xorder, yorder, xterms, xmin, xmax,
+	ymin, ymax)
+.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)
+real	xmin		# minimum x value
+real	xmax		# maximum x value
+real	ymin		# minimum y value
+real	ymax		# maximum y value
+.fi
+.ih
+ARGUMENTS
+.ls sf    
+Pointer to the surface descriptor structure.
+.le
+.ls surface_type
+Fitting function. Permitted values are GS_LEGENDRE and GS_CHEBYSHEV for
+Legendre and Chebyshev polynomials.
+.le
+.ls xorder, yorder
+Order of the polynomial to be fit. The order must be greater than or
+equal to 1. If xorder == 1 and yorder == 1 a constant is fit to the data.
+.le
+.ls xterms
+Set the cross-terms type? The options are GS_XNONE (the old NO option) for
+no cross terms, GS_XHALF for diagonal cross terms (new option), and GS_XFULL
+for full cross terms (the old YES option).
+.le
+.ls xmin, xmax
+Minimum and maximum x values. All the x values of interest including the
+data x values and the x values of any surface to be evaluated must
+fall in the range xmin <= x <= xmax.
+.le
+.ls ymin, ymax
+Minimum and maximum y values. All the y values of interest including the
+data y values and the y values of any surface to be evaluated must
+fall in the range ymin <= y <= ymax.
+.le
+.ih
+DESCRIPTION
+GSINIT allocates space for the surface descriptor and the arrays and vectors
+used by the numerical routines. It initializes all arrays and vectors to zero
+and returns the surface descriptor to the calling routine.
+.ih
+NOTES
+GSINIT must be the first GSURFIT routine called.
+.ih
+SEE ALSO
+gsfree
+.endhelp
diff --git a/math/gsurfit/doc/gsrefit.hlp b/math/gsurfit/doc/gsrefit.hlp
new file mode 100644
index 00000000..629697e8
--- /dev/null
+++ b/math/gsurfit/doc/gsrefit.hlp
@@ -0,0 +1,55 @@
+.help gsrefit Aug85 "Gsurfit Package"
+.ih
+NAME
+gsrefit -- refit with new z vector using old x, y and weight vector
+.ih
+SYNOPSIS
+include < math/gsurfit.h>
+
+gsrefit (sf, x, y, z, w, ier)
+
+.nf
+pointer	sf		# surface descriptor
+real	x[ARB]		# x array, xmin <= x[i] <= xmax
+real	y[ARB]		# y array, ymin <= y[i] <= ymax
+real	z[ARB]		# array of data values
+real	w[ARB]		# array of weights
+int	ier		# error code
+.fi
+.ih
+ARGUMENTS
+.ls sf          
+Pointer to the surface descriptor structure.
+.le
+.ls x, y
+Array of x and y values.
+.le
+.ls z    
+Array of data values.
+.le
+.ls w    
+Array of weights.
+.le
+.ls ier    
+Error code. The options are OK, SINGULAR and NO_DEG_FREEDOM. If ier =
+SINGULAR a solution is computed but one or more coefficients may be zero.
+If ier equals NO_DEG_FREEDOM, there are insufficient data points to
+compute a solution and GSREFIT returns without solving for the coefficients.
+.le
+.ih
+DESCRIPTION
+In some applications the x, y and weight 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. GSREFIT
+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
+In the first call to GSREFIT space is allocated for the non-zero basis
+functions. Subsequent calls to GSREFIT reference this array to avoid
+recalculating basis functions at every call.
+.ih
+SEE ALSO
+gsfit, gsaccum, gsacpts, gssolve
+.endhelp
diff --git a/math/gsurfit/doc/gsreject.hlp b/math/gsurfit/doc/gsreject.hlp
new file mode 100644
index 00000000..96344ac2
--- /dev/null
+++ b/math/gsurfit/doc/gsreject.hlp
@@ -0,0 +1,44 @@
+.help gsreject Aug85 "Gsurfit Package"
+.ih
+NAME
+gsreject -- reject a data point from the fit
+.ih
+SYNOPSIS
+gsreject (sf, x, y, z, weight)
+
+.nf
+pointer	sf		# surface descriptor
+real	x		# x value, xmin <= x <= xmax
+real	y		# y value, ymin <= y <= ymax
+real	z		# data value
+real	weight		# weight
+.fi
+.ih
+ARGUMENTS
+.ls sf     
+Pointer to the surface descriptor structure.
+.le
+.ls x, y
+X and y values.
+.le
+.ls z      
+Data value.
+.le
+.ls weight
+Weight value.
+.le
+.ih
+DESCRIPTION
+GSREJECT removes a data point from the fit. The non-zero basis functions for
+each x and y are calculated, and the contribution of the point to the normal
+equations is computed and subtracted from the appropriate arrays and vectors.
+An array of points can be removed from the fit by repeated calls to GSREJECT
+followed by a single call to GSSOLVE to calculate a new set of coefficients.
+.ih
+NOTES
+Out of bounds x and y values and INDEF valued data values are the responsibility
+of the calling program.
+.ih
+SEE ALSO
+gsaccum, gsacpts
+.endhelp
diff --git a/math/gsurfit/doc/gsrestore.hlp b/math/gsurfit/doc/gsrestore.hlp
new file mode 100644
index 00000000..f71aab56
--- /dev/null
+++ b/math/gsurfit/doc/gsrestore.hlp
@@ -0,0 +1,36 @@
+.help gsrestore Aug85 "Gsurfit Package"
+.ih
+NAME
+gsrestore -- restore fit parameters
+.ih
+SYNOPSIS
+gsrestore (sf, fit)
+
+.nf
+pointer	sf		# surface descriptor
+real	fit[ARB]	# fit array
+.fi
+.ih
+ARGUMENTS
+.ls sf   
+Pointer to the surface descriptor structure. Returned by GSRESTORE.
+.le
+.ls fit    
+Array containing the surface parameters. The size of the fit array
+can be determined by a call gsstati.
+
+.nf
+	    len_fit = gsstati (gs, GSNSAVE)
+.fi
+.le
+.ih
+DESCRIPTION
+GSRESTORE returns the surface descriptor to the calling program and
+restores the surface parameters and fit ready for use by GSEVAL or
+GSVECTOR.
+.ih
+NOTES
+.ih
+SEE ALSO
+gssave
+.endhelp
diff --git a/math/gsurfit/doc/gssave.hlp b/math/gsurfit/doc/gssave.hlp
new file mode 100644
index 00000000..a6faf568
--- /dev/null
+++ b/math/gsurfit/doc/gssave.hlp
@@ -0,0 +1,39 @@
+.help gssave Aug85 "Gsurfit Package"
+.ih
+NAME
+gssave -- save parameters of the fit
+.ih
+SYNOPSIS
+call gssave (sf, fit)
+
+.nf
+pointer	sf		# surface descriptor
+real	fit[ARB]	# fit array
+.fi
+.ih
+ARGUMENTS
+.ls sf     
+Pointer to the surface descriptor structure.
+.le
+.ls fit    
+Array containing fit parameters. The size of the fit array can be determined
+by a call to gsstati.
+.le
+
+.nf
+	    len_fit = gsstati (sf, GSNSAVE)
+.fi
+.ih
+DESCRIPTION
+GSSAVE saves the surface parameters in the real array fit. The first eight
+elements of fit contain the surface_type, xorder, yorder, xterms, xmin,
+xmax, ymin and ymax. The coefficients are stored in the remaining array
+elements.
+.ih
+NOTES
+GSSAVE does not preserve the matrices and vectors used by the fitting
+routines.
+.ih
+SEE ALSO
+gsrestore
+.endhelp
diff --git a/math/gsurfit/doc/gsscoeff.hlp b/math/gsurfit/doc/gsscoeff.hlp
new file mode 100644
index 00000000..5f5e4a6a
--- /dev/null
+++ b/math/gsurfit/doc/gsscoeff.hlp
@@ -0,0 +1,35 @@
+.help gsscoeff Aug85 "Gsurfit Package"
+.ih
+NAME
+gsscoeff -- Procedure to set a coefficient
+.ih
+SYNOPSIS
+gsscoeff (sf, xorder, yorder, coeff)
+
+.nf
+    pointer	sf	# surface descriptor	
+    int		xorder	# x order of desired coefficient
+    int		yorder	# y order of desired coefficient
+    real	coeff	# coefficient value
+.fi
+.ih
+ARGUMENTS
+.ls sf    
+Pointer to the surface descriptor.
+.le
+.ls xorder, yorder
+The x and y order of the desired coefficient.
+.le
+.ls coeff
+The value of the coefficient to be set.
+.le
+.ih
+DESCRIPTION
+GSSCOEFF sets the coefficient of x ** (xorder - 1) * y ** (yorder - 1).
+GSSCOEFF returns if xorder and yorder are out of range.
+.ih
+NOTES
+.ih
+SEE ALSO
+gsgcoeff
+.endhelp
diff --git a/math/gsurfit/doc/gssolve.hlp b/math/gsurfit/doc/gssolve.hlp
new file mode 100644
index 00000000..8ddf42fc
--- /dev/null
+++ b/math/gsurfit/doc/gssolve.hlp
@@ -0,0 +1,40 @@
+.help gssolve Aug85 "Gsurfit Package"
+.ih
+NAME
+gssolve -- solve a linear system of equations by the Cholesky method
+.ih
+SYNOPSIS
+include <math/gsurfit.h>
+
+gssolve (sf, ier)
+
+.nf
+pointer	sf		# surface descriptor
+int	ier		# error code
+.fi
+.ih
+ARGUMENTS
+.ls sf     
+Pointer to the surface descriptor structure.
+.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, GSSOLVE
+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. GSSOLVE returns without
+fitting the data.
+.le
+.ih
+DESCRIPTION
+GSSOLVE computes the Cholesky factorization of the data matrix and
+solves for the coefficients
+of the fitting function by forward and back substitution. An error code is
+returned by GSSOLVE if it is unable to solve the normal equations as
+formulated.
+.ih
+NOTES
+.ih
+SEE ALSO
+gsfit, gsrefit, gsaccum, gsacpts
+.endhelp
diff --git a/math/gsurfit/doc/gsstati.hlp b/math/gsurfit/doc/gsstati.hlp
new file mode 100644
index 00000000..7b1b7f2e
--- /dev/null
+++ b/math/gsurfit/doc/gsstati.hlp
@@ -0,0 +1,35 @@
+.help gsstati Aug85 "Gsurfit Package"
+.ih
+NAME
+include <math/gsurfit.h>
+
+gsstati -- get integer parameter
+.ih
+SYNOPSIS
+ival = gsstati (sf, parameter)
+
+.nf
+pointer	sf		# surface descriptor
+int	parameter	# integer parameter to be returned
+.fi
+.ih
+ARGUMENTS
+.ls sf     
+The pointer to the surface descriptor structure.
+.le
+.ls parameter
+Parameter to be returned. The options are GSTYPE, GSXORDER, GSYORDER,
+GSNXCOEFF, GSNYCOEFF, and GSNSAVE. The parameter definitions are
+found in the package header file math/gsurfit.h.
+.le
+.ih
+DESCRIPTION
+The values of the integer parameters are returned. The parameters include
+the surface_type, the x and y orders, the number of x and y coefficients
+and the length of the buffer required by GSSAVE.
+.ih
+NOTES
+.ih
+SEE ALSO
+gsstatr
+.endhelp
diff --git a/math/gsurfit/doc/gsstatr.hlp b/math/gsurfit/doc/gsstatr.hlp
new file mode 100644
index 00000000..5a5578f2
--- /dev/null
+++ b/math/gsurfit/doc/gsstatr.hlp
@@ -0,0 +1,34 @@
+.help gsstatr Aug85 "Gsurfit Package"
+.ih
+NAME
+gsstatr -- get real parameter
+.ih
+SYNOPSIS
+include <math/gsurfit.h>
+
+rval = gsstatr (sf, parameter)
+
+.nf
+pointer	sf		# surface descriptor
+real	parameter	# real parameter to be returned
+.fi
+.ih
+ARGUMENTS
+.ls sf     
+The pointer to the surface descriptor structure.
+.le
+.ls parameter
+Parameter to be returned. The options are GSXMIN, GSXMAX, GSYMIN and
+and GSYMAX. The parameter definitions are
+found in the package header file math/gsurfit.h.
+.le
+.ih
+DESCRIPTION
+The values of the integer parameters are returned. The parameters include
+the minimum and maximum x values and the minimum and maximum y values.
+.ih
+NOTES
+.ih
+SEE ALSO
+gsstati
+.endhelp
diff --git a/math/gsurfit/doc/gssub.hlp b/math/gsurfit/doc/gssub.hlp
new file mode 100644
index 00000000..2c771612
--- /dev/null
+++ b/math/gsurfit/doc/gssub.hlp
@@ -0,0 +1,35 @@
+.help gssub Aug85 "Gsurfit Package"
+.ih
+NAME
+gssub -- subtract surface 1 from surface 2
+.ih
+SYNOPSIS
+gssub (sf1, sf2, sf3)
+
+.nf
+pointer	sf1		# first surface descriptor
+pointer	sf2		# second surface descriptor
+pointer	sf3		# resultant surface descriptor
+.fi
+.ih
+ARGUMENTS
+.ls sf1    
+Pointer to the first surface descriptor.
+.le
+.ls sf2   
+Pointer to the second surface descriptor.
+.le
+.ls sf3    
+Pointer to the resultant surface descriptor.
+.le
+.ih
+DESCRIPTION
+The coefficients of surface 2 are subtracted from surface 1. GSSUB checks
+that the surface_types are the same and that the fits are normalized over
+the same range of data.
+.ih
+NOTES
+.ih
+SEE ALSO
+gscopy, gsadd
+.endhelp
diff --git a/math/gsurfit/doc/gsurfit.hd b/math/gsurfit/doc/gsurfit.hd
new file mode 100644
index 00000000..169961f7
--- /dev/null
+++ b/math/gsurfit/doc/gsurfit.hd
@@ -0,0 +1,25 @@
+# Help directory for the GSURFIT (surface fitting) package.
+
+$gsurfit	= "math$gsurfit/"
+
+gsaccum		hlp = gsaccum.hlp,	src = gsurfit$gsaccum.x
+gsacpts		hlp = gsacpts.hlp,	src = gsurfit$gsacpts.x
+gsadd		hlp = gsadd.hlp,	src = gsurfit$gsadd.x
+gscoeff		hlp = gscoeff.hlp,	src = gsurfit$gscoeff.x
+gscopy		hlp = gscopy.hlp,	src = gsurfit$gscopy.x
+gsder		hlp = gsder.hlp,	src = gsurfit$gsder.x
+gserrors	hlp = gserrors.hlp,	src = gsurfit$gserrors.x
+gseval		hlp = gseval.hlp,	src = gsurfit$gseval.x
+gsinit		hlp = gsinit.hlp,	src = gsurfit$gsinit.x
+gsfit		hlp = gsfit.hlp,	src = gsurfit$gsfit.x
+gsfree		hlp = gsfree.hlp,	src = gsurfit$gsfree.x
+gsrefit		hlp = gsrefit.hlp,	src = gsurfit$gsrefit.x
+gsreject	hlp = gsreject.hlp,	src = gsurfit$gsreject.x
+gsrestore	hlp = gsrestore.hlp,	src = gsurfit$gsrestore.x
+gssave		hlp = gssave.hlp,	src = gsurfit$gssave.x
+gssolve		hlp = gssolve.hlp,	src = gsurfit$gssolve.x
+gsstati		hlp = gsstati.hlp,	src = gsurfit$gsstat.x
+gsstatr		hlp = gsstatr.hlp,	src = gsurfit$gsstat.x
+gssub		hlp = gssub.hlp,	src = gsurfit$gssub.x
+gsvector	hlp = gsvector.hlp,	src = gsurfit$gsvector.x
+gszero		hlp = gszero.hlp,	src = gsurfit$gszero.x
diff --git a/math/gsurfit/doc/gsurfit.hlp b/math/gsurfit/doc/gsurfit.hlp
new file mode 100644
index 00000000..99f42444
--- /dev/null
+++ b/math/gsurfit/doc/gsurfit.hlp
@@ -0,0 +1,169 @@
+.help gsurfit Aug85 "Math Package"
+.ih
+NAME
+gsurfit -- surface fitting package
+.ih
+SYNOPSIS
+
+.nf
+        gsinit (sf, surf_type, xorder, yorder, xterms, xmin, xmax, ymin, ymax)
+       gsaccum (sf, x, y, z, w, wtflag)
+       gsacpts (sf, x, y, z, w, npts, wtflag)
+      gsreject (sf, x, y, z, w)
+       gssolve (sf, ier)
+         gsfit (sf, x, y, z, w, npts, wtflag, ier)
+       gsrefit (sf, x, y, z, w, ier)
+    y = gseval (sf, x, y)
+      gsvector (sf, x, y, zfit, npts)
+         gsder (sf, x, y, zfit, npts, nxder, nyder)
+       gscoeff (sf, coeff, ncoeff)
+      gserrors (sf, z, w, zfit, rms, errors)
+        gssave (sf, fit)
+     gsrestore (sf, fit)
+ival = gsstati (sf, param)
+rval = gsstatr (sf, param)
+         gsadd (sf1, sf2, sf3)
+         gssub (sf1, sf2, sf3)
+        gscopy (sf1, sf2)
+        gsfree (sf)
+.fi
+.ih
+DESCRIPTION
+The gsurfit package provides a set of routines for fitting data to functions
+of two variables, linear in their coefficients, using least squares
+techniques. The numerical technique employed is the solution of the normal
+equations by the Cholesky method.
+.ih
+NOTES
+The fitting function is chosen at run time from the following list.
+
+.nf
+	GS_LEGENDRE		# Lengendre polynomials in x and y
+	GS_CHEBYSHEV		# Chebyshev polynomials in x and y
+.fi
+
+The gsurfit package performs a weighted fit. The weighting options are
+WTS_USER and WTS_UNIFORM. The user must supply a weight array. In the
+WTS_UNIFORM mode the gsurfit routines set the weights to 1. In WTS_USER
+mode the user must supply an array of weight values.
+
+In WTS_UNIFORM mode the reduced chi-square returned by GSERRORS is the
+variance of the fit and the errors in the coefficients are scaled
+by the square root of this variance. Otherwise the weights are
+interpreted as one over the variance of the data and the true reduced
+chi-square is returned.
+
+The routines assume that all the x and y values of interest lie in the
+region xmin <= x <= xmax and ymin <= y <= ymax. Checking for out of
+bounds x and y values is the responsibility of the calling program.
+The package routines assume that INDEF valued data points have been removed
+from the data set prior to entering the package routines.
+
+In order to make the package definitions available to the calling program
+an include <math/gsurfit.h> statement must be inserted in the calling program.
+GSINIT must be called before each fit. GSFREE frees the space used by
+the GSURFIT package.
+.ih
+EXAMPLES
+.nf
+Example 1: Fit surface to data, uniform weighting
+
+include <math/gsurfit.h>
+
+...
+
+call gsinit (sf, GS_CHEBYSHEV, 4, 4, NO, 1., 512., 1., 512.)
+
+call gsfit (sf, x, y, z, w, npts, WTS_UNIFORM, ier)
+if (ier != OK)
+    call error (...)
+
+do i = 1, 512 {
+    call printf ("x = %g y = %g z = %g zfit = %g\n")
+	call pargr (x[i])
+	call pargr (y[i])
+	call pargr (z[i])
+	call pargr (gseval (sf, x[i], y[i]))
+}
+
+call gsfree (sf)
+
+...
+
+Example 2: Fit a surface using accumulate mode, uniform weighting
+
+include <math/gsurfit.h>
+
+...
+
+do i = 1, 512 {
+    if (y[i] != INDEF)
+	call gsaccum (sf, x[i], y[i], z[i], weight[i], WTS_UNIFORM)
+}
+
+call gssolve (sf, ier)
+if (ier != OK)
+    call error (...)
+
+...
+
+call gsfree (sf)
+
+...
+
+Example 3: Fit and subtract a smooth surface from image lines
+
+include <math/gsurfit.h>
+
+...
+
+call gsinit (gs, GS_CHEBYSHEV, xorder, yorder, YES, 1., 512., 1., 512.)
+
+call gsfit (sf, xpts, ypts, zpts, w, WTS_UNIFORM, ier)
+if (ier != OK)
+    call error (...)
+
+do i = 1, 512
+    Memr[x+i-1] = i
+
+do line = 1, 512 {
+
+    inpix = imgl2r (im, line)
+    outpix = imgl2r (im, line)
+
+    yval = line
+    call amovkr (yval, Memr[y], 512)
+    call gsvector (sf, Memr[x], Memr[y], Memr[outpix], 512)
+    call asubr (Memr[inpix], Memr[outpix], Memr[outpix, 512) 
+}
+
+call gsfree (sf)
+
+...
+
+Example 4: Fit curve, save fit for later use for GSEVAL
+
+include <math/gsurfit.h>
+
+call gsinit (sf, GS_LEGENDRE, xorder, yorder, YES, xmin, xmax, ymin, ymax)
+
+call gsfit (sf, x, y, z, w, npts, WTS_UNIFORM, ier)
+if (ier != OK)
+    ...
+
+nsave = gsstati (sf, GSNSAVE)
+call salloc (fit, nsave, TY_REAL)
+call gssave (sf, Memr[fit])
+call gsfree (sf)
+
+...
+
+call gsrestore (sf, Memr[fit])
+do i = 1, npts
+    zfit[i] = gseval (sf, x[i], y[i])
+
+call gsfree (sf)
+
+...
+.fi
+.endhelp
diff --git a/math/gsurfit/doc/gsurfit.men b/math/gsurfit/doc/gsurfit.men
new file mode 100644
index 00000000..0b938ac8
--- /dev/null
+++ b/math/gsurfit/doc/gsurfit.men
@@ -0,0 +1,21 @@
+	gsaccum	- Accumulate point into data set
+	gsacpts - Accumulate points into a data set
+	  gsadd - Add two surfaces
+	gscoeff - Get coefficients 
+	 gscopy - Copy one surface to another
+	  gsder - Evaluate the derivatives of a surface
+       gserrors	- Calculate chi-square and errors in coefficients
+	 gseval	- Evaluate surface at x and y
+	  gsfit	- Fit surface
+	 gsfree - Free space allocated by gsinit
+	 gsinit - Make ready to fit a surface; set up parameters of fit
+	gsrefit	- Refit surface, same x, y and weight, different z
+       gsreject	- Reject point from data set
+      gsrestore - Restore surface parameters and coefficients
+	 gssave - Save surface parameters and coefficients
+	gssolve	- Solve matrix for coefficients
+	gsstati - Get integer parameter
+	gsstatr - Get real parameter
+	  gssub - Subtract one surface from another
+       gsvector - Evaluate surface at an array of x and y
+	 gszero - Zero arrays for new fit
diff --git a/math/gsurfit/doc/gsvector.hlp b/math/gsurfit/doc/gsvector.hlp
new file mode 100644
index 00000000..16003101
--- /dev/null
+++ b/math/gsurfit/doc/gsvector.hlp
@@ -0,0 +1,41 @@
+.help gsvector Aug85 "Gsurfit Package"
+.ih
+NAME
+gsvector -- evaluate the fitted surface at a set of points
+.ih
+SYNOPSIS
+gsvector (sf, x, y, zfit, npts)
+
+.nf
+pointer	sf		# surface descriptor
+real	x[npts]		# x array, xmin <= x <= xmax
+real	y[npts]		# y array, ymin <= y <= ymax
+real	zfit[npts]	# data values
+int	npts		# number of data points
+.fi
+.ih
+ARGUMENTS
+.ls sf    
+Pointer to the surface descriptor structure.
+.le
+.ls x, y
+Array 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 data points. GSVECTOR 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
+gseval, gsder
+.endhelp
diff --git a/math/gsurfit/doc/gszero.hlp b/math/gsurfit/doc/gszero.hlp
new file mode 100644
index 00000000..c7d411dd
--- /dev/null
+++ b/math/gsurfit/doc/gszero.hlp
@@ -0,0 +1,27 @@
+.help gszero Aug85 "Gsurfit Package"
+.ih
+NAME
+gszero -- set up for a new surface fit
+.ih
+SYNOPSIS
+gszero (sf)
+
+.nf
+pointer	sf		# surface descriptor
+.fi
+.ih
+ARGUMENTS
+.ls sf     
+Pointer to the surface descriptor structure.
+.le
+.ih
+DESCRIPTION
+GSZERO zeros the matrix and right side of the matrix equation.
+.ih
+NOTES
+GSZERO can be used to reinitialize the matrix and right side of the
+equation to begin a new fit in accumulate mode.
+.ih
+SEE ALSO
+gsinit, gsfit, gsrefit, gsaccum, gsacpts
+.endhelp