Initial commit

61 files changed, 3988 insertions, 0 deletions

diff --git a/math/nlfit/README b/math/nlfit/README
new file mode 100644
index 00000000..c47732af
--- /dev/null
+++ b/math/nlfit/README
@@ -0,0 +1,81 @@
+                      THE NLFIT PACKAGE
+
+This subdirectory contains the routines in the non-linear least squares
+fitting package NLFIT. NLFIT uses the Levenberg-Marquardt method  to
+to solve for the parameters of a user specified non-linear equation.
+The user must supply two routines. The first routine evaluates the function
+in terms of its parameters. The second routine evaluates the function and its
+derivatives in terms of its parameters.  The user must also supply initial
+guesses for the parameters and parameter increments, the list of
+parameters to be varied during the fitting processa, a fitting
+tolerance and the maximum number of iterations.
+
+The principal entry points into the NLFIT are listed below.
+
+	nlinit - Initialize the fitting routines
+       nlstati - Get the value of an integer NLFIT parameter
+        nlstat - Get the value of floating point NLFIT parameter
+	nlpget - Get the values of the fitted parameters
+	nlfree - Free memory allocated by nlinit
+	 nlfit - Fit the function
+	nleval - Evaluate the curve at a given point
+      nlvector - Evaluate the curve at an array of points 
+      nlerrors - Compute the fits statistics and errors in the parameters
+
+The calling sequences for the above routines are listed below. The [ird]
+stand for integer, real and double precision versions of each routine
+respectively.
+
+               nlinit[rd] (nl, address(fnc), address(dfnc), params, dparams,
+	            nparams, plist, nfparams, tol, itmax)
+        ival = nlstati (nl, param)
+     [rd]val = nlstat[rd] (nl, param)
+               nlpget[rd] (nl, params, nparams)
+               nlfree[rd] (nl)
+               nlfit[rd] (nl, x, z, w, npts, nvars, wtflag, ier)
+     [rd]val = nleval[rd] (nl, x, nvars)
+               nlvector[rd] (nl, x, zfit, npts, nvars)
+               nlerrors[rd] (nl, z, zfit, wts, npts, variance, chisqr, errors)
+
+The user supplied functions fnc and dfnc have the following calling
+sequences.
+	
+	       fnc (x, nvars, params, nparams, zfit)
+	      dfnc (x, nvars, params, dparams, nparams, zfit, derivs)
+
+    The address of the user supplied function can be determined with a
+call to locpr as in
+
+	      address = locpr (fnc)
+
+The user definitions for the NLFIT package can be found in the file 
+lib$math/nlfit.h and can be made available to user applications programs
+with the statement "include <math/nlfit.h>".
+
+The permitted values for the param argument in nlstat[ird] are the following.
+
+define	NLNPARAMS	integer		# Number of parameters
+define	NLNFPARAMS	integer         # Number of fitted parameters
+define	NLITMAX		integer		# Maximum number of iterations
+define	NLITER		integer		# Current number of iterations
+define	NLNPTS		integer         # Number of points
+define	NLSUMSQ		real/double	# Current reduced chi-squared
+define	NLOLDSQ		real/double     # Previous reduced chi-squared
+define	NLLAMBDA	real/double     # Value of lambda factor
+define	NLTOL		real/double	# Fitting tolerance in %chi-squared
+define	NLSCATTER	real/double     # Mean scatter in the fit
+
+The permitted values of the wtflag argument in nlfit[rd] are the following.
+
+define	WTS_USER			# User enters weights
+define	WTS_UNIFORM			# Equal weights
+define	WTS_CHISQ 			# Chi-squared weights
+define	WTS_SCATTER			# Weights include scatter term
+
+The permitted error values returned from nlfit[rd] are the following.
+
+define	DONE			0	# Solution converged
+define	SINGULAR		1	# Singular matrix
+define	NO_DEG_FREEDOM		2	# Too few points
+define	NOT_DONE		3	# Solution did not converge
+
diff --git a/math/nlfit/doc/nlerrors.hlp b/math/nlfit/doc/nlerrors.hlp
new file mode 100644
index 00000000..3b463c14
--- /dev/null
+++ b/math/nlfit/doc/nlerrors.hlp
@@ -0,0 +1,67 @@
+.help nlerrors Feb91 "Nlfit Package"
+.ih
+NAME
+nlerrors -- compute the fit statistics and errors in the parameters
+.ih
+SYNOPSIS
+nlerrors[rd] (nl, z, zfit, w, npts, variance, chisqr, errors)
+
+.nf
+pointer		nl		# curve descriptor
+real/double	z[npts]		# array of input function values
+real/double	zfit[npts]	# array of fitted function values
+real/double	w[npts]		# array of weights
+int		npts		# number of data points
+real/double	variance	# the computed variance of the fit
+real/double	chisqr		# the computed reduced chi-square of the fit
+real/double	errors[*]	# errors in the fitted parameters
+.fi
+.ih
+ARGUMENTS
+.ls nl     
+Pointer to the curve descriptor structure.
+.le
+.ls z   
+Array of function values.
+.le
+.ls zfit  
+Array of fitted function values.
+.le
+.ls w
+Array of weights.
+.le
+.ls npts
+The number of data points.
+.le
+.ls variance
+The computed variance of the fit.
+.le
+.ls chisqr
+The computed reduced chi-squared of the fit.
+.le
+.ls errors
+Array of errors in the computed parameters.
+.le
+.ih
+DESCRIPTION
+Compute the variance and reduced chi-squared of the fit and the
+errors in the fitted parameters.
+.ih
+NOTES
+The reduced chi-squared 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.
+The variance of the fit is the square root of the sum of the
+squares of the residuals divided by the number of degrees of freedom.
+If the weighting is uniform, then the reduced chi-squared is equal to the
+variance of the fit.
+The error of the j-th parameter is the square root of the j-th diagonal
+element of the inverse of the data matrix. If the weighting is uniform,
+then the errors are scaled by the square root of the variance of the data.
+
+The zfit array can be computed by a call to the nlvector[rd] routine.
+The size of the array required to hold the output error array can be
+determined by a call to nlstati.
+.ih
+SEE ALSO
+nlvector,nlstat
+.endhelp
diff --git a/math/nlfit/doc/nleval.hlp b/math/nlfit/doc/nleval.hlp
new file mode 100644
index 00000000..f0bbb400
--- /dev/null
+++ b/math/nlfit/doc/nleval.hlp
@@ -0,0 +1,35 @@
+.help nleval Feb91 "Nlfit Package"
+.ih
+NAME
+nleval -- evaluate the fitted function at a single point
+.ih
+SYNOPSIS
+z = nleval[rd] (nl, x, nvars)
+
+.nf
+pointer		nl		# curve descriptor
+real/double	x[nvars]	# array of variables
+int		nvars		# the number of variables
+.fi
+.ih
+ARGUMENTS
+.ls nl    
+Pointer to the curve descriptor structure.
+.le
+.ls x   
+Array of variable values at which the curve is to be evaluated.
+.le
+.ls nvars
+The number of variables.
+.le
+.ih
+DESCRIPTION
+Evaluate the curve at the specified point. NLEVAL is a real or double
+function which returns the fitted z value.
+.ih
+NOTES
+NLEVAL uses the parameter array stored in the curve descriptor structure.
+.ih
+SEE ALSO
+nlvector
+.endhelp
diff --git a/math/nlfit/doc/nlfit.hd b/math/nlfit/doc/nlfit.hd
new file mode 100644
index 00000000..13bb3065
--- /dev/null
+++ b/math/nlfit/doc/nlfit.hd
@@ -0,0 +1,12 @@
+# Help directory for the NLFIT (non-linear least-squares fitting) package.
+
+$nlfit	= "math$nlfit/"
+
+nlerrors	hlp = nlerrors.hlp,	src = nlfit$nlerrorsr.x
+nleval		hlp = nleval.hlp,	src = nlfit$nlevalr.x
+nlinit		hlp = nlinit.hlp,	src = nlfit$nlinitr.x
+nlfit		hlp = nlfit.hlp,	src = nlfit$nlfitr.x
+nlfree		hlp = nlfree.hlp,	src = nlfit$nlfreer.x
+nlpget		hlp = nlpget.hlp,	src = nlfit$nlpgetr.x
+nlstat		hlp = nlstat.hlp,	src = nlfit$nlstatr.x
+nlvector	hlp = nlvector.hlp,	src = nlfit$nlvectorr.x
diff --git a/math/nlfit/doc/nlfit.hlp b/math/nlfit/doc/nlfit.hlp
new file mode 100644
index 00000000..8ed82512
--- /dev/null
+++ b/math/nlfit/doc/nlfit.hlp
@@ -0,0 +1,81 @@
+.help nlfit Feb91 "Nlfit Package"
+.ih
+NAME
+nlfit -- fit a curve to a set of data values
+.ih
+SYNOPSIS
+nlfit (nl, x, z, w, npts, nvars, wtflag, ier)
+
+.nf
+pointer		nl		# curve descriptor
+real/double	x[npts*nvars]	# variable values (stored as x[nvars,npts])
+real/double	z[npts]		# array of z values
+real/double	w[npts]		# array of weights
+int		npts		# number of data points
+int		wtflag		# type of weighting
+int		ier		# error code
+.fi
+.ih
+ARGUMENTS
+.ls nl    
+Pointer to the curve descriptor structure.
+.le
+.ls x    
+Array of x values.
+.le
+.ls z
+Array of function values.
+.le
+.ls w
+Array of weights.
+.le
+.ls npts
+The number of data points.
+.le
+.ls wtflag
+Type of weighting. The options are WTS_USER, WTS_UNIFORM, WTS_CHISQ and
+WTS_SCATTER. 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. If wtflag = WTS_CHISQ the weights are set equal to the
+reciprocal of the function value. If wtflag = WTS_SCATTER the fitting
+routine adds a scatter term to the weights if the reduced chi-squared of
+the fit is significantly greater than 1.
+.le
+.ls ier
+Error code for the fit. The options are DONE, SINGULAR,
+NO_DEG_FREEDOM and NOT_DONE. 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.
+If ier = NOT_DONE, the fit could not achieve the desired tolerance with
+the given input data in the specified maximum number of iterations.
+.le
+.ih
+DESCRIPTION
+NLFIT accumulate the data into the appropriate internal matrices and vectors
+and does the fit.
+.ih
+NOTES
+The permitted values of the input wtflag argument to nlfit[rd] are the
+following.
+
+.nf
+define	WTS_USER			# User enters weights
+define	WTS_UNIFORM			# Equal weights
+define	WTS_CHISQ 			# Chi-squared weights
+define	WTS_SCATTER			# Weights include scatter term
+.fi
+
+The permitted error values returned from nlfit[rd] are the following.
+
+.nf
+define	DONE			0	# Solution converged
+define	SINGULAR		1	# Singular matrix
+define	NO_DEG_FREEDOM		2	# Too few points
+define	NOT_DONE		3	# Solution did not converge
+.fi
+.ih
+SEE ALSO
+.endhelp
diff --git a/math/nlfit/doc/nlfit.men b/math/nlfit/doc/nlfit.men
new file mode 100644
index 00000000..5fb68623
--- /dev/null
+++ b/math/nlfit/doc/nlfit.men
@@ -0,0 +1,8 @@
+      nlerrors - Compute the fits statistics and errors in the parameters
+	 nlfit - Fit the function
+	nlfree - Free memory allocated by nlinit
+	nlinit - Initialize the fitting routines
+	nlpget - Get the values of the fitted parameters
+        nlstat - Get the value of an NLFIT parameter
+	nleval - Evaluate the curve at a given point
+      nlvector - Evaluate the curve at an array of points 
diff --git a/math/nlfit/doc/nlfree.hlp b/math/nlfit/doc/nlfree.hlp
new file mode 100644
index 00000000..9688d5e7
--- /dev/null
+++ b/math/nlfit/doc/nlfree.hlp
@@ -0,0 +1,26 @@
+.help nlfree Feb91 "Nlfit Package"
+.ih
+NAME
+nlfree -- free the curve descriptor structure
+.ih
+SYNOPSIS
+nlfree[rd] (nl)
+
+.nf
+pointer	nl	# curve descriptor
+.fi
+.ih
+ARGUMENTS
+.ls nl     
+Pointer to the curve descriptor structure.
+.le
+.ih
+DESCRIPTION
+Frees the curve descriptor structure.
+.ih
+NOTES
+NLFREE should be called after each curve fit.
+.ih
+SEE ALSO
+nlinit
+.endhelp
diff --git a/math/nlfit/doc/nlinit.hlp b/math/nlfit/doc/nlinit.hlp
new file mode 100644
index 00000000..94e5e5a9
--- /dev/null
+++ b/math/nlfit/doc/nlinit.hlp
@@ -0,0 +1,86 @@
+.help nlinit Feb91 "Nlfit Package"
+.ih
+NAME
+nlinit -- initialise curve descriptor
+.ih
+SYNOPSIS
+include <math/nlfit.h>
+
+nlinit[rd] (nl, fnc, dfnc, params, dparams, nparams, plist, nfparams,
+                tol, itmax)
+
+.nf
+pointer		nl		# curve descriptor
+int		fnc		# address of first user-supplied function
+int		dfnc		# address of second user-supplied function
+real/double	params[nparams]	# list of initial parameter values
+real/double	dparams[nparams]# list of parameter increments
+int		nparams		# number of parameters
+int		plist[nparams]	# list of parameters to be fit
+int		nfparams	# number of parameters to be fit
+real/double	tol		# fitting tolerance
+int		itmax		# maximum number of iterations
+.fi
+.ih
+ARGUMENTS
+.ls nl             
+Pointer to the curve descriptor structure.
+.le
+.ls fnc       
+The address of the user-supplied subroutine fncname for evaluating the function
+to be fit. The calling sequence of fncname is the following.
+
+.nf
+fncname (x, nvars, params, nparams, zfit)
+.fi
+.le
+.ls dfnc       
+The address of the user-supplied subroutine dfncname for evaluating the
+function to be fit and its derivatives with respect to the parameters.
+The calling sequence of dfncname is the following.
+
+.nf
+dfncname (x, nvars, params, dparams, nparams, zfit, derivs)
+.fi
+.le
+.ls params
+An array containing initial values of all the parameters, including both
+those to be fit and those that are to be held constant.
+.le
+.ls dparams
+An array parameter increment values over which the derivatives are to be
+computed empirically. If equations are supplied for the derivatives
+inside dfnc, the dparams array is not used.
+.le
+.ls nparams
+The total number of parameters in the user-supplied function.
+.le
+.ls plist
+The list of parameters to be fit. Plist is an integer array containing
+the indices of those parameters in params to be fit.
+.le
+.ls nfparams
+The total number of parameters to be fit.
+.le
+.ls tol   
+The fitting tolerance. If the difference in the chi-squared
+from one iteration to the next is less than the fitting tolerance
+times the current chi-squared, the fit has converged.
+.le
+.ls itmax
+The maximum number of fitting iterations. Itmax must be greater than
+or equal to 3.
+.le
+.ih
+DESCRIPTION
+Allocate space for the curve descriptor structure and the arrays and
+vectors used by the numerical routines. Initialize all arrays and vectors
+to zero. Return the curve descriptor to the calling routine.
+.ih
+NOTES
+NLINIT must be the first NLFIT routine called. NLINIT returns a NULL pointer
+if it encounters an illegal parameter list.
+.ih
+SEE ALSO
+nlfree
+.endhelp
diff --git a/math/nlfit/doc/nllmfit.hlp b/math/nlfit/doc/nllmfit.hlp
new file mode 100644
index 00000000..59309133
--- /dev/null
+++ b/math/nlfit/doc/nllmfit.hlp
@@ -0,0 +1,172 @@
+.help nlfit Feb91 "Math Package"
+.ih
+NAME
+nlfit -- Levenberg-Marquardt non-linear least-squares fitting package
+.ih
+SYNOPSIS
+The principal entry points into the NLFIT package are listed below.
+
+.nf
+	nlinit - Initialize the fitting routines
+        nlstat - Get the value of an NLFIT parameter
+	nlpget - Get the values of the fitted parameters
+	nlfree - Free memory allocated by nlinit
+	 nlfit - Fit the function
+	nleval - Evaluate the curve at a given point
+      nlvector - Evaluate the curve at an array of points 
+      nlerrors - Compute the fits statistics and errors in the parameters
+.fi
+
+The calling sequences for the above routines are listed below. The [ird]
+stand for integer, real and double precision versions of each routine
+respectively.
+
+.nf
+            nlinit[rd] (nl, address(fnc), address(dfnc), params, dparams,
+	        nparams, plist, nfparams, tol, itmax)
+        ival = nlstati (nl, param)
+  [rd]val = nlstat[rd] (nl, param)
+            nlpget[rd] (nl, params, nparams)
+            nlfree[rd] (nl)
+             nlfit[rd] (nl, x, z, w, npts, nvars, wtflag, ier)
+  [rd]val = nleval[rd] (nl, x, nvars)
+          nlvector[rd] (nl, x, zfit, npts, nvars)
+          nlerrors[rd] (nl, z, zfit, wts, npts, variance, chisqr, errors)
+.fi
+
+The user supplied functions fnc and dfnc must have the following form.
+	
+.nf
+	       fnc (x, nvars, params, nparams, zfit)
+	      dfnc (x, nvars, params, dparams, nparams, zfit, derivs)
+.fi
+
+The addresses of these functions can be obtained by a call to lcopr as
+follows.
+
+.nf
+		address = locpr (fnc)
+.fi
+
+.ih
+DESCRIPTION
+The NLFIT package provides a set of routines for fitting data to non-linear
+functions of several variables using least squares techniques.
+NLFIT uses the Levenberg-Marquardt
+method to solve for the parameters of a user specified non-linear equation.
+The user must supply two subroutines. The first subroutine evaluates the
+function in terms of its parameters. The second subroutine evaluates the
+function and its derivatives in terms of its parameters.  The user must
+also supply initial
+guesses for the parameters and parameter increments, the list of
+parameters to be varied during the fitting process, a fitting
+tolerance and the maximum number of iterations. 
+.ih
+NOTES
+The poackage definitions for NLFIT can be found in the file 
+lib$math/nlfit.h. In order to make these definitions available
+to the calling program, the user must insert the statement
+"include <math/nlfit.h>" inside the calling program.
+
+The permitted values for the param argument in nlstat[ird] are the following.
+
+.nf
+define	NLNPARAMS	integer		# Number of parameters
+define	NLNFPARAMS	integer         # Number of fitted parameters
+define	NLITMAX		integer		# Maximum number of iterations
+define	NLITER		integer		# Current number of iterations
+define	NLNPTS		integer         # Number of points
+define	NLSUMSQ		real/double	# Current reduced chi-squared
+define	NLOLDSQ		real/double     # Previous reduced chi-squared
+define	NLLAMBDA	real/double     # Value of lambda factor
+define	NLTOL		real/double	# Fitting tolerance in %chi-squared
+define	NLSCATTER	real/double     # Mean scatter in the fit
+.fi
+
+The permitted values of the wtflag argument in nlfit[rd] are the following.
+
+.nf
+define	WTS_USER			# User enters weights
+define	WTS_UNIFORM			# Equal weights
+define	WTS_CHISQ 			# Chi-squared weights
+define	WTS_SCATTER			# Weights include scatter term
+.fi
+
+The permitted error values returned from nlfit[rd] are the following.
+
+.nf
+define	DONE			0	# Solution converged
+define	SINGULAR		1	# Singular matrix
+define	NO_DEG_FREEDOM		2	# Too few points
+define	NOT_DONE		3	# Solution did not converge
+.fi
+.ih
+REFERENCES
+1. Bevington,P.R., 1969, Data Reduction and Error Analysis for the Physical
+Sciences, Chapter 11, page 235.
+
+2. Press, W.H. et al., 1986, Numerical Recipes: The Art of Scientific
+Computing, Chapter 14, page 523
+.ih
+EXAMPLES
+.nf
+Example 1: Fit a curve to the data using uniform weighting. Fit all the
+parameters.
+
+    # Include nlfit definitions
+    include <math/nlfit.h>
+
+    # Declare the variables
+    int	nparams, nfparams, itmax, npts, nvars, ier
+    int	plist[nparams]
+    real tol, variance, chisqr
+    real params[nparams], dparams[nparams], vars[nvars,npts], z[npts]
+    real w[npts], zfit[npts], errors[nparams]
+
+    # Declare the functions
+    extern fncname(), dfncname()
+    int	locpr()
+
+    begin
+            ... get data, set initial values of all parameters, parameter
+	    ... increments, tol and itmax
+
+            # Define list of variables to be fitted
+            nfparams = nparams
+            do i = 1, nparams
+	        plist[i] = i
+
+	    # Fit only parameters 1,3 and 5
+	    #nfparams = 3
+	    #plist[1] = 1
+	    #plist[2] = 3
+	    #plist[3] = 5
+
+            # Initialize the fitting routines
+            call nlinitr (nl, locpr (fncname), locpr (dfncname),
+	        params, dparams, nparams, plist, nfparams, tol, itmax)
+
+            # Fit the data
+            call nlfitr (nl, x, z, w, npts, nvars, WTS_UNIFORM, ier)
+            if (ier != DONE)
+	        ... take error action
+
+            # Compute the statistics of the fit
+            call nlvectorr (nl, x, zfit, npts, nvars)
+            call nlerrorsr (nl, z, zfit, wts, npts, variance, chisqr,
+	        errors)
+
+            # Get the values of the fitted parameters
+            call nlpgetr (nl, params, nparams)
+            call nlfreer (nl)
+
+            # Print the parameters and their errors.
+            do i = 1, nparams {
+	        call printf ("%d  %g  %g\n")
+	        call pargi (i)
+	        call pargr (params[i])
+	        call pargr (errors[i])
+	    }
+    end
+.fi
+.endhelp
diff --git a/math/nlfit/doc/nlpget.hlp b/math/nlfit/doc/nlpget.hlp
new file mode 100644
index 00000000..07f2ef5e
--- /dev/null
+++ b/math/nlfit/doc/nlpget.hlp
@@ -0,0 +1,38 @@
+.help nlpget Feb91 "Nlfit Package"
+
+.ih
+NAME
+nlpget -- get the number and values of the fitted parameters
+.ih
+SYNOPSIS
+nlpget[rd] (nl, params, nparams)
+
+.nf
+pointer		nl		# curve descriptor
+real/double	params[nparams]	# the parameters array
+int		nparams		# the number of parameters
+.fi
+.ih
+ARGUMENTS
+.ls nl
+Pointer to the curve descriptor.
+.le
+.ls params
+Array of fitted parameters.
+.le
+.ls nparameters
+The total number of parameters.
+.le
+.ih
+DESCRIPTION
+NLPGET fetches the parameters array and the number of parameters from the
+curve descriptor structure. All the parameters, both those that were
+fit and those that were held constant are output.
+.ih
+NOTES
+The array size required to hold the output parameters can be determined
+by a call to nlstati.
+.ih
+SEE ALSO
+nlstat
+.endhelp
diff --git a/math/nlfit/doc/nlstat.hlp b/math/nlfit/doc/nlstat.hlp
new file mode 100644
index 00000000..a745a663
--- /dev/null
+++ b/math/nlfit/doc/nlstat.hlp
@@ -0,0 +1,57 @@
+.help nlstat Feb91 "Nlfit Package"
+.ih
+NAME
+nlstat[ird] -- get an NLFIT parameter
+.ih
+SYNOPSIS
+include <math/nlfit.h>
+
+[ird]val = nlstat[ird] (nl, parameter)
+
+.nf
+pointer	nl		# curve descriptor
+int	parameter	# parameter to be returned
+.fi
+.ih
+ARGUMENTS
+.ls nl     
+The pointer to the curve descriptor structure.
+.le
+.ls parameter
+Parameter to be return.  Definitions in nlfit.h are:
+
+.nf
+NLNPARAMS	integer		# Number of parameters
+NLNFPARAMS	integer         # Number of fitted parameters
+NLITMAX		integer		# Maximum number of iterations
+NLITER		integer		# Current number of iterations
+NLNPTS		integer         # Number of points
+NLSUMSQ		real/double	# Current reduced chi-squared
+NLOLDSQ		real/double     # Previous reduced chi-squared
+NLLAMBDA	real/double     # Value of lambda factor
+NLTOL		real/double	# Fitting tolerance in %chi-squared
+NLSCATTER	real/double     # Mean scatter in the fit
+.fi
+.le
+.ih
+DESCRIPTION
+The values of integer, real or double  parameters are returned.
+The parameters include the number of parameters, the number of fitted
+parameters, the maximum number of iterations, the number of iterations,
+the total number of points, the current chi-squared, the previous
+chi-squared, the lambda factor, the fitting tolerance and the fitted
+scatter term.
+.ih
+EXAMPLES
+.nf
+include <math/nlfit.h>
+
+int	nlstati()
+
+call malloc (buf, nlstati (nl, NLNPARAMS), TY_REAL)
+call nlpgetr (nl, Memr[buf], nparams)
+.fi
+.ih
+SEE ALSO
+nlpget
+.endhelp
diff --git a/math/nlfit/doc/nlvector.hlp b/math/nlfit/doc/nlvector.hlp
new file mode 100644
index 00000000..98308cc3
--- /dev/null
+++ b/math/nlfit/doc/nlvector.hlp
@@ -0,0 +1,43 @@
+.help nlvector Feb91 "Nlfit Package"
+.ih
+NAME
+nlvector -- evaluate the fitted curve at a set of points
+.ih
+SYNOPSIS
+nlvector[rd] (nl, x, zfit, npts, nvars)
+
+.nf
+pointer		nl		# curve descriptor
+real/double	x[nvars*npts]	# array of variables stored as x[nvars,npts]
+real/double	zfit[npts]	# array of fitted function values
+int		npts		# number of data points
+int		nvars		# number of variables
+.fi
+.ih
+ARGUMENTS
+.ls nl   
+Pointer to the curve descriptor structure.
+.le
+.ls x     
+Array of variable values for all the data points.
+.le
+.ls zfit   
+Array of fitted function values.
+.le
+.ls npts   
+The number of data points at which the curve is to be evaluated.
+.le
+.ls nvars
+The number of variables in the fitted function.
+.le
+.ih
+DESCRIPTION
+Fit the curve to an array of data points.
+.ih
+NOTES
+NLVECTOR uses the coefficient array stored in the curfit descriptor
+structure.
+.ih
+SEE ALSO
+nleval
+.endhelp
diff --git a/math/nlfit/mkpkg b/math/nlfit/mkpkg
new file mode 100644
index 00000000..4573dfbf
--- /dev/null
+++ b/math/nlfit/mkpkg
@@ -0,0 +1,63 @@
+# The Non-linear Least-squares Fitting Package
+
+$checkout	libnlfit.a lib$
+$update		libnlfit.a
+$checkin	libnlfit.a lib$
+$exit
+
+tfiles:
+	$set	GEN = "$$generic -k -t rd"
+
+	$ifnewer (nlacpts.gx,  nlacptsr.x)	$(GEN) nlacpts.gx  $endif
+	$ifnewer (nlchomat.gx, nlchomatr.x)	$(GEN) nlchomat.gx $endif
+	$ifnewer (nldump.gx,   nldumpr.x)	$(GEN) nldump.gx   $endif
+	$ifnewer (nlerrors.gx, nlerrorsr.x)	$(GEN) nlerrors.gx $endif
+	$ifnewer (nleval.gx,   nlevalr.x)	$(GEN) nleval.gx   $endif
+	$ifnewer (nlfit.gx,    nlfitr.x)	$(GEN) nlfit.gx    $endif
+	$ifnewer (nlfitdef.gh, nlfitdefr.h)	$(GEN) nlfitdef.gh $endif
+	$ifnewer (nlfree.gx,   nlfreer.x)	$(GEN) nlfree.gx   $endif
+	$ifnewer (nlinit.gx,   nlinitr.x)	$(GEN) nlinit.gx   $endif
+	$ifnewer (nliter.gx,   nliterr.x)	$(GEN) nliter.gx   $endif
+	$ifnewer (nlpget.gx,   nlpgetr.x)	$(GEN) nlpget.gx   $endif
+	$ifnewer (nlsolve.gx,  nlsolver.x)	$(GEN) nlsolve.gx  $endif
+	$ifnewer (nlstat.gx,   nlstatr.x)	$(GEN) nlstat.gx   $endif
+	$ifnewer (nlvector.gx, nlvectorr.x)	$(GEN) nlvector.gx $endif
+	$ifnewer (nlzero.gx,   nlzeror.x)	$(GEN) nlzero.gx   $endif
+	;
+
+libnlfit.a:
+
+	$ifeq (USE_GENERIC, yes) $call tfiles $endif
+
+	nlacptsd.x	<math/nlfit.h>   "nlfitdefd.h"
+	nlacptsr.x	<math/nlfit.h>   "nlfitdefr.h"
+	nlchomatd.x	<math/nlfit.h>   "nlfitdefd.h" <mach.h>
+	nlchomatr.x	<math/nlfit.h>   "nlfitdefr.h" <mach.h>
+	nldumpd.x	"nlfitdefd.h"
+	nldumpr.x	"nlfitdefr.h"
+	nlerrmsg.x	<math/nlfit.h>
+	nlerrorsd.x	<math/nlfit.h>   "nlfitdefd.h" <mach.h>
+	nlerrorsr.x	<math/nlfit.h>   "nlfitdefr.h" <mach.h>
+	nlevald.x	<math/nlfit.h>   "nlfitdefd.h"
+	nlevalr.x	<math/nlfit.h>   "nlfitdefr.h"
+	nlfitd.x	<math/nlfit.h>   "nlfitdefd.h" <mach.h>
+	nlfitr.x	<math/nlfit.h>   "nlfitdefr.h" <mach.h>
+	nlfreed.x	"nlfitdefd.h"
+	nlfreer.x	"nlfitdefr.h"
+	nlinitd.x	"nlfitdefd.h"
+	nlinitr.x	"nlfitdefr.h"
+	nliterd.x	<math/nlfit.h>   "nlfitdefd.h" <mach.h>
+	nliterr.x	<math/nlfit.h>   "nlfitdefr.h" <mach.h>
+	nllist.x	
+	nlpgetd.x	"nlfitdefd.h"
+	nlpgetr.x	"nlfitdefr.h"
+	nlsolved.x	<math/nlfit.h>   "nlfitdefd.h"
+	nlsolver.x	<math/nlfit.h>   "nlfitdefr.h"
+	nlstati.x	<math/nlfit.h>   "nlfitdefr.h"
+	nlstatd.x	<math/nlfit.h>   "nlfitdefd.h"
+	nlstatr.x	<math/nlfit.h>   "nlfitdefr.h"
+	nlvectord.x	<math/nlfit.h>   "nlfitdefd.h"
+	nlvectorr.x	<math/nlfit.h>   "nlfitdefr.h"
+	nlzerod.x	"nlfitdefd.h"
+	nlzeror.x	"nlfitdefr.h"
+	;
diff --git a/math/nlfit/nlacpts.gx b/math/nlfit/nlacpts.gx
new file mode 100644
index 00000000..a9f947ba
--- /dev/null
+++ b/math/nlfit/nlacpts.gx
@@ -0,0 +1,111 @@
+include	<math/nlfit.h>
+$if (datatype == r)
+include "nlfitdefr.h"
+$else
+include "nlfitdefd.h"
+$endif
+
+define	VAR	(($1 - 1) * $2 + 1)
+
+# NLACPTS - Accumulate a series of data points.
+
+PIXEL procedure nlacpts$t (nl, x, z, w, npts, nvars)
+
+pointer	nl		# pointer to nl fitting structure
+PIXEL	x[ARB]		# independent variables (npts * nvars)
+PIXEL	z[ARB]		# function values (npts)
+PIXEL	w[ARB]		# weights (npts)
+int	npts		# number of points
+int	nvars		# number of independent variables
+
+int	i, nfree
+PIXEL	sum, z0, dz
+
+begin
+	# Zero the accumulators.
+	call aclr$t (ALPHA(NL_ALPHA(nl)), NL_NFPARAMS(nl) ** 2)
+	call aclr$t (BETA(NL_BETA(nl)), NL_NFPARAMS(nl))
+
+	# Accumulate the points into the fit.
+	NL_NPTS (nl) = npts
+	sum = PIXEL(0.0)
+	do i = 1, npts {
+	    call zcall7 (NL_DFUNC(nl), x[VAR (i, nvars)], nvars,
+		PARAM(NL_PARAM(nl)), DPARAM(NL_DELPARAM(nl)), NL_NPARAMS(nl),
+		z0, DERIV(NL_DERIV(nl)))
+	    dz = z[i] - z0
+	    call nl_accum$t (DERIV(NL_DERIV(nl)), PLIST(NL_PLIST(nl)),
+		w[i], dz, NL_NFPARAMS(nl), ALPHA(NL_ALPHA(nl)), 
+		BETA(NL_BETA(nl)))
+	    sum = sum + w[i] * dz * dz
+	}
+
+	# Return the reduced chisqr.
+	nfree = NL_NPTS(nl) - NL_NFPARAMS(nl)
+	if (nfree <= 0)
+	    return (PIXEL (0.0))
+	else
+	    return (sum / nfree)
+end
+
+
+# NLRESID -- Recompute the residuals
+
+PIXEL procedure nlresid$t (nl, x, z, w, npts, nvars)
+
+pointer	nl		# pointer to nl fitting structure
+PIXEL	x[ARB]		# independent variables (npts * nvars)
+PIXEL	z[ARB]		# function values (npts)
+PIXEL	w[ARB]		# weights (npts)
+int	npts		# number of points
+int	nvars		# number of independent variables
+
+int	i, nfree
+PIXEL	sum, z0, dz
+
+begin
+	# Accumulate the residuals.
+	NL_NPTS(nl) = npts
+	sum = PIXEL (0.0)
+	do i = 1, npts {
+	    call zcall5 (NL_FUNC(nl), x[VAR (i, nvars)], nvars, TRY(NL_TRY(nl)),
+		NL_NPARAMS(nl), z0)
+	    dz = z[i] - z0
+	    sum = sum + dz * dz * w[i]
+	}
+
+	# Compute the reduced chisqr.
+	nfree = NL_NPTS(nl) - NL_NFPARAMS(nl)
+	if (nfree <= 0)
+	    return (PIXEL (0.0))
+	else
+	    return (sum / nfree)
+end
+
+
+# NL_ACCUM -- Accumulate a single point into the fit
+
+procedure nl_accum$t (deriv, list, w, dz, nfit, alpha, beta)
+
+PIXEL	deriv[ARB]	# derivatives
+int	list[ARB]	# list of active parameters
+PIXEL	w		# weight
+PIXEL	dz		# difference between data and model
+int	nfit		# number of fitted parameters
+PIXEL	alpha[nfit,ARB]	# alpha matrix
+PIXEL	beta[nfit]	# beta matrix
+
+int	i, j, k
+PIXEL	wt
+
+begin
+	do i = 1, nfit {
+	    wt = deriv[list[i]] * w
+	    k = 1
+	    do j = i, nfit {
+	        alpha[k,i] = alpha[k,i] + wt * deriv[list[j]]
+		k = k + 1
+	    }
+	    beta[i] = beta[i] + dz * wt
+	}
+end
diff --git a/math/nlfit/nlacptsd.x b/math/nlfit/nlacptsd.x
new file mode 100644
index 00000000..b08a897f
--- /dev/null
+++ b/math/nlfit/nlacptsd.x
@@ -0,0 +1,107 @@
+include	<math/nlfit.h>
+include "nlfitdefd.h"
+
+define	VAR	(($1 - 1) * $2 + 1)
+
+# NLACPTS - Accumulate a series of data points.
+
+double procedure nlacptsd (nl, x, z, w, npts, nvars)
+
+pointer	nl		# pointer to nl fitting structure
+double	x[ARB]		# independent variables (npts * nvars)
+double	z[ARB]		# function values (npts)
+double	w[ARB]		# weights (npts)
+int	npts		# number of points
+int	nvars		# number of independent variables
+
+int	i, nfree
+double	sum, z0, dz
+
+begin
+	# Zero the accumulators.
+	call aclrd (ALPHA(NL_ALPHA(nl)), NL_NFPARAMS(nl) ** 2)
+	call aclrd (BETA(NL_BETA(nl)), NL_NFPARAMS(nl))
+
+	# Accumulate the points into the fit.
+	NL_NPTS (nl) = npts
+	sum = double(0.0)
+	do i = 1, npts {
+	    call zcall7 (NL_DFUNC(nl), x[VAR (i, nvars)], nvars,
+		PARAM(NL_PARAM(nl)), DPARAM(NL_DELPARAM(nl)), NL_NPARAMS(nl),
+		z0, DERIV(NL_DERIV(nl)))
+	    dz = z[i] - z0
+	    call nl_accumd (DERIV(NL_DERIV(nl)), PLIST(NL_PLIST(nl)),
+		w[i], dz, NL_NFPARAMS(nl), ALPHA(NL_ALPHA(nl)), 
+		BETA(NL_BETA(nl)))
+	    sum = sum + w[i] * dz * dz
+	}
+
+	# Return the reduced chisqr.
+	nfree = NL_NPTS(nl) - NL_NFPARAMS(nl)
+	if (nfree <= 0)
+	    return (double (0.0))
+	else
+	    return (sum / nfree)
+end
+
+
+# NLRESID -- Recompute the residuals
+
+double procedure nlresidd (nl, x, z, w, npts, nvars)
+
+pointer	nl		# pointer to nl fitting structure
+double	x[ARB]		# independent variables (npts * nvars)
+double	z[ARB]		# function values (npts)
+double	w[ARB]		# weights (npts)
+int	npts		# number of points
+int	nvars		# number of independent variables
+
+int	i, nfree
+double	sum, z0, dz
+
+begin
+	# Accumulate the residuals.
+	NL_NPTS(nl) = npts
+	sum = double (0.0)
+	do i = 1, npts {
+	    call zcall5 (NL_FUNC(nl), x[VAR (i, nvars)], nvars, TRY(NL_TRY(nl)),
+		NL_NPARAMS(nl), z0)
+	    dz = z[i] - z0
+	    sum = sum + dz * dz * w[i]
+	}
+
+	# Compute the reduced chisqr.
+	nfree = NL_NPTS(nl) - NL_NFPARAMS(nl)
+	if (nfree <= 0)
+	    return (double (0.0))
+	else
+	    return (sum / nfree)
+end
+
+
+# NL_ACCUM -- Accumulate a single point into the fit
+
+procedure nl_accumd (deriv, list, w, dz, nfit, alpha, beta)
+
+double	deriv[ARB]	# derivatives
+int	list[ARB]	# list of active parameters
+double	w		# weight
+double	dz		# difference between data and model
+int	nfit		# number of fitted parameters
+double	alpha[nfit,ARB]	# alpha matrix
+double	beta[nfit]	# beta matrix
+
+int	i, j, k
+double	wt
+
+begin
+	do i = 1, nfit {
+	    wt = deriv[list[i]] * w
+	    k = 1
+	    do j = i, nfit {
+	        alpha[k,i] = alpha[k,i] + wt * deriv[list[j]]
+		k = k + 1
+	    }
+	    beta[i] = beta[i] + dz * wt
+	}
+end
diff --git a/math/nlfit/nlacptsr.x b/math/nlfit/nlacptsr.x
new file mode 100644
index 00000000..d68daf72
--- /dev/null
+++ b/math/nlfit/nlacptsr.x
@@ -0,0 +1,107 @@
+include	<math/nlfit.h>
+include "nlfitdefr.h"
+
+define	VAR	(($1 - 1) * $2 + 1)
+
+# NLACPTS - Accumulate a series of data points.
+
+real procedure nlacptsr (nl, x, z, w, npts, nvars)
+
+pointer	nl		# pointer to nl fitting structure
+real	x[ARB]		# independent variables (npts * nvars)
+real	z[ARB]		# function values (npts)
+real	w[ARB]		# weights (npts)
+int	npts		# number of points
+int	nvars		# number of independent variables
+
+int	i, nfree
+real	sum, z0, dz
+
+begin
+	# Zero the accumulators.
+	call aclrr (ALPHA(NL_ALPHA(nl)), NL_NFPARAMS(nl) ** 2)
+	call aclrr (BETA(NL_BETA(nl)), NL_NFPARAMS(nl))
+
+	# Accumulate the points into the fit.
+	NL_NPTS (nl) = npts
+	sum = real(0.0)
+	do i = 1, npts {
+	    call zcall7 (NL_DFUNC(nl), x[VAR (i, nvars)], nvars,
+		PARAM(NL_PARAM(nl)), DPARAM(NL_DELPARAM(nl)), NL_NPARAMS(nl),
+		z0, DERIV(NL_DERIV(nl)))
+	    dz = z[i] - z0
+	    call nl_accumr (DERIV(NL_DERIV(nl)), PLIST(NL_PLIST(nl)),
+		w[i], dz, NL_NFPARAMS(nl), ALPHA(NL_ALPHA(nl)), 
+		BETA(NL_BETA(nl)))
+	    sum = sum + w[i] * dz * dz
+	}
+
+	# Return the reduced chisqr.
+	nfree = NL_NPTS(nl) - NL_NFPARAMS(nl)
+	if (nfree <= 0)
+	    return (real (0.0))
+	else
+	    return (sum / nfree)
+end
+
+
+# NLRESID -- Recompute the residuals
+
+real procedure nlresidr (nl, x, z, w, npts, nvars)
+
+pointer	nl		# pointer to nl fitting structure
+real	x[ARB]		# independent variables (npts * nvars)
+real	z[ARB]		# function values (npts)
+real	w[ARB]		# weights (npts)
+int	npts		# number of points
+int	nvars		# number of independent variables
+
+int	i, nfree
+real	sum, z0, dz
+
+begin
+	# Accumulate the residuals.
+	NL_NPTS(nl) = npts
+	sum = real (0.0)
+	do i = 1, npts {
+	    call zcall5 (NL_FUNC(nl), x[VAR (i, nvars)], nvars, TRY(NL_TRY(nl)),
+		NL_NPARAMS(nl), z0)
+	    dz = z[i] - z0
+	    sum = sum + dz * dz * w[i]
+	}
+
+	# Compute the reduced chisqr.
+	nfree = NL_NPTS(nl) - NL_NFPARAMS(nl)
+	if (nfree <= 0)
+	    return (real (0.0))
+	else
+	    return (sum / nfree)
+end
+
+
+# NL_ACCUM -- Accumulate a single point into the fit
+
+procedure nl_accumr (deriv, list, w, dz, nfit, alpha, beta)
+
+real	deriv[ARB]	# derivatives
+int	list[ARB]	# list of active parameters
+real	w		# weight
+real	dz		# difference between data and model
+int	nfit		# number of fitted parameters
+real	alpha[nfit,ARB]	# alpha matrix
+real	beta[nfit]	# beta matrix
+
+int	i, j, k
+real	wt
+
+begin
+	do i = 1, nfit {
+	    wt = deriv[list[i]] * w
+	    k = 1
+	    do j = i, nfit {
+	        alpha[k,i] = alpha[k,i] + wt * deriv[list[j]]
+		k = k + 1
+	    }
+	    beta[i] = beta[i] + dz * wt
+	}
+end
diff --git a/math/nlfit/nlchomat.gx b/math/nlfit/nlchomat.gx
new file mode 100644
index 00000000..5a36c63c
--- /dev/null
+++ b/math/nlfit/nlchomat.gx
@@ -0,0 +1,130 @@
+include <mach.h>
+include <math/nlfit.h>
+$if (datatype == r)
+include "nlfitdefr.h"
+$else
+include "nlfitdefd.h"
+$endif
+
+
+# NL_CHFAC -- 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 nl_chfac$t (matrix, nbands, nrows, matfac, ier)
+
+PIXEL   matrix[nbands, nrows]	# data matrix
+int	nbands			# number of bands
+int	nrows			# number of rows
+PIXEL   matfac[nbands, nrows]	# Cholesky factorization
+int	ier			# error code
+
+int	i, n, j, imax, jmax
+PIXEL   ratio
+
+begin
+	# Test for a single element matrix.
+	if (nrows == 1) {
+	    if (matrix[1,1] > PIXEL (0.0))
+	        matfac[1,1] = 1. / matrix[1,1]
+	    return
+	}
+		
+	# Copy the original matrix into matfac.
+	do n = 1, nrows {
+	    do j = 1, nbands
+		matfac[j,n] = matrix[j,n]
+	}
+
+	# Compute the factorization of the matrix.
+	do n = 1, nrows {
+
+	    # Test to see if matrix is singular.
+	    if (((matfac[1,n] + matrix[1,n]) - matrix[1,n]) <= EPSILON$T) {
+	    #if (((matfac[1,n] + matrix[1,n]) - matrix[1,n]) <= PIXEL(0.0)) {
+		do j = 1, nbands
+		    matfac[j,n] = PIXEL (0.0)
+		ier = SINGULAR
+		next
+	    }
+
+	    matfac[1,n] = PIXEL (1.0) / 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
+
+
+# NL_CHSLV -- Solve the matrix whose Cholesky factorization was calculated in
+# NL_CHFAC for the coefficients. This routine was adapted from bchslv.f
+# described in "A Practical Guide to Splines", by Carl de Boor (1978).
+
+procedure nl_chslv$t (matfac, nbands, nrows, vector, coeff)
+
+PIXEL	matfac[nbands,nrows] 		# Cholesky factorization
+int	nbands				# number of bands
+int	nrows				# number of rows
+PIXEL	vector[nrows]			# right side of matrix equation
+PIXEL	coeff[nrows]			# coefficients
+
+int	i, n, j, jmax, nbndm1
+
+begin
+	# Test for a single element matrix.
+	if (nrows == 1) {
+	    coeff[1] = vector[1] * matfac[1,1]
+	    return
+	}
+
+	# Copy input vector to coefficients vector.
+	do i = 1, nrows
+	    coeff[i] = vector[i]
+
+	# Perform 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]
+	    }
+	}
+
+	# Perform backward 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
+
+
+# NL_DAMP -- Procedure to add damping to matrix
+
+procedure nl_damp$t (inmatrix, outmatrix, constant, nbands, nrows)
+
+PIXEL	inmatrix[nbands,ARB]		# input matrix
+PIXEL	outmatrix[nbands,ARB]		# output matrix
+PIXEL	constant			# damping constant
+int	nbands, nrows			# dimensions of matrix
+
+int	i
+
+begin
+	do i = 1, nrows
+	    outmatrix[1,i] = inmatrix[1,i] * constant
+end
diff --git a/math/nlfit/nlchomatd.x b/math/nlfit/nlchomatd.x
new file mode 100644
index 00000000..775793ac
--- /dev/null
+++ b/math/nlfit/nlchomatd.x
@@ -0,0 +1,126 @@
+include <mach.h>
+include <math/nlfit.h>
+include "nlfitdefd.h"
+
+
+# NL_CHFAC -- 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 nl_chfacd (matrix, nbands, nrows, matfac, ier)
+
+double   matrix[nbands, nrows]	# data matrix
+int	nbands			# number of bands
+int	nrows			# number of rows
+double   matfac[nbands, nrows]	# Cholesky factorization
+int	ier			# error code
+
+int	i, n, j, imax, jmax
+double   ratio
+
+begin
+	# Test for a single element matrix.
+	if (nrows == 1) {
+	    if (matrix[1,1] > double (0.0))
+	        matfac[1,1] = 1. / matrix[1,1]
+	    return
+	}
+		
+	# Copy the original matrix into matfac.
+	do n = 1, nrows {
+	    do j = 1, nbands
+		matfac[j,n] = matrix[j,n]
+	}
+
+	# Compute the factorization of the matrix.
+	do n = 1, nrows {
+
+	    # Test to see if matrix is singular.
+	    if (((matfac[1,n] + matrix[1,n]) - matrix[1,n]) <= EPSILOND) {
+	    #if (((matfac[1,n] + matrix[1,n]) - matrix[1,n]) <= PIXEL(0.0)) {
+		do j = 1, nbands
+		    matfac[j,n] = double (0.0)
+		ier = SINGULAR
+		next
+	    }
+
+	    matfac[1,n] = double (1.0) / 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
+
+
+# NL_CHSLV -- Solve the matrix whose Cholesky factorization was calculated in
+# NL_CHFAC for the coefficients. This routine was adapted from bchslv.f
+# described in "A Practical Guide to Splines", by Carl de Boor (1978).
+
+procedure nl_chslvd (matfac, nbands, nrows, vector, coeff)
+
+double	matfac[nbands,nrows] 		# Cholesky factorization
+int	nbands				# number of bands
+int	nrows				# number of rows
+double	vector[nrows]			# right side of matrix equation
+double	coeff[nrows]			# coefficients
+
+int	i, n, j, jmax, nbndm1
+
+begin
+	# Test for a single element matrix.
+	if (nrows == 1) {
+	    coeff[1] = vector[1] * matfac[1,1]
+	    return
+	}
+
+	# Copy input vector to coefficients vector.
+	do i = 1, nrows
+	    coeff[i] = vector[i]
+
+	# Perform 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]
+	    }
+	}
+
+	# Perform backward 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
+
+
+# NL_DAMP -- Procedure to add damping to matrix
+
+procedure nl_dampd (inmatrix, outmatrix, constant, nbands, nrows)
+
+double	inmatrix[nbands,ARB]		# input matrix
+double	outmatrix[nbands,ARB]		# output matrix
+double	constant			# damping constant
+int	nbands, nrows			# dimensions of matrix
+
+int	i
+
+begin
+	do i = 1, nrows
+	    outmatrix[1,i] = inmatrix[1,i] * constant
+end
diff --git a/math/nlfit/nlchomatr.x b/math/nlfit/nlchomatr.x
new file mode 100644
index 00000000..09f3a7b3
--- /dev/null
+++ b/math/nlfit/nlchomatr.x
@@ -0,0 +1,126 @@
+include <mach.h>
+include <math/nlfit.h>
+include "nlfitdefr.h"
+
+
+# NL_CHFAC -- 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 nl_chfacr (matrix, nbands, nrows, matfac, ier)
+
+real   matrix[nbands, nrows]	# data matrix
+int	nbands			# number of bands
+int	nrows			# number of rows
+real   matfac[nbands, nrows]	# Cholesky factorization
+int	ier			# error code
+
+int	i, n, j, imax, jmax
+real   ratio
+
+begin
+	# Test for a single element matrix.
+	if (nrows == 1) {
+	    if (matrix[1,1] > real (0.0))
+	        matfac[1,1] = 1. / matrix[1,1]
+	    return
+	}
+		
+	# Copy the original matrix into matfac.
+	do n = 1, nrows {
+	    do j = 1, nbands
+		matfac[j,n] = matrix[j,n]
+	}
+
+	# Compute the factorization of the matrix.
+	do n = 1, nrows {
+
+	    # Test to see if matrix is singular.
+	    if (((matfac[1,n] + matrix[1,n]) - matrix[1,n]) <= EPSILONR) {
+	    #if (((matfac[1,n] + matrix[1,n]) - matrix[1,n]) <= PIXEL(0.0)) {
+		do j = 1, nbands
+		    matfac[j,n] = real (0.0)
+		ier = SINGULAR
+		next
+	    }
+
+	    matfac[1,n] = real (1.0) / 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
+
+
+# NL_CHSLV -- Solve the matrix whose Cholesky factorization was calculated in
+# NL_CHFAC for the coefficients. This routine was adapted from bchslv.f
+# described in "A Practical Guide to Splines", by Carl de Boor (1978).
+
+procedure nl_chslvr (matfac, nbands, nrows, vector, coeff)
+
+real	matfac[nbands,nrows] 		# Cholesky factorization
+int	nbands				# number of bands
+int	nrows				# number of rows
+real	vector[nrows]			# right side of matrix equation
+real	coeff[nrows]			# coefficients
+
+int	i, n, j, jmax, nbndm1
+
+begin
+	# Test for a single element matrix.
+	if (nrows == 1) {
+	    coeff[1] = vector[1] * matfac[1,1]
+	    return
+	}
+
+	# Copy input vector to coefficients vector.
+	do i = 1, nrows
+	    coeff[i] = vector[i]
+
+	# Perform 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]
+	    }
+	}
+
+	# Perform backward 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
+
+
+# NL_DAMP -- Procedure to add damping to matrix
+
+procedure nl_dampr (inmatrix, outmatrix, constant, nbands, nrows)
+
+real	inmatrix[nbands,ARB]		# input matrix
+real	outmatrix[nbands,ARB]		# output matrix
+real	constant			# damping constant
+int	nbands, nrows			# dimensions of matrix
+
+int	i
+
+begin
+	do i = 1, nrows
+	    outmatrix[1,i] = inmatrix[1,i] * constant
+end
diff --git a/math/nlfit/nldump.gx b/math/nlfit/nldump.gx
new file mode 100644
index 00000000..aa5899bd
--- /dev/null
+++ b/math/nlfit/nldump.gx
@@ -0,0 +1,164 @@
+$if (datatype == r)
+include "nlfitdefr.h"
+$else
+include "nlfitdefd.h"
+$endif
+
+# NL_DUMP -- Dump NLFIT structure to a file
+
+procedure nl_dump$t (fd, nl)
+
+int	fd		# file descriptor
+pointer	nl		# NLFIT descriptor
+
+int	i, npars, nfpars
+
+begin
+	# Test NLFIT pointer
+	if (nl == NULL) {
+	    call fprintf (fd, "\n****** nl_dump: Null NLFIT pointer\n")
+	    call flush (fd)
+	    return
+	}
+
+	# File and NLFIT descriptors
+	call fprintf (fd, "\n****** nl_dump: (fd=%d), (nl=%d)\n")
+	    call pargi (fd)
+	    call pargi (nl)
+	call flush (fd)
+
+	# Dump function and derivative addresses
+	call fprintf (fd, "Fitting function pointer    = %d\n")
+	    call pargi (NL_FUNC (nl))
+	call fprintf (fd, "Derivative function pointer = %d\n")
+	    call pargi (NL_DFUNC (nl))
+	call flush (fd)
+
+	# Number of parameters
+	npars = NL_NPARAMS (nl)
+	nfpars = NL_NFPARAMS (nl)
+	call fprintf (fd, "Number of parameters        = %d\n")
+	    call pargi (npars)
+	call fprintf (fd, "Number of fitted parameters = %d\n")
+	    call pargi (nfpars)
+	call flush (fd)
+
+	# Fit parameters
+	call fprintf (fd, "Max number of iterations    = %d\n")
+	    call pargi (NL_ITMAX (nl))
+	call fprintf (fd, "Tolerance for convergence   = %g\n")
+	    call parg$t (NL_TOL (nl))
+	call flush (fd)
+
+	# Sums
+	call fprintf (fd, "Damping factor = %g\n")
+	    call parg$t (NL_LAMBDA (nl))
+	call fprintf (fd, "Sum of residuals squared last iteration = %g\n")
+	    call parg$t (NL_OLDSQ (nl))
+	call fprintf (fd, "Sum of residuals squared                = %g\n")
+	    call parg$t (NL_SUMSQ (nl))
+	call flush (fd)
+
+	# Counters
+	call fprintf (fd, "Iteration counter       = %d\n")
+	    call pargi (NL_ITER (nl))
+	call fprintf (fd, "Number of points in fit = %d\n")
+	    call pargi (NL_NPTS (nl))
+	call flush (fd)
+
+	# Parameter values
+	call fprintf (fd, "Parameter values (%d):\n")
+	    call pargi (NL_PARAM (nl))
+	if (NL_PARAM (nl) != NULL) {
+	    do i = 1, npars {
+		call fprintf (fd, "%d -> %g\n")
+		    call pargi (i)
+		    call parg$t (Mem$t[NL_PARAM (nl) + i - 1])
+	    }
+	} else
+	    call fprintf (fd, "	Null pointer\n")
+	call flush (fd)
+
+	# Parameter errors
+	call fprintf (fd, "Parameter errors (%d):\n")
+	    call pargi (NL_DPARAM (nl))
+	if (NL_DPARAM (nl) != NULL) {
+	    do i = 1, npars {
+		call fprintf (fd, "%d -> %g\n")
+		    call pargi (i)
+		    call parg$t (Mem$t[NL_DPARAM (nl) + i - 1])
+	    }
+	} else
+	    call fprintf (fd, "	Null pointer\n")
+	call flush (fd)
+
+	# Parameter list
+	call fprintf (fd, "Parameter list (%d):\n")
+	    call pargi (NL_PLIST (nl))
+	if (NL_PLIST (nl) != NULL) {
+	    do i = 1, npars {
+		call fprintf (fd, "%d -> %d\n")
+		    call pargi (i)
+		    call pargi (Memi[NL_PLIST (nl) + i - 1])
+	    }
+	} else
+	    call fprintf (fd, "	Null pointer\n")
+	call flush (fd)
+
+	# Alpha matrix
+	call fprintf (fd, "Alpha matrix (%d):\n")
+	    call pargi (NL_ALPHA (nl))
+	if (NL_ALPHA (nl) != NULL)
+	    call nl_adump$t (fd, Mem$t[NL_ALPHA (nl)], nfpars, nfpars)
+	else
+	    call fprintf (fd, "	Null pointer\n")
+	call flush (fd)
+
+	# Beta matrix
+	call fprintf (fd, "Beta matrix (%d):\n")
+	    call pargi (NL_BETA (nl))
+	if (NL_BETA (nl) != NULL)
+	    call nl_adump$t (fd, Mem$t[NL_BETA (nl)], nfpars, 1)
+	else
+	    call fprintf (fd, "	Null pointer\n")
+	call flush (fd)
+
+	# Covariance matrix
+	call fprintf (fd, "Covariance matrix (%d):\n")
+	    call pargi (NL_COVAR (nl))
+	if (NL_COVAR (nl) != NULL)
+	    call nl_adump$t (fd, Mem$t[NL_COVAR (nl)], nfpars, nfpars)
+	else
+	    call fprintf (fd, "	Null pointer\n")
+	call flush (fd)
+
+	# Cholesky factorization
+	call fprintf (fd, "Cholesky factorization matrix (%d):\n")
+	    call pargi (NL_CHOFAC (nl))
+	if (NL_CHOFAC (nl) != NULL)
+	    call nl_adump$t (fd, Mem$t[NL_CHOFAC (nl)], nfpars, nfpars)
+	else
+	    call fprintf (fd, "	Null pointer\n")
+	call flush (fd)
+end
+
+
+# NL_ADUMP -- Dump array to file
+
+procedure nl_adump$t (fd, a, nrows, ncols)
+
+int	fd			# file descriptor
+PIXEL	a[nrows, ncols]		# array
+int	nrows, ncols		# dimension
+
+int	i, j
+
+begin
+	do i = 1, nrows {
+	    do j = 1, ncols {
+		call fprintf (fd, "%g  ")
+		    call parg$t (a[i, j])
+	    }
+	    call fprintf (fd, "\n")
+	}
+end
diff --git a/math/nlfit/nldumpd.x b/math/nlfit/nldumpd.x
new file mode 100644
index 00000000..f4c31d7d
--- /dev/null
+++ b/math/nlfit/nldumpd.x
@@ -0,0 +1,160 @@
+include "nlfitdefd.h"
+
+# NL_DUMP -- Dump NLFIT structure to a file
+
+procedure nl_dumpd (fd, nl)
+
+int	fd		# file descriptor
+pointer	nl		# NLFIT descriptor
+
+int	i, npars, nfpars
+
+begin
+	# Test NLFIT pointer
+	if (nl == NULL) {
+	    call fprintf (fd, "\n****** nl_dump: Null NLFIT pointer\n")
+	    call flush (fd)
+	    return
+	}
+
+	# File and NLFIT descriptors
+	call fprintf (fd, "\n****** nl_dump: (fd=%d), (nl=%d)\n")
+	    call pargi (fd)
+	    call pargi (nl)
+	call flush (fd)
+
+	# Dump function and derivative addresses
+	call fprintf (fd, "Fitting function pointer    = %d\n")
+	    call pargi (NL_FUNC (nl))
+	call fprintf (fd, "Derivative function pointer = %d\n")
+	    call pargi (NL_DFUNC (nl))
+	call flush (fd)
+
+	# Number of parameters
+	npars = NL_NPARAMS (nl)
+	nfpars = NL_NFPARAMS (nl)
+	call fprintf (fd, "Number of parameters        = %d\n")
+	    call pargi (npars)
+	call fprintf (fd, "Number of fitted parameters = %d\n")
+	    call pargi (nfpars)
+	call flush (fd)
+
+	# Fit parameters
+	call fprintf (fd, "Max number of iterations    = %d\n")
+	    call pargi (NL_ITMAX (nl))
+	call fprintf (fd, "Tolerance for convergence   = %g\n")
+	    call pargd (NL_TOL (nl))
+	call flush (fd)
+
+	# Sums
+	call fprintf (fd, "Damping factor = %g\n")
+	    call pargd (NL_LAMBDA (nl))
+	call fprintf (fd, "Sum of residuals squared last iteration = %g\n")
+	    call pargd (NL_OLDSQ (nl))
+	call fprintf (fd, "Sum of residuals squared                = %g\n")
+	    call pargd (NL_SUMSQ (nl))
+	call flush (fd)
+
+	# Counters
+	call fprintf (fd, "Iteration counter       = %d\n")
+	    call pargi (NL_ITER (nl))
+	call fprintf (fd, "Number of points in fit = %d\n")
+	    call pargi (NL_NPTS (nl))
+	call flush (fd)
+
+	# Parameter values
+	call fprintf (fd, "Parameter values (%d):\n")
+	    call pargi (NL_PARAM (nl))
+	if (NL_PARAM (nl) != NULL) {
+	    do i = 1, npars {
+		call fprintf (fd, "%d -> %g\n")
+		    call pargi (i)
+		    call pargd (Memd[NL_PARAM (nl) + i - 1])
+	    }
+	} else
+	    call fprintf (fd, "	Null pointer\n")
+	call flush (fd)
+
+	# Parameter errors
+	call fprintf (fd, "Parameter errors (%d):\n")
+	    call pargi (NL_DPARAM (nl))
+	if (NL_DPARAM (nl) != NULL) {
+	    do i = 1, npars {
+		call fprintf (fd, "%d -> %g\n")
+		    call pargi (i)
+		    call pargd (Memd[NL_DPARAM (nl) + i - 1])
+	    }
+	} else
+	    call fprintf (fd, "	Null pointer\n")
+	call flush (fd)
+
+	# Parameter list
+	call fprintf (fd, "Parameter list (%d):\n")
+	    call pargi (NL_PLIST (nl))
+	if (NL_PLIST (nl) != NULL) {
+	    do i = 1, npars {
+		call fprintf (fd, "%d -> %d\n")
+		    call pargi (i)
+		    call pargi (Memi[NL_PLIST (nl) + i - 1])
+	    }
+	} else
+	    call fprintf (fd, "	Null pointer\n")
+	call flush (fd)
+
+	# Alpha matrix
+	call fprintf (fd, "Alpha matrix (%d):\n")
+	    call pargi (NL_ALPHA (nl))
+	if (NL_ALPHA (nl) != NULL)
+	    call nl_adumpd (fd, Memd[NL_ALPHA (nl)], nfpars, nfpars)
+	else
+	    call fprintf (fd, "	Null pointer\n")
+	call flush (fd)
+
+	# Beta matrix
+	call fprintf (fd, "Beta matrix (%d):\n")
+	    call pargi (NL_BETA (nl))
+	if (NL_BETA (nl) != NULL)
+	    call nl_adumpd (fd, Memd[NL_BETA (nl)], nfpars, 1)
+	else
+	    call fprintf (fd, "	Null pointer\n")
+	call flush (fd)
+
+	# Covariance matrix
+	call fprintf (fd, "Covariance matrix (%d):\n")
+	    call pargi (NL_COVAR (nl))
+	if (NL_COVAR (nl) != NULL)
+	    call nl_adumpd (fd, Memd[NL_COVAR (nl)], nfpars, nfpars)
+	else
+	    call fprintf (fd, "	Null pointer\n")
+	call flush (fd)
+
+	# Cholesky factorization
+	call fprintf (fd, "Cholesky factorization matrix (%d):\n")
+	    call pargi (NL_CHOFAC (nl))
+	if (NL_CHOFAC (nl) != NULL)
+	    call nl_adumpd (fd, Memd[NL_CHOFAC (nl)], nfpars, nfpars)
+	else
+	    call fprintf (fd, "	Null pointer\n")
+	call flush (fd)
+end
+
+
+# NL_ADUMP -- Dump array to file
+
+procedure nl_adumpd (fd, a, nrows, ncols)
+
+int	fd			# file descriptor
+double	a[nrows, ncols]		# array
+int	nrows, ncols		# dimension
+
+int	i, j
+
+begin
+	do i = 1, nrows {
+	    do j = 1, ncols {
+		call fprintf (fd, "%g  ")
+		    call pargd (a[i, j])
+	    }
+	    call fprintf (fd, "\n")
+	}
+end
diff --git a/math/nlfit/nldumpr.x b/math/nlfit/nldumpr.x
new file mode 100644
index 00000000..a81c6366
--- /dev/null
+++ b/math/nlfit/nldumpr.x
@@ -0,0 +1,160 @@
+include "nlfitdefr.h"
+
+# NL_DUMP -- Dump NLFIT structure to a file
+
+procedure nl_dumpr (fd, nl)
+
+int	fd		# file descriptor
+pointer	nl		# NLFIT descriptor
+
+int	i, npars, nfpars
+
+begin
+	# Test NLFIT pointer
+	if (nl == NULL) {
+	    call fprintf (fd, "\n****** nl_dump: Null NLFIT pointer\n")
+	    call flush (fd)
+	    return
+	}
+
+	# File and NLFIT descriptors
+	call fprintf (fd, "\n****** nl_dump: (fd=%d), (nl=%d)\n")
+	    call pargi (fd)
+	    call pargi (nl)
+	call flush (fd)
+
+	# Dump function and derivative addresses
+	call fprintf (fd, "Fitting function pointer    = %d\n")
+	    call pargi (NL_FUNC (nl))
+	call fprintf (fd, "Derivative function pointer = %d\n")
+	    call pargi (NL_DFUNC (nl))
+	call flush (fd)
+
+	# Number of parameters
+	npars = NL_NPARAMS (nl)
+	nfpars = NL_NFPARAMS (nl)
+	call fprintf (fd, "Number of parameters        = %d\n")
+	    call pargi (npars)
+	call fprintf (fd, "Number of fitted parameters = %d\n")
+	    call pargi (nfpars)
+	call flush (fd)
+
+	# Fit parameters
+	call fprintf (fd, "Max number of iterations    = %d\n")
+	    call pargi (NL_ITMAX (nl))
+	call fprintf (fd, "Tolerance for convergence   = %g\n")
+	    call pargr (NL_TOL (nl))
+	call flush (fd)
+
+	# Sums
+	call fprintf (fd, "Damping factor = %g\n")
+	    call pargr (NL_LAMBDA (nl))
+	call fprintf (fd, "Sum of residuals squared last iteration = %g\n")
+	    call pargr (NL_OLDSQ (nl))
+	call fprintf (fd, "Sum of residuals squared                = %g\n")
+	    call pargr (NL_SUMSQ (nl))
+	call flush (fd)
+
+	# Counters
+	call fprintf (fd, "Iteration counter       = %d\n")
+	    call pargi (NL_ITER (nl))
+	call fprintf (fd, "Number of points in fit = %d\n")
+	    call pargi (NL_NPTS (nl))
+	call flush (fd)
+
+	# Parameter values
+	call fprintf (fd, "Parameter values (%d):\n")
+	    call pargi (NL_PARAM (nl))
+	if (NL_PARAM (nl) != NULL) {
+	    do i = 1, npars {
+		call fprintf (fd, "%d -> %g\n")
+		    call pargi (i)
+		    call pargr (Memr[NL_PARAM (nl) + i - 1])
+	    }
+	} else
+	    call fprintf (fd, "	Null pointer\n")
+	call flush (fd)
+
+	# Parameter errors
+	call fprintf (fd, "Parameter errors (%d):\n")
+	    call pargi (NL_DPARAM (nl))
+	if (NL_DPARAM (nl) != NULL) {
+	    do i = 1, npars {
+		call fprintf (fd, "%d -> %g\n")
+		    call pargi (i)
+		    call pargr (Memr[NL_DPARAM (nl) + i - 1])
+	    }
+	} else
+	    call fprintf (fd, "	Null pointer\n")
+	call flush (fd)
+
+	# Parameter list
+	call fprintf (fd, "Parameter list (%d):\n")
+	    call pargi (NL_PLIST (nl))
+	if (NL_PLIST (nl) != NULL) {
+	    do i = 1, npars {
+		call fprintf (fd, "%d -> %d\n")
+		    call pargi (i)
+		    call pargi (Memi[NL_PLIST (nl) + i - 1])
+	    }
+	} else
+	    call fprintf (fd, "	Null pointer\n")
+	call flush (fd)
+
+	# Alpha matrix
+	call fprintf (fd, "Alpha matrix (%d):\n")
+	    call pargi (NL_ALPHA (nl))
+	if (NL_ALPHA (nl) != NULL)
+	    call nl_adumpr (fd, Memr[NL_ALPHA (nl)], nfpars, nfpars)
+	else
+	    call fprintf (fd, "	Null pointer\n")
+	call flush (fd)
+
+	# Beta matrix
+	call fprintf (fd, "Beta matrix (%d):\n")
+	    call pargi (NL_BETA (nl))
+	if (NL_BETA (nl) != NULL)
+	    call nl_adumpr (fd, Memr[NL_BETA (nl)], nfpars, 1)
+	else
+	    call fprintf (fd, "	Null pointer\n")
+	call flush (fd)
+
+	# Covariance matrix
+	call fprintf (fd, "Covariance matrix (%d):\n")
+	    call pargi (NL_COVAR (nl))
+	if (NL_COVAR (nl) != NULL)
+	    call nl_adumpr (fd, Memr[NL_COVAR (nl)], nfpars, nfpars)
+	else
+	    call fprintf (fd, "	Null pointer\n")
+	call flush (fd)
+
+	# Cholesky factorization
+	call fprintf (fd, "Cholesky factorization matrix (%d):\n")
+	    call pargi (NL_CHOFAC (nl))
+	if (NL_CHOFAC (nl) != NULL)
+	    call nl_adumpr (fd, Memr[NL_CHOFAC (nl)], nfpars, nfpars)
+	else
+	    call fprintf (fd, "	Null pointer\n")
+	call flush (fd)
+end
+
+
+# NL_ADUMP -- Dump array to file
+
+procedure nl_adumpr (fd, a, nrows, ncols)
+
+int	fd			# file descriptor
+real	a[nrows, ncols]		# array
+int	nrows, ncols		# dimension
+
+int	i, j
+
+begin
+	do i = 1, nrows {
+	    do j = 1, ncols {
+		call fprintf (fd, "%g  ")
+		    call pargr (a[i, j])
+	    }
+	    call fprintf (fd, "\n")
+	}
+end
diff --git a/math/nlfit/nlerrmsg.x b/math/nlfit/nlerrmsg.x
new file mode 100644
index 00000000..e2b04b73
--- /dev/null
+++ b/math/nlfit/nlerrmsg.x
@@ -0,0 +1,24 @@
+include	<math/nlfit.h>
+
+# NLERRMSG -- Convert NLFIT error code into an error message.
+
+procedure nlerrmsg (ier, errmsg, maxch)
+
+int	ier			# NLFIT error code
+char	errmsg[maxch]		# output error message 
+int	maxch			# maximum number of chars
+
+begin
+	switch (ier) {
+	case DONE:
+	    call strcpy ("Solution converged", errmsg, maxch)
+	case SINGULAR:
+	    call strcpy ("Singular matrix", errmsg, maxch)
+	case NO_DEG_FREEDOM:
+	    call strcpy ("Too few points", errmsg, maxch)
+	case NOT_DONE:
+	    call strcpy ("Solution did not converge", errmsg, maxch)
+	default:
+	    call strcpy ("Unknown error code", errmsg, maxch)
+	}
+end
diff --git a/math/nlfit/nlerrors.gx b/math/nlfit/nlerrors.gx
new file mode 100644
index 00000000..061d2214
--- /dev/null
+++ b/math/nlfit/nlerrors.gx
@@ -0,0 +1,111 @@
+include <mach.h>
+include <math/nlfit.h>
+$if (datatype == r)
+include "nlfitdefr.h"
+$else
+include "nlfitdefd.h"
+$endif
+
+define	COV	Mem$t[P2P($1)]	# element of COV
+
+
+# NLERRORS -- Procedure to calculate the reduced chi-squared of the fit
+# and the errors of the coefficients. First the variance
+# and the reduced chi-squared of the fit are estimated. If these two
+# quantities are identical the variance is used to scale the errors
+# in the coefficients. The errors in the coefficients are proportional
+# to the inverse diagonal elements of MATRIX.
+
+procedure nlerrors$t (nl, z, zfit, w, npts, variance, chisqr, errors)
+
+pointer	nl		# curve descriptor
+PIXEL	z[ARB]		# data points
+PIXEL	zfit[ARB]	# fitted data points
+PIXEL	w[ARB]		# array of weights
+int	npts		# number of points
+PIXEL	variance	# variance of the fit
+PIXEL	chisqr		# reduced chi-squared of fit (output)
+PIXEL	errors[ARB]	# errors in coefficients (output)
+
+int	i, n, nfree
+pointer	sp, covptr
+PIXEL   factor
+
+begin
+	# Allocate space for covariance vector.
+	call smark (sp)
+	call salloc (covptr, NL_NPARAMS(nl), TY_PIXEL)
+
+	# Estimate the variance and reduce chi-squared of the fit.
+	n = 0
+	variance = PIXEL (0.0)
+	chisqr = PIXEL (0.0)
+	do i = 1, npts {
+	    if (w[i] <= PIXEL (0.0))
+	        next
+	    factor = (z[i] - zfit[i]) ** 2
+	    variance = variance + factor
+	    chisqr = chisqr + factor * w[i]
+	    n = n + 1
+	}
+
+	# Calculate the reduced chi-squared.
+	nfree = n - NL_NFPARAMS(nl)
+	if (nfree  > 0) {
+	    variance = variance / nfree
+	    chisqr = chisqr / nfree
+	} else {
+	    variance = PIXEL (0.0)
+	    chisqr = PIXEL (0.0)
+	}
+
+	# If the variance equals the reduced chi_squared as in the case
+	# of uniform weights then scale the errors in the coefficients
+	# by the variance and not the reduced chi-squared
+
+	if (abs (chisqr - variance) <= EPSILON$T) {
+	    if (nfree > 0)
+		factor = chisqr
+	    else
+		factor = PIXEL (0.0)
+	} else
+	    factor = 1.
+
+	# Calculate the  errors in the coefficients.
+	call aclr$t (errors, NL_NPARAMS(nl))
+	call nlinvert$t (ALPHA(NL_ALPHA(nl)), CHOFAC(NL_CHOFAC(nl)),
+	    COV(covptr), errors, PLIST(NL_PLIST(nl)), NL_NFPARAMS(nl), factor)
+
+	call sfree (sp)
+end
+
+
+# NLINVERT -- Procedure to invert matrix and compute errors
+
+procedure nlinvert$t (alpha, chofac, cov, errors, list, nfit, variance)
+
+PIXEL	alpha[nfit,ARB]		# data matrix
+PIXEL	chofac[nfit, ARB]	# cholesky factorization
+PIXEL	cov[ARB]		# covariance vector
+PIXEL	errors[ARB]		# computed errors
+int	list[ARB]		# list of active parameters
+int	nfit			# number of fitted parameters
+PIXEL	variance		# variance of the fit
+
+int	i, ier
+
+begin
+	# Factorize the data matrix to determine the errors.
+	call nl_chfac$t (alpha, nfit, nfit, chofac, ier)
+
+	# Estimate the errors.
+	do i = 1, nfit {
+	    call aclr$t (cov, nfit)
+	    cov[i] = PIXEL (1.0)
+	    call nl_chslv$t (chofac, nfit, nfit, cov, cov)
+	    if (cov[i] >= PIXEL (0.0))
+	        errors[list[i]] = sqrt (cov[i] * variance)
+	    else
+		errors[list[i]] = PIXEL (0.0)
+	}
+end
diff --git a/math/nlfit/nlerrorsd.x b/math/nlfit/nlerrorsd.x
new file mode 100644
index 00000000..b23a7b36
--- /dev/null
+++ b/math/nlfit/nlerrorsd.x
@@ -0,0 +1,107 @@
+include <mach.h>
+include <math/nlfit.h>
+include "nlfitdefd.h"
+
+define	COV	Memd[P2P($1)]	# element of COV
+
+
+# NLERRORS -- Procedure to calculate the reduced chi-squared of the fit
+# and the errors of the coefficients. First the variance
+# and the reduced chi-squared of the fit are estimated. If these two
+# quantities are identical the variance is used to scale the errors
+# in the coefficients. The errors in the coefficients are proportional
+# to the inverse diagonal elements of MATRIX.
+
+procedure nlerrorsd (nl, z, zfit, w, npts, variance, chisqr, errors)
+
+pointer	nl		# curve descriptor
+double	z[ARB]		# data points
+double	zfit[ARB]	# fitted data points
+double	w[ARB]		# array of weights
+int	npts		# number of points
+double	variance	# variance of the fit
+double	chisqr		# reduced chi-squared of fit (output)
+double	errors[ARB]	# errors in coefficients (output)
+
+int	i, n, nfree
+pointer	sp, covptr
+double   factor
+
+begin
+	# Allocate space for covariance vector.
+	call smark (sp)
+	call salloc (covptr, NL_NPARAMS(nl), TY_DOUBLE)
+
+	# Estimate the variance and reduce chi-squared of the fit.
+	n = 0
+	variance = double (0.0)
+	chisqr = double (0.0)
+	do i = 1, npts {
+	    if (w[i] <= double (0.0))
+	        next
+	    factor = (z[i] - zfit[i]) ** 2
+	    variance = variance + factor
+	    chisqr = chisqr + factor * w[i]
+	    n = n + 1
+	}
+
+	# Calculate the reduced chi-squared.
+	nfree = n - NL_NFPARAMS(nl)
+	if (nfree  > 0) {
+	    variance = variance / nfree
+	    chisqr = chisqr / nfree
+	} else {
+	    variance = double (0.0)
+	    chisqr = double (0.0)
+	}
+
+	# If the variance equals the reduced chi_squared as in the case
+	# of uniform weights then scale the errors in the coefficients
+	# by the variance and not the reduced chi-squared
+
+	if (abs (chisqr - variance) <= EPSILOND) {
+	    if (nfree > 0)
+		factor = chisqr
+	    else
+		factor = double (0.0)
+	} else
+	    factor = 1.
+
+	# Calculate the  errors in the coefficients.
+	call aclrd (errors, NL_NPARAMS(nl))
+	call nlinvertd (ALPHA(NL_ALPHA(nl)), CHOFAC(NL_CHOFAC(nl)),
+	    COV(covptr), errors, PLIST(NL_PLIST(nl)), NL_NFPARAMS(nl), factor)
+
+	call sfree (sp)
+end
+
+
+# NLINVERT -- Procedure to invert matrix and compute errors
+
+procedure nlinvertd (alpha, chofac, cov, errors, list, nfit, variance)
+
+double	alpha[nfit,ARB]		# data matrix
+double	chofac[nfit, ARB]	# cholesky factorization
+double	cov[ARB]		# covariance vector
+double	errors[ARB]		# computed errors
+int	list[ARB]		# list of active parameters
+int	nfit			# number of fitted parameters
+double	variance		# variance of the fit
+
+int	i, ier
+
+begin
+	# Factorize the data matrix to determine the errors.
+	call nl_chfacd (alpha, nfit, nfit, chofac, ier)
+
+	# Estimate the errors.
+	do i = 1, nfit {
+	    call aclrd (cov, nfit)
+	    cov[i] = double (1.0)
+	    call nl_chslvd (chofac, nfit, nfit, cov, cov)
+	    if (cov[i] >= double (0.0))
+	        errors[list[i]] = sqrt (cov[i] * variance)
+	    else
+		errors[list[i]] = double (0.0)
+	}
+end
diff --git a/math/nlfit/nlerrorsr.x b/math/nlfit/nlerrorsr.x
new file mode 100644
index 00000000..40057b9c
--- /dev/null
+++ b/math/nlfit/nlerrorsr.x
@@ -0,0 +1,107 @@
+include <mach.h>
+include <math/nlfit.h>
+include "nlfitdefr.h"
+
+define	COV	Memr[P2P($1)]	# element of COV
+
+
+# NLERRORS -- Procedure to calculate the reduced chi-squared of the fit
+# and the errors of the coefficients. First the variance
+# and the reduced chi-squared of the fit are estimated. If these two
+# quantities are identical the variance is used to scale the errors
+# in the coefficients. The errors in the coefficients are proportional
+# to the inverse diagonal elements of MATRIX.
+
+procedure nlerrorsr (nl, z, zfit, w, npts, variance, chisqr, errors)
+
+pointer	nl		# curve descriptor
+real	z[ARB]		# data points
+real	zfit[ARB]	# fitted data points
+real	w[ARB]		# array of weights
+int	npts		# number of points
+real	variance	# variance of the fit
+real	chisqr		# reduced chi-squared of fit (output)
+real	errors[ARB]	# errors in coefficients (output)
+
+int	i, n, nfree
+pointer	sp, covptr
+real   factor
+
+begin
+	# Allocate space for covariance vector.
+	call smark (sp)
+	call salloc (covptr, NL_NPARAMS(nl), TY_REAL)
+
+	# Estimate the variance and reduce chi-squared of the fit.
+	n = 0
+	variance = real (0.0)
+	chisqr = real (0.0)
+	do i = 1, npts {
+	    if (w[i] <= real (0.0))
+	        next
+	    factor = (z[i] - zfit[i]) ** 2
+	    variance = variance + factor
+	    chisqr = chisqr + factor * w[i]
+	    n = n + 1
+	}
+
+	# Calculate the reduced chi-squared.
+	nfree = n - NL_NFPARAMS(nl)
+	if (nfree  > 0) {
+	    variance = variance / nfree
+	    chisqr = chisqr / nfree
+	} else {
+	    variance = real (0.0)
+	    chisqr = real (0.0)
+	}
+
+	# If the variance equals the reduced chi_squared as in the case
+	# of uniform weights then scale the errors in the coefficients
+	# by the variance and not the reduced chi-squared
+
+	if (abs (chisqr - variance) <= EPSILONR) {
+	    if (nfree > 0)
+		factor = chisqr
+	    else
+		factor = real (0.0)
+	} else
+	    factor = 1.
+
+	# Calculate the  errors in the coefficients.
+	call aclrr (errors, NL_NPARAMS(nl))
+	call nlinvertr (ALPHA(NL_ALPHA(nl)), CHOFAC(NL_CHOFAC(nl)),
+	    COV(covptr), errors, PLIST(NL_PLIST(nl)), NL_NFPARAMS(nl), factor)
+
+	call sfree (sp)
+end
+
+
+# NLINVERT -- Procedure to invert matrix and compute errors
+
+procedure nlinvertr (alpha, chofac, cov, errors, list, nfit, variance)
+
+real	alpha[nfit,ARB]		# data matrix
+real	chofac[nfit, ARB]	# cholesky factorization
+real	cov[ARB]		# covariance vector
+real	errors[ARB]		# computed errors
+int	list[ARB]		# list of active parameters
+int	nfit			# number of fitted parameters
+real	variance		# variance of the fit
+
+int	i, ier
+
+begin
+	# Factorize the data matrix to determine the errors.
+	call nl_chfacr (alpha, nfit, nfit, chofac, ier)
+
+	# Estimate the errors.
+	do i = 1, nfit {
+	    call aclrr (cov, nfit)
+	    cov[i] = real (1.0)
+	    call nl_chslvr (chofac, nfit, nfit, cov, cov)
+	    if (cov[i] >= real (0.0))
+	        errors[list[i]] = sqrt (cov[i] * variance)
+	    else
+		errors[list[i]] = real (0.0)
+	}
+end
diff --git a/math/nlfit/nleval.gx b/math/nlfit/nleval.gx
new file mode 100644
index 00000000..e116982b
--- /dev/null
+++ b/math/nlfit/nleval.gx
@@ -0,0 +1,25 @@
+include <math/nlfit.h>
+$if (datatype == r)
+include "nlfitdefr.h"
+$else
+include "nlfitdefd.h"
+$endif
+
+
+# NLEVAL -- Evaluate the fit at a point.
+
+PIXEL procedure nleval$t (nl, x, nvars)
+
+pointer	nl		# nlfit descriptor
+PIXEL	x[ARB]		# x values
+int	nvars		# number of variables
+
+real	zfit
+
+begin
+	# Evaluate the function.
+	call zcall5 (NL_FUNC(nl), x, nvars, PARAM(NL_PARAM(nl)),
+	    NL_NPARAMS(nl), zfit)
+
+	return (zfit)
+end
diff --git a/math/nlfit/nlevald.x b/math/nlfit/nlevald.x
new file mode 100644
index 00000000..c7710acc
--- /dev/null
+++ b/math/nlfit/nlevald.x
@@ -0,0 +1,21 @@
+include <math/nlfit.h>
+include "nlfitdefd.h"
+
+
+# NLEVAL -- Evaluate the fit at a point.
+
+double procedure nlevald (nl, x, nvars)
+
+pointer	nl		# nlfit descriptor
+double	x[ARB]		# x values
+int	nvars		# number of variables
+
+real	zfit
+
+begin
+	# Evaluate the function.
+	call zcall5 (NL_FUNC(nl), x, nvars, PARAM(NL_PARAM(nl)),
+	    NL_NPARAMS(nl), zfit)
+
+	return (zfit)
+end
diff --git a/math/nlfit/nlevalr.x b/math/nlfit/nlevalr.x
new file mode 100644
index 00000000..4c68d865
--- /dev/null
+++ b/math/nlfit/nlevalr.x
@@ -0,0 +1,21 @@
+include <math/nlfit.h>
+include "nlfitdefr.h"
+
+
+# NLEVAL -- Evaluate the fit at a point.
+
+real procedure nlevalr (nl, x, nvars)
+
+pointer	nl		# nlfit descriptor
+real	x[ARB]		# x values
+int	nvars		# number of variables
+
+real	zfit
+
+begin
+	# Evaluate the function.
+	call zcall5 (NL_FUNC(nl), x, nvars, PARAM(NL_PARAM(nl)),
+	    NL_NPARAMS(nl), zfit)
+
+	return (zfit)
+end
diff --git a/math/nlfit/nlfit.gx b/math/nlfit/nlfit.gx
new file mode 100644
index 00000000..50dc0b21
--- /dev/null
+++ b/math/nlfit/nlfit.gx
@@ -0,0 +1,171 @@
+include <mach.h>
+include <math/nlfit.h>
+$if (datatype == r)
+include "nlfitdefr.h"
+$else
+include "nlfitdefd.h"
+$endif
+
+# NLFIT -- Routine to perform a non-linear least squares fit. At least MINITER
+# iterations must be performed before the fit is complete.
+
+procedure nlfit$t (nl, x, z, w, npts, nvars, wtflag, stat)
+
+pointer	nl		# pointer to nlfit structure
+PIXEL	x[ARB]		# independent variables (npts * nvars)
+PIXEL	z[ARB]		# function values (npts)
+PIXEL	w[ARB]		# weights (npts)
+int	npts		# number of points
+int	nvars		# number of independent variables
+int	wtflag		# weighting type
+int	stat		# error code
+
+int	i, miniter, ier
+PIXEL	scatter, dscatter
+PIXEL	nlscatter$t()
+
+begin
+	# Initialize.
+	NL_ITER(nl) = 0
+	NL_LAMBDA(nl) = PIXEL (.001)
+	NL_REFSQ(nl) = PIXEL (0.0)
+	NL_SCATTER(nl) = PIXEL(0.0)
+
+	# Initialize the weights.
+	switch (wtflag) {
+	case WTS_UNIFORM:
+	    do i = 1, npts
+		w[i] = PIXEL (1.0)
+	case WTS_SCATTER:
+	    ;
+	case WTS_USER:
+	    ;
+	case WTS_CHISQ:
+	    do i = 1, npts {
+		if (z[i] > PIXEL (0.0))
+		    w[i] = PIXEL (1.0) / z[i]
+		else if (z[i] < PIXEL (0.0))
+		    w[i] = PIXEL (-1.0) / z[i]
+	        else
+	            w[i] = PIXEL (0.0)
+	    }
+	default:
+	    do i = 1, npts
+		w[i] = PIXEL (1.0)
+	}
+
+	# Initialize.
+	scatter = PIXEL(0.0)
+	if (wtflag == WTS_SCATTER)
+	    miniter = MINITER + 1
+	else
+	    miniter = MINITER
+
+	repeat {
+
+	    # Perform a single iteration.
+	    call nliter$t (nl, x, z, w, npts, nvars, ier)
+	    NL_ITER(nl) = NL_ITER(nl) + 1
+	    #call eprintf ("niter=%d refsq=%g oldsq=%g sumsq=%g\n")
+		#call pargi (NL_ITER(nl))
+		#call parg$t (NL_REFSQ(nl))
+		#call parg$t (NL_OLDSQ(nl))
+		#call parg$t (NL_SUMSQ(nl))
+	    stat = ier
+	    if (stat == NO_DEG_FREEDOM)
+	        break
+
+	    # Make the convergence checks.
+	    if (NL_ITER(nl) < miniter)
+		stat = NOT_DONE
+	    else if (NL_SUMSQ(nl) <= PIXEL (10.0) * EPSILON$T)
+		stat = DONE
+	    else if (NL_REFSQ(nl) < NL_SUMSQ(nl))
+		stat = NOT_DONE
+	    else if (((NL_REFSQ(nl) - NL_SUMSQ(nl)) / NL_SUMSQ(nl)) <
+	        NL_TOL(nl))
+		stat = DONE
+	    else
+		stat = NOT_DONE
+
+	    # Check for a singular solution.
+	    if (stat == DONE) {
+		if (ier == SINGULAR)
+		    stat = ier
+		break
+	    }
+
+	    # Quit if the lambda parameter goes to zero.
+	    if (NL_LAMBDA(nl) <= 0.0)
+		break
+
+	    # Check the number of iterations.
+	    if ((NL_ITER(nl) >= miniter) && (NL_ITER(nl) >= NL_ITMAX(nl)))
+		break
+
+	    # Adjust the weights if necessary.
+	    switch (wtflag) {
+	    case WTS_SCATTER:
+	        dscatter = nlscatter$t (nl, x, z, w, npts, nvars)
+		if ((NL_ITER(nl) >= MINITER) && (dscatter >= PIXEL (10.0) *
+		    EPSILON$T)) {
+		    do i = 1, npts {
+			if (w[i] <= PIXEL(0.0))
+			    w[i] = PIXEL(0.0)
+			else {
+		            w[i] = PIXEL(1.0) / (PIXEL(1.0) / w[i] + dscatter)
+			}
+		    }
+		    scatter = scatter + dscatter
+		}
+	    default:
+		;
+	    }
+
+	    # Get ready for next iteration.
+	    NL_REFSQ(nl) = min (NL_OLDSQ(nl), NL_SUMSQ(nl))
+	}
+
+	NL_SCATTER(nl) = NL_SCATTER(nl) + scatter
+end
+
+
+# NLSCATTER -- Routine to estimate the original scatter in the fit.
+
+PIXEL procedure nlscatter$t (nl, x, z, w, npts, nvars)
+
+pointer	nl		# Pointer to nl fitting structure
+PIXEL	x[ARB]		# independent variables (npts * nvars)
+PIXEL	z[ARB]		# function values (npts)
+PIXEL	w[ARB]		# weights (npts)
+int	npts		# number of points
+int	nvars		# number of independent variables
+
+pointer	sp, zfit, errors
+PIXEL	scatter, variance, chisqr
+int	nlstati()
+
+begin
+	# Allocate working memory.
+	call smark (sp)
+	call salloc (zfit, npts, TY_PIXEL)
+	call salloc (errors, nlstati (nl, NLNPARAMS), TY_PIXEL)
+
+	# Initialize
+	scatter = PIXEL (0.0)
+
+	# Compute the fit and the errors.
+	call nlvector$t (nl, x, Mem$t[zfit], npts, nvars)
+	call nlerrors$t (nl, z, Mem$t[zfit], w, npts, variance, chisqr,
+	    Mem$t[errors])
+
+	# Estimate the scatter.
+	if (chisqr <= PIXEL(0.0) || variance <= PIXEL(0.0))
+	    scatter = PIXEL (0.0)
+	else
+	    scatter = PIXEL(0.5) * variance * (chisqr - PIXEL(1.0)) / chisqr
+
+	call sfree (sp)
+
+	return (scatter)
+end
diff --git a/math/nlfit/nlfitd.x b/math/nlfit/nlfitd.x
new file mode 100644
index 00000000..975284ab
--- /dev/null
+++ b/math/nlfit/nlfitd.x
@@ -0,0 +1,167 @@
+include <mach.h>
+include <math/nlfit.h>
+include "nlfitdefd.h"
+
+# NLFIT -- Routine to perform a non-linear least squares fit. At least MINITER
+# iterations must be performed before the fit is complete.
+
+procedure nlfitd (nl, x, z, w, npts, nvars, wtflag, stat)
+
+pointer	nl		# pointer to nlfit structure
+double	x[ARB]		# independent variables (npts * nvars)
+double	z[ARB]		# function values (npts)
+double	w[ARB]		# weights (npts)
+int	npts		# number of points
+int	nvars		# number of independent variables
+int	wtflag		# weighting type
+int	stat		# error code
+
+int	i, miniter, ier
+double	scatter, dscatter
+double	nlscatterd()
+
+begin
+	# Initialize.
+	NL_ITER(nl) = 0
+	NL_LAMBDA(nl) = double (.001)
+	NL_REFSQ(nl) = double (0.0)
+	NL_SCATTER(nl) = double(0.0)
+
+	# Initialize the weights.
+	switch (wtflag) {
+	case WTS_UNIFORM:
+	    do i = 1, npts
+		w[i] = double (1.0)
+	case WTS_SCATTER:
+	    ;
+	case WTS_USER:
+	    ;
+	case WTS_CHISQ:
+	    do i = 1, npts {
+		if (z[i] > double (0.0))
+		    w[i] = double (1.0) / z[i]
+		else if (z[i] < double (0.0))
+		    w[i] = double (-1.0) / z[i]
+	        else
+	            w[i] = double (0.0)
+	    }
+	default:
+	    do i = 1, npts
+		w[i] = double (1.0)
+	}
+
+	# Initialize.
+	scatter = double(0.0)
+	if (wtflag == WTS_SCATTER)
+	    miniter = MINITER + 1
+	else
+	    miniter = MINITER
+
+	repeat {
+
+	    # Perform a single iteration.
+	    call nliterd (nl, x, z, w, npts, nvars, ier)
+	    NL_ITER(nl) = NL_ITER(nl) + 1
+	    #call eprintf ("niter=%d refsq=%g oldsq=%g sumsq=%g\n")
+		#call pargi (NL_ITER(nl))
+		#call parg$t (NL_REFSQ(nl))
+		#call parg$t (NL_OLDSQ(nl))
+		#call parg$t (NL_SUMSQ(nl))
+	    stat = ier
+	    if (stat == NO_DEG_FREEDOM)
+	        break
+
+	    # Make the convergence checks.
+	    if (NL_ITER(nl) < miniter)
+		stat = NOT_DONE
+	    else if (NL_SUMSQ(nl) <= double (10.0) * EPSILOND)
+		stat = DONE
+	    else if (NL_REFSQ(nl) < NL_SUMSQ(nl))
+		stat = NOT_DONE
+	    else if (((NL_REFSQ(nl) - NL_SUMSQ(nl)) / NL_SUMSQ(nl)) <
+	        NL_TOL(nl))
+		stat = DONE
+	    else
+		stat = NOT_DONE
+
+	    # Check for a singular solution.
+	    if (stat == DONE) {
+		if (ier == SINGULAR)
+		    stat = ier
+		break
+	    }
+
+	    # Quit if the lambda parameter goes to zero.
+	    if (NL_LAMBDA(nl) <= 0.0)
+		break
+
+	    # Check the number of iterations.
+	    if ((NL_ITER(nl) >= miniter) && (NL_ITER(nl) >= NL_ITMAX(nl)))
+		break
+
+	    # Adjust the weights if necessary.
+	    switch (wtflag) {
+	    case WTS_SCATTER:
+	        dscatter = nlscatterd (nl, x, z, w, npts, nvars)
+		if ((NL_ITER(nl) >= MINITER) && (dscatter >= double (10.0) *
+		    EPSILOND)) {
+		    do i = 1, npts {
+			if (w[i] <= double(0.0))
+			    w[i] = double(0.0)
+			else {
+		            w[i] = double(1.0) / (double(1.0) / w[i] + dscatter)
+			}
+		    }
+		    scatter = scatter + dscatter
+		}
+	    default:
+		;
+	    }
+
+	    # Get ready for next iteration.
+	    NL_REFSQ(nl) = min (NL_OLDSQ(nl), NL_SUMSQ(nl))
+	}
+
+	NL_SCATTER(nl) = NL_SCATTER(nl) + scatter
+end
+
+
+# NLSCATTER -- Routine to estimate the original scatter in the fit.
+
+double procedure nlscatterd (nl, x, z, w, npts, nvars)
+
+pointer	nl		# Pointer to nl fitting structure
+double	x[ARB]		# independent variables (npts * nvars)
+double	z[ARB]		# function values (npts)
+double	w[ARB]		# weights (npts)
+int	npts		# number of points
+int	nvars		# number of independent variables
+
+pointer	sp, zfit, errors
+double	scatter, variance, chisqr
+int	nlstati()
+
+begin
+	# Allocate working memory.
+	call smark (sp)
+	call salloc (zfit, npts, TY_DOUBLE)
+	call salloc (errors, nlstati (nl, NLNPARAMS), TY_DOUBLE)
+
+	# Initialize
+	scatter = double (0.0)
+
+	# Compute the fit and the errors.
+	call nlvectord (nl, x, Memd[zfit], npts, nvars)
+	call nlerrorsd (nl, z, Memd[zfit], w, npts, variance, chisqr,
+	    Memd[errors])
+
+	# Estimate the scatter.
+	if (chisqr <= double(0.0) || variance <= double(0.0))
+	    scatter = double (0.0)
+	else
+	    scatter = double(0.5) * variance * (chisqr - double(1.0)) / chisqr
+
+	call sfree (sp)
+
+	return (scatter)
+end
diff --git a/math/nlfit/nlfitdef.gh b/math/nlfit/nlfitdef.gh
new file mode 100644
index 00000000..62a542a8
--- /dev/null
+++ b/math/nlfit/nlfitdef.gh
@@ -0,0 +1,51 @@
+# The NLFIT data structure.
+
+# Structure length
+define	LEN_NLSTRUCT	35
+
+
+# Structure definition
+define  NL_TOL          Mem$t[P2$T($1+0)]   # Tolerance for convergence
+define  NL_LAMBDA       Mem$t[P2$T($1+2)]   # Damping factor
+define	NL_OLDSQ	Mem$t[P2$T($1+4)]   # Sum resid. squared last iter.
+define	NL_SUMSQ	Mem$t[P2$T($1+6)]   # Sum of residuals squared 
+define	NL_REFSQ	Mem$t[P2$T($1+8)]   # Reference sum of squares
+define	NL_SCATTER	Mem$t[P2$T($1+10)]  # Additional scatter
+
+define	NL_FUNC		Memi[$1+12]	# Fitting function
+define	NL_DFUNC	Memi[$1+13]	# Derivative function
+define	NL_NPARAMS	Memi[$1+14]	# Number of parameters
+define	NL_NFPARAMS	Memi[$1+15]	# Number of fitted parameters
+define	NL_ITMAX	Memi[$1+16]	# Max number of iterations
+define	NL_ITER		Memi[$1+17]	# Iteration counter
+define	NL_NPTS		Memi[$1+18]	# Number of points in fit
+
+define	NL_PARAM	Memi[$1+19]	# Pointer to parameter vector
+define	NL_OPARAM	Memi[$1+20]	# Pointer to orignal parameter vector
+define	NL_DPARAM	Memi[$1+21]	# Pointer to parameter change vector
+define	NL_DELPARAM	Memi[$1+22]	# Pointer to delta param vector
+define	NL_PLIST	Memi[$1+23]	# Pointer to parameter list
+define	NL_ALPHA	Memi[$1+24]	# Pointer to alpha matrix
+define	NL_COVAR	Memi[$1+25]	# Pointer to covariance matrix
+define	NL_BETA		Memi[$1+26]	# Pointer to beta matrix
+define	NL_TRY		Memi[$1+27]	# Pointer to trial vector
+define	NL_DERIV	Memi[$1+28]	# Pointer to derivatives
+define	NL_CHOFAC	Memi[$1+29]	# Pointer to Cholesky factorization
+
+# next free location	($1 + 30)
+
+# Access to buffers
+define	PLIST		Memi[$1]	# Parameter list
+define	OPARAM		Mem$t[P2P($1)]	# Original parameter vector
+define	PARAM		Mem$t[P2P($1)]	# Parameter vector
+define	DPARAM		Mem$t[P2P($1)]	# Parameter change vector
+define	ALPHA		Mem$t[P2P($1)]	# Alpha matrix
+define	BETA		Mem$t[P2P($1)]	# Beta matrix
+define	TRY		Mem$t[P2P($1)]	# Trial vector
+define	DERIV		Mem$t[P2P($1)]	# Derivatives
+define	CHOFAC		Mem$t[P2P($1)]	# Cholesky factorization
+define	COVAR		Mem$t[P2P($1)]	# Covariance matrix
+
+# Defined constants alter for tricky problems
+define	LAMBDAMAX	1000.0
+define	MINITER		3
diff --git a/math/nlfit/nlfitdefd.h b/math/nlfit/nlfitdefd.h
new file mode 100644
index 00000000..b5ed3bf3
--- /dev/null
+++ b/math/nlfit/nlfitdefd.h
@@ -0,0 +1,52 @@
+# The NLFIT data structure.
+
+# Structure length
+define	LEN_NLSTRUCT	35
+
+
+# Structure definition
+define  NL_TOL          Memd[P2D($1+0)]   # Tolerance for convergence
+define  NL_LAMBDA       Memd[P2D($1+2)]   # Damping factor
+define	NL_OLDSQ	Memd[P2D($1+4)]   # Sum resid. squared last iter.
+define	NL_SUMSQ	Memd[P2D($1+6)]   # Sum of residuals squared 
+define	NL_REFSQ	Memd[P2D($1+8)]   # Reference sum of squares
+define	NL_SCATTER	Memd[P2D($1+10)]  # Additional scatter
+
+define	NL_FUNC		Memi[$1+12]	# Fitting function
+define	NL_DFUNC	Memi[$1+13]	# Derivative function
+define	NL_NPARAMS	Memi[$1+14]	# Number of parameters
+define	NL_NFPARAMS	Memi[$1+15]	# Number of fitted parameters
+define	NL_ITMAX	Memi[$1+16]	# Max number of iterations
+define	NL_ITER		Memi[$1+17]	# Iteration counter
+define	NL_NPTS		Memi[$1+18]	# Number of points in fit
+
+define	NL_PARAM	Memi[$1+19]	# Pointer to parameter vector
+define	NL_OPARAM	Memi[$1+20]	# Pointer to orignal parameter vector
+define	NL_DPARAM	Memi[$1+21]	# Pointer to parameter change vector
+define	NL_DELPARAM	Memi[$1+22]	# Pointer to delta param vector
+define	NL_PLIST	Memi[$1+23]	# Pointer to parameter list
+define	NL_ALPHA	Memi[$1+24]	# Pointer to alpha matrix
+define	NL_COVAR	Memi[$1+25]	# Pointer to covariance matrix
+define	NL_BETA		Memi[$1+26]	# Pointer to beta matrix
+define	NL_TRY		Memi[$1+27]	# Pointer to trial vector
+define	NL_DERIV	Memi[$1+28]	# Pointer to derivatives
+define	NL_CHOFAC	Memi[$1+29]	# Pointer to Cholesky factorization
+
+# next free location	($1 + 30)
+
+# Access to buffers
+define	PLIST		Memi[$1]	# Parameter list
+define	OPARAM		Memd[$1]	# Original parameter vector
+define	PARAM		Memd[$1]	# Parameter vector
+define	DPARAM		Memd[$1]	# Parameter change vector
+define	ALPHA		Memd[$1]	# Alpha matrix
+define	BETA		Memd[$1]	# Beta matrix
+define	TRY		Memd[$1]	# Trial vector
+define	DERIV		Memd[$1]	# Derivatives
+define	CHOFAC		Memd[$1]	# Cholesky factorization
+define	COVAR		Memd[$1]	# Covariance matrix
+
+
+# Defined constants alter for tricky problems
+define	LAMBDAMAX	1000.0
+define	MINITER		3
diff --git a/math/nlfit/nlfitdefr.h b/math/nlfit/nlfitdefr.h
new file mode 100644
index 00000000..e3bf3ae0
--- /dev/null
+++ b/math/nlfit/nlfitdefr.h
@@ -0,0 +1,52 @@
+# The NLFIT data structure.
+
+# Structure length
+define	LEN_NLSTRUCT	35
+
+
+# Structure definition
+define  NL_TOL          Memr[P2R($1+0)]   # Tolerance for convergence
+define  NL_LAMBDA       Memr[P2R($1+2)]   # Damping factor
+define	NL_OLDSQ	Memr[P2R($1+4)]   # Sum resid. squared last iter.
+define	NL_SUMSQ	Memr[P2R($1+6)]   # Sum of residuals squared 
+define	NL_REFSQ	Memr[P2R($1+8)]   # Reference sum of squares
+define	NL_SCATTER	Memr[P2R($1+10)]  # Additional scatter
+
+define	NL_FUNC		Memi[$1+12]	# Fitting function
+define	NL_DFUNC	Memi[$1+13]	# Derivative function
+define	NL_NPARAMS	Memi[$1+14]	# Number of parameters
+define	NL_NFPARAMS	Memi[$1+15]	# Number of fitted parameters
+define	NL_ITMAX	Memi[$1+16]	# Max number of iterations
+define	NL_ITER		Memi[$1+17]	# Iteration counter
+define	NL_NPTS		Memi[$1+18]	# Number of points in fit
+
+define	NL_PARAM	Memi[$1+19]	# Pointer to parameter vector
+define	NL_OPARAM	Memi[$1+20]	# Pointer to orignal parameter vector
+define	NL_DPARAM	Memi[$1+21]	# Pointer to parameter change vector
+define	NL_DELPARAM	Memi[$1+22]	# Pointer to delta param vector
+define	NL_PLIST	Memi[$1+23]	# Pointer to parameter list
+define	NL_ALPHA	Memi[$1+24]	# Pointer to alpha matrix
+define	NL_COVAR	Memi[$1+25]	# Pointer to covariance matrix
+define	NL_BETA		Memi[$1+26]	# Pointer to beta matrix
+define	NL_TRY		Memi[$1+27]	# Pointer to trial vector
+define	NL_DERIV	Memi[$1+28]	# Pointer to derivatives
+define	NL_CHOFAC	Memi[$1+29]	# Pointer to Cholesky factorization
+
+# next free location	($1 + 30)
+
+# Access to buffers
+define	PLIST		Memi[$1]	# Parameter list
+define	OPARAM		Memr[$1]	# Original parameter vector
+define	PARAM		Memr[$1]	# Parameter vector
+define	DPARAM		Memr[$1]	# Parameter change vector
+define	ALPHA		Memr[$1]	# Alpha matrix
+define	BETA		Memr[$1]	# Beta matrix
+define	TRY		Memr[$1]	# Trial vector
+define	DERIV		Memr[$1]	# Derivatives
+define	CHOFAC		Memr[$1]	# Cholesky factorization
+define	COVAR		Memr[$1]	# Covariance matrix
+
+
+# Defined constants alter for tricky problems
+define	LAMBDAMAX	1000.0
+define	MINITER		3
diff --git a/math/nlfit/nlfitr.x b/math/nlfit/nlfitr.x
new file mode 100644
index 00000000..dc06b720
--- /dev/null
+++ b/math/nlfit/nlfitr.x
@@ -0,0 +1,167 @@
+include <mach.h>
+include <math/nlfit.h>
+include "nlfitdefr.h"
+
+# NLFIT -- Routine to perform a non-linear least squares fit. At least MINITER
+# iterations must be performed before the fit is complete.
+
+procedure nlfitr (nl, x, z, w, npts, nvars, wtflag, stat)
+
+pointer	nl		# pointer to nlfit structure
+real	x[ARB]		# independent variables (npts * nvars)
+real	z[ARB]		# function values (npts)
+real	w[ARB]		# weights (npts)
+int	npts		# number of points
+int	nvars		# number of independent variables
+int	wtflag		# weighting type
+int	stat		# error code
+
+int	i, miniter, ier
+real	scatter, dscatter
+real	nlscatterr()
+
+begin
+	# Initialize.
+	NL_ITER(nl) = 0
+	NL_LAMBDA(nl) = real (.001)
+	NL_REFSQ(nl) = real (0.0)
+	NL_SCATTER(nl) = real(0.0)
+
+	# Initialize the weights.
+	switch (wtflag) {
+	case WTS_UNIFORM:
+	    do i = 1, npts
+		w[i] = real (1.0)
+	case WTS_SCATTER:
+	    ;
+	case WTS_USER:
+	    ;
+	case WTS_CHISQ:
+	    do i = 1, npts {
+		if (z[i] > real (0.0))
+		    w[i] = real (1.0) / z[i]
+		else if (z[i] < real (0.0))
+		    w[i] = real (-1.0) / z[i]
+	        else
+	            w[i] = real (0.0)
+	    }
+	default:
+	    do i = 1, npts
+		w[i] = real (1.0)
+	}
+
+	# Initialize.
+	scatter = real(0.0)
+	if (wtflag == WTS_SCATTER)
+	    miniter = MINITER + 1
+	else
+	    miniter = MINITER
+
+	repeat {
+
+	    # Perform a single iteration.
+	    call nliterr (nl, x, z, w, npts, nvars, ier)
+	    NL_ITER(nl) = NL_ITER(nl) + 1
+	    #call eprintf ("niter=%d refsq=%g oldsq=%g sumsq=%g\n")
+		#call pargi (NL_ITER(nl))
+		#call parg$t (NL_REFSQ(nl))
+		#call parg$t (NL_OLDSQ(nl))
+		#call parg$t (NL_SUMSQ(nl))
+	    stat = ier
+	    if (stat == NO_DEG_FREEDOM)
+	        break
+
+	    # Make the convergence checks.
+	    if (NL_ITER(nl) < miniter)
+		stat = NOT_DONE
+	    else if (NL_SUMSQ(nl) <= real (10.0) * EPSILONR)
+		stat = DONE
+	    else if (NL_REFSQ(nl) < NL_SUMSQ(nl))
+		stat = NOT_DONE
+	    else if (((NL_REFSQ(nl) - NL_SUMSQ(nl)) / NL_SUMSQ(nl)) <
+	        NL_TOL(nl))
+		stat = DONE
+	    else
+		stat = NOT_DONE
+
+	    # Check for a singular solution.
+	    if (stat == DONE) {
+		if (ier == SINGULAR)
+		    stat = ier
+		break
+	    }
+
+	    # Quit if the lambda parameter goes to zero.
+	    if (NL_LAMBDA(nl) <= 0.0)
+		break
+
+	    # Check the number of iterations.
+	    if ((NL_ITER(nl) >= miniter) && (NL_ITER(nl) >= NL_ITMAX(nl)))
+		break
+
+	    # Adjust the weights if necessary.
+	    switch (wtflag) {
+	    case WTS_SCATTER:
+	        dscatter = nlscatterr (nl, x, z, w, npts, nvars)
+		if ((NL_ITER(nl) >= MINITER) && (dscatter >= real (10.0) *
+		    EPSILONR)) {
+		    do i = 1, npts {
+			if (w[i] <= real(0.0))
+			    w[i] = real(0.0)
+			else {
+		            w[i] = real(1.0) / (real(1.0) / w[i] + dscatter)
+			}
+		    }
+		    scatter = scatter + dscatter
+		}
+	    default:
+		;
+	    }
+
+	    # Get ready for next iteration.
+	    NL_REFSQ(nl) = min (NL_OLDSQ(nl), NL_SUMSQ(nl))
+	}
+
+	NL_SCATTER(nl) = NL_SCATTER(nl) + scatter
+end
+
+
+# NLSCATTER -- Routine to estimate the original scatter in the fit.
+
+real procedure nlscatterr (nl, x, z, w, npts, nvars)
+
+pointer	nl		# Pointer to nl fitting structure
+real	x[ARB]		# independent variables (npts * nvars)
+real	z[ARB]		# function values (npts)
+real	w[ARB]		# weights (npts)
+int	npts		# number of points
+int	nvars		# number of independent variables
+
+pointer	sp, zfit, errors
+real	scatter, variance, chisqr
+int	nlstati()
+
+begin
+	# Allocate working memory.
+	call smark (sp)
+	call salloc (zfit, npts, TY_REAL)
+	call salloc (errors, nlstati (nl, NLNPARAMS), TY_REAL)
+
+	# Initialize
+	scatter = real (0.0)
+
+	# Compute the fit and the errors.
+	call nlvectorr (nl, x, Memr[zfit], npts, nvars)
+	call nlerrorsr (nl, z, Memr[zfit], w, npts, variance, chisqr,
+	    Memr[errors])
+
+	# Estimate the scatter.
+	if (chisqr <= real(0.0) || variance <= real(0.0))
+	    scatter = real (0.0)
+	else
+	    scatter = real(0.5) * variance * (chisqr - real(1.0)) / chisqr
+
+	call sfree (sp)
+
+	return (scatter)
+end
diff --git a/math/nlfit/nlfree.gx b/math/nlfit/nlfree.gx
new file mode 100644
index 00000000..507dce57
--- /dev/null
+++ b/math/nlfit/nlfree.gx
@@ -0,0 +1,41 @@
+$if (datatype == r)
+include "nlfitdefr.h"
+$else
+include "nlfitdefd.h"
+$endif
+
+# NLFREE -- Deallocate all assigned space
+
+procedure nlfree$t (nl)
+
+pointer	nl		# pointer to non-linear fitting structure
+
+errchk	mfree
+
+begin
+	if (nl == NULL)
+	    return
+	if (NL_PARAM(nl) != NULL)
+	    call mfree (NL_PARAM(nl), TY_PIXEL)
+	if (NL_OPARAM(nl) != NULL)
+	    call mfree (NL_OPARAM(nl), TY_PIXEL)
+	if (NL_DPARAM(nl) != NULL)
+	    call mfree (NL_DPARAM(nl), TY_PIXEL)
+	if (NL_DELPARAM(nl) != NULL)
+	    call mfree (NL_DELPARAM(nl), TY_PIXEL)
+	if (NL_PLIST(nl) != NULL)
+	    call mfree (NL_PLIST(nl), TY_INT)
+	if (NL_ALPHA(nl) != NULL)
+	    call mfree (NL_ALPHA(nl), TY_PIXEL)
+	if (NL_CHOFAC(nl) != NULL)
+	    call mfree (NL_CHOFAC(nl), TY_PIXEL)
+	if (NL_COVAR(nl) != NULL)
+	    call mfree (NL_COVAR(nl), TY_PIXEL)
+	if (NL_BETA(nl) != NULL)
+	    call mfree (NL_BETA(nl), TY_PIXEL)
+	if (NL_TRY(nl) != NULL)
+	    call mfree (NL_TRY(nl), TY_PIXEL)
+	if (NL_DERIV(nl) != NULL)
+	    call mfree (NL_DERIV(nl), TY_PIXEL)
+	call mfree (nl, TY_STRUCT)
+end
diff --git a/math/nlfit/nlfreed.x b/math/nlfit/nlfreed.x
new file mode 100644
index 00000000..56bcda62
--- /dev/null
+++ b/math/nlfit/nlfreed.x
@@ -0,0 +1,37 @@
+include "nlfitdefd.h"
+
+# NLFREE -- Deallocate all assigned space
+
+procedure nlfreed (nl)
+
+pointer	nl		# pointer to non-linear fitting structure
+
+errchk	mfree
+
+begin
+	if (nl == NULL)
+	    return
+	if (NL_PARAM(nl) != NULL)
+	    call mfree (NL_PARAM(nl), TY_DOUBLE)
+	if (NL_OPARAM(nl) != NULL)
+	    call mfree (NL_OPARAM(nl), TY_DOUBLE)
+	if (NL_DPARAM(nl) != NULL)
+	    call mfree (NL_DPARAM(nl), TY_DOUBLE)
+	if (NL_DELPARAM(nl) != NULL)
+	    call mfree (NL_DELPARAM(nl), TY_DOUBLE)
+	if (NL_PLIST(nl) != NULL)
+	    call mfree (NL_PLIST(nl), TY_INT)
+	if (NL_ALPHA(nl) != NULL)
+	    call mfree (NL_ALPHA(nl), TY_DOUBLE)
+	if (NL_CHOFAC(nl) != NULL)
+	    call mfree (NL_CHOFAC(nl), TY_DOUBLE)
+	if (NL_COVAR(nl) != NULL)
+	    call mfree (NL_COVAR(nl), TY_DOUBLE)
+	if (NL_BETA(nl) != NULL)
+	    call mfree (NL_BETA(nl), TY_DOUBLE)
+	if (NL_TRY(nl) != NULL)
+	    call mfree (NL_TRY(nl), TY_DOUBLE)
+	if (NL_DERIV(nl) != NULL)
+	    call mfree (NL_DERIV(nl), TY_DOUBLE)
+	call mfree (nl, TY_STRUCT)
+end
diff --git a/math/nlfit/nlfreer.x b/math/nlfit/nlfreer.x
new file mode 100644
index 00000000..639c4ecc
--- /dev/null
+++ b/math/nlfit/nlfreer.x
@@ -0,0 +1,37 @@
+include "nlfitdefr.h"
+
+# NLFREE -- Deallocate all assigned space
+
+procedure nlfreer (nl)
+
+pointer	nl		# pointer to non-linear fitting structure
+
+errchk	mfree
+
+begin
+	if (nl == NULL)
+	    return
+	if (NL_PARAM(nl) != NULL)
+	    call mfree (NL_PARAM(nl), TY_REAL)
+	if (NL_OPARAM(nl) != NULL)
+	    call mfree (NL_OPARAM(nl), TY_REAL)
+	if (NL_DPARAM(nl) != NULL)
+	    call mfree (NL_DPARAM(nl), TY_REAL)
+	if (NL_DELPARAM(nl) != NULL)
+	    call mfree (NL_DELPARAM(nl), TY_REAL)
+	if (NL_PLIST(nl) != NULL)
+	    call mfree (NL_PLIST(nl), TY_INT)
+	if (NL_ALPHA(nl) != NULL)
+	    call mfree (NL_ALPHA(nl), TY_REAL)
+	if (NL_CHOFAC(nl) != NULL)
+	    call mfree (NL_CHOFAC(nl), TY_REAL)
+	if (NL_COVAR(nl) != NULL)
+	    call mfree (NL_COVAR(nl), TY_REAL)
+	if (NL_BETA(nl) != NULL)
+	    call mfree (NL_BETA(nl), TY_REAL)
+	if (NL_TRY(nl) != NULL)
+	    call mfree (NL_TRY(nl), TY_REAL)
+	if (NL_DERIV(nl) != NULL)
+	    call mfree (NL_DERIV(nl), TY_REAL)
+	call mfree (nl, TY_STRUCT)
+end
diff --git a/math/nlfit/nlinit.gx b/math/nlfit/nlinit.gx
new file mode 100644
index 00000000..be6e436a
--- /dev/null
+++ b/math/nlfit/nlinit.gx
@@ -0,0 +1,66 @@
+$if (datatype == r)
+include "nlfitdefr.h"
+$else
+include "nlfitdefd.h"
+$endif
+
+# NLINIT --  Initialize for non-linear fitting
+
+procedure nlinit$t (nl, fnc, dfnc, params, dparams, nparams, plist, nfparams,
+		    tol, itmax)
+
+pointer	nl		# pointer to nl fitting structure
+int	fnc		# fitting function address
+int	dfnc		# derivative function address
+PIXEL	params[ARB]	# initial values for the parameters
+PIXEL	dparams[ARB]	# initial guess at uncertainties in parameters
+int	nparams		# number of parameters
+int	plist[ARB]	# list of active parameters
+int	nfparams	# number of fitted parameters
+PIXEL	tol		# fitting tolerance
+int	itmax		# maximum number of iterations
+
+errchk	malloc, calloc, nl_list
+
+begin
+	# Allocate space for the non-linear package structure.
+	call calloc (nl, LEN_NLSTRUCT, TY_STRUCT)
+
+	# Store the addresses of the non-linear functions.
+	NL_FUNC(nl) = fnc
+	NL_DFUNC(nl) = dfnc
+
+	# Allocate temporary space for arrays.
+	call calloc (NL_ALPHA(nl), nfparams * nfparams, TY_PIXEL)
+	call calloc (NL_COVAR(nl), nfparams * nfparams, TY_PIXEL)
+	call calloc (NL_CHOFAC(nl), nfparams * nfparams, TY_PIXEL)
+	call calloc (NL_BETA(nl), nfparams, TY_PIXEL)
+
+	# Allocate space for parameter and trial parameter vectors.
+	call calloc (NL_DERIV(nl), nparams, TY_PIXEL)
+	call calloc (NL_PARAM(nl), nparams, TY_PIXEL)
+	call calloc (NL_OPARAM(nl), nparams, TY_PIXEL)
+	call calloc (NL_TRY(nl), nparams, TY_PIXEL)
+	call calloc (NL_DPARAM(nl), nparams, TY_PIXEL)
+	call calloc (NL_DELPARAM(nl), nparams, TY_PIXEL)
+	call calloc (NL_PLIST(nl), nparams, TY_INT)
+
+	# Initialize the parameters.
+	NL_NPARAMS(nl) = nparams
+	NL_NFPARAMS(nl) = nfparams
+	call amov$t (params, PARAM(NL_PARAM(nl)), nparams)
+	call amov$t (params, PARAM(NL_OPARAM(nl)), nparams)
+	call amov$t (dparams, DPARAM(NL_DELPARAM(nl)), nparams)
+	call amovi (plist, PLIST(NL_PLIST(nl)), nfparams)
+	NL_TOL(nl) = tol
+	NL_ITMAX(nl) = itmax
+	NL_SCATTER(nl) = PIXEL(0.0)
+
+	# Set up the parameter list.
+	iferr {
+	    call nl_list (PLIST(NL_PLIST(nl)), NL_NPARAMS(nl), NL_NFPARAMS(nl))
+	} then {
+	    call nlfree$t (nl)
+	    nl = NULL
+	}
+end
diff --git a/math/nlfit/nlinitd.x b/math/nlfit/nlinitd.x
new file mode 100644
index 00000000..16804764
--- /dev/null
+++ b/math/nlfit/nlinitd.x
@@ -0,0 +1,62 @@
+include "nlfitdefd.h"
+
+# NLINIT --  Initialize for non-linear fitting
+
+procedure nlinitd (nl, fnc, dfnc, params, dparams, nparams, plist, nfparams,
+		    tol, itmax)
+
+pointer	nl		# pointer to nl fitting structure
+int	fnc		# fitting function address
+int	dfnc		# derivative function address
+double	params[ARB]	# initial values for the parameters
+double	dparams[ARB]	# initial guess at uncertainties in parameters
+int	nparams		# number of parameters
+int	plist[ARB]	# list of active parameters
+int	nfparams	# number of fitted parameters
+double	tol		# fitting tolerance
+int	itmax		# maximum number of iterations
+
+errchk	malloc, calloc, nl_list
+
+begin
+	# Allocate space for the non-linear package structure.
+	call calloc (nl, LEN_NLSTRUCT, TY_STRUCT)
+
+	# Store the addresses of the non-linear functions.
+	NL_FUNC(nl) = fnc
+	NL_DFUNC(nl) = dfnc
+
+	# Allocate temporary space for arrays.
+	call calloc (NL_ALPHA(nl), nfparams * nfparams, TY_DOUBLE)
+	call calloc (NL_COVAR(nl), nfparams * nfparams, TY_DOUBLE)
+	call calloc (NL_CHOFAC(nl), nfparams * nfparams, TY_DOUBLE)
+	call calloc (NL_BETA(nl), nfparams, TY_DOUBLE)
+
+	# Allocate space for parameter and trial parameter vectors.
+	call calloc (NL_DERIV(nl), nparams, TY_DOUBLE)
+	call calloc (NL_PARAM(nl), nparams, TY_DOUBLE)
+	call calloc (NL_OPARAM(nl), nparams, TY_DOUBLE)
+	call calloc (NL_TRY(nl), nparams, TY_DOUBLE)
+	call calloc (NL_DPARAM(nl), nparams, TY_DOUBLE)
+	call calloc (NL_DELPARAM(nl), nparams, TY_DOUBLE)
+	call calloc (NL_PLIST(nl), nparams, TY_INT)
+
+	# Initialize the parameters.
+	NL_NPARAMS(nl) = nparams
+	NL_NFPARAMS(nl) = nfparams
+	call amovd (params, PARAM(NL_PARAM(nl)), nparams)
+	call amovd (params, PARAM(NL_OPARAM(nl)), nparams)
+	call amovd (dparams, DPARAM(NL_DELPARAM(nl)), nparams)
+	call amovi (plist, PLIST(NL_PLIST(nl)), nfparams)
+	NL_TOL(nl) = tol
+	NL_ITMAX(nl) = itmax
+	NL_SCATTER(nl) = double(0.0)
+
+	# Set up the parameter list.
+	iferr {
+	    call nl_list (PLIST(NL_PLIST(nl)), NL_NPARAMS(nl), NL_NFPARAMS(nl))
+	} then {
+	    call nlfreed (nl)
+	    nl = NULL
+	}
+end
diff --git a/math/nlfit/nlinitr.x b/math/nlfit/nlinitr.x
new file mode 100644
index 00000000..cb97f269
--- /dev/null
+++ b/math/nlfit/nlinitr.x
@@ -0,0 +1,62 @@
+include "nlfitdefr.h"
+
+# NLINIT --  Initialize for non-linear fitting
+
+procedure nlinitr (nl, fnc, dfnc, params, dparams, nparams, plist, nfparams,
+		    tol, itmax)
+
+pointer	nl		# pointer to nl fitting structure
+int	fnc		# fitting function address
+int	dfnc		# derivative function address
+real	params[ARB]	# initial values for the parameters
+real	dparams[ARB]	# initial guess at uncertainties in parameters
+int	nparams		# number of parameters
+int	plist[ARB]	# list of active parameters
+int	nfparams	# number of fitted parameters
+real	tol		# fitting tolerance
+int	itmax		# maximum number of iterations
+
+errchk	malloc, calloc, nl_list
+
+begin
+	# Allocate space for the non-linear package structure.
+	call calloc (nl, LEN_NLSTRUCT, TY_STRUCT)
+
+	# Store the addresses of the non-linear functions.
+	NL_FUNC(nl) = fnc
+	NL_DFUNC(nl) = dfnc
+
+	# Allocate temporary space for arrays.
+	call calloc (NL_ALPHA(nl), nfparams * nfparams, TY_REAL)
+	call calloc (NL_COVAR(nl), nfparams * nfparams, TY_REAL)
+	call calloc (NL_CHOFAC(nl), nfparams * nfparams, TY_REAL)
+	call calloc (NL_BETA(nl), nfparams, TY_REAL)
+
+	# Allocate space for parameter and trial parameter vectors.
+	call calloc (NL_DERIV(nl), nparams, TY_REAL)
+	call calloc (NL_PARAM(nl), nparams, TY_REAL)
+	call calloc (NL_OPARAM(nl), nparams, TY_REAL)
+	call calloc (NL_TRY(nl), nparams, TY_REAL)
+	call calloc (NL_DPARAM(nl), nparams, TY_REAL)
+	call calloc (NL_DELPARAM(nl), nparams, TY_REAL)
+	call calloc (NL_PLIST(nl), nparams, TY_INT)
+
+	# Initialize the parameters.
+	NL_NPARAMS(nl) = nparams
+	NL_NFPARAMS(nl) = nfparams
+	call amovr (params, PARAM(NL_PARAM(nl)), nparams)
+	call amovr (params, PARAM(NL_OPARAM(nl)), nparams)
+	call amovr (dparams, DPARAM(NL_DELPARAM(nl)), nparams)
+	call amovi (plist, PLIST(NL_PLIST(nl)), nfparams)
+	NL_TOL(nl) = tol
+	NL_ITMAX(nl) = itmax
+	NL_SCATTER(nl) = real(0.0)
+
+	# Set up the parameter list.
+	iferr {
+	    call nl_list (PLIST(NL_PLIST(nl)), NL_NPARAMS(nl), NL_NFPARAMS(nl))
+	} then {
+	    call nlfreer (nl)
+	    nl = NULL
+	}
+end
diff --git a/math/nlfit/nliter.gx b/math/nlfit/nliter.gx
new file mode 100644
index 00000000..7857f852
--- /dev/null
+++ b/math/nlfit/nliter.gx
@@ -0,0 +1,59 @@
+include <mach.h>
+include	<math/nlfit.h>
+$if (datatype == r)
+include "nlfitdefr.h"
+$else
+include "nlfitdefd.h"
+$endif
+
+# NLITER - Routine to compute one iteration of the fitting process
+
+procedure nliter$t (nl, x, z, w, npts, nvars, ier)
+
+pointer	nl		# pointer to nl fitting structure
+PIXEL	x[ARB]		# independent variables (npts * nvars)
+PIXEL	z[ARB]		# function values (npts)
+PIXEL	w[ARB]		# weights (npts)
+int	npts		# number of points
+int	nvars		# number of independent variables
+int	ier		# error code
+
+int	i, index
+PIXEL	nlacpts$t(), nlresid$t()
+
+begin
+	# Do the initial accumulation.
+   	NL_OLDSQ(nl) = nlacpts$t (nl, x, z, w, npts, nvars)
+
+	# Set up temporary parameter array.
+	call amov$t (PARAM(NL_PARAM(nl)), TRY(NL_TRY(nl)), NL_NPARAMS(nl))
+
+	repeat {
+
+	    # Solve the matrix.
+	    call nlsolve$t (nl, ier)
+	    if (ier == NO_DEG_FREEDOM)
+	        return
+
+	    # Increment the parameters.
+	    do i = 1, NL_NFPARAMS(nl) {
+	        index = PLIST(NL_PLIST(nl)+i-1)
+	        TRY(NL_TRY(nl)+index-1) = TRY(NL_TRY(nl)+index-1) +
+		    DPARAM(NL_DPARAM(nl)+i-1)
+	    }
+
+	    # Reaccumulate the residuals and increment lambda.
+	    NL_SUMSQ(nl) = nlresid$t (nl, x, z, w, npts, nvars)
+	    #if (NL_OLDSQ(nl) > (NL_SUMSQ(nl) + NL_TOL(nl) * NL_SUMSQ(nl))) {
+	    if (NL_OLDSQ(nl) >= NL_SUMSQ(nl)) {
+		call amov$t (TRY(NL_TRY(nl)), PARAM(NL_PARAM(nl)),
+		    NL_NPARAMS(nl))
+		NL_LAMBDA(nl) = PIXEL (0.10) * NL_LAMBDA(nl)
+		break
+	    } else
+	        NL_LAMBDA(nl) = min (PIXEL(LAMBDAMAX),
+		    PIXEL (10.0) * NL_LAMBDA(nl))
+
+	} until (NL_LAMBDA(nl) <= PIXEL(0.0) ||
+	    NL_LAMBDA(nl) >= PIXEL(LAMBDAMAX))
+end
diff --git a/math/nlfit/nliterd.x b/math/nlfit/nliterd.x
new file mode 100644
index 00000000..418ab8fa
--- /dev/null
+++ b/math/nlfit/nliterd.x
@@ -0,0 +1,55 @@
+include <mach.h>
+include	<math/nlfit.h>
+include "nlfitdefd.h"
+
+# NLITER - Routine to compute one iteration of the fitting process
+
+procedure nliterd (nl, x, z, w, npts, nvars, ier)
+
+pointer	nl		# pointer to nl fitting structure
+double	x[ARB]		# independent variables (npts * nvars)
+double	z[ARB]		# function values (npts)
+double	w[ARB]		# weights (npts)
+int	npts		# number of points
+int	nvars		# number of independent variables
+int	ier		# error code
+
+int	i, index
+double	nlacptsd(), nlresidd()
+
+begin
+	# Do the initial accumulation.
+   	NL_OLDSQ(nl) = nlacptsd (nl, x, z, w, npts, nvars)
+
+	# Set up temporary parameter array.
+	call amovd (PARAM(NL_PARAM(nl)), TRY(NL_TRY(nl)), NL_NPARAMS(nl))
+
+	repeat {
+
+	    # Solve the matrix.
+	    call nlsolved (nl, ier)
+	    if (ier == NO_DEG_FREEDOM)
+	        return
+
+	    # Increment the parameters.
+	    do i = 1, NL_NFPARAMS(nl) {
+	        index = PLIST(NL_PLIST(nl)+i-1)
+	        TRY(NL_TRY(nl)+index-1) = TRY(NL_TRY(nl)+index-1) +
+		    DPARAM(NL_DPARAM(nl)+i-1)
+	    }
+
+	    # Reaccumulate the residuals and increment lambda.
+	    NL_SUMSQ(nl) = nlresidd (nl, x, z, w, npts, nvars)
+	    #if (NL_OLDSQ(nl) > (NL_SUMSQ(nl) + NL_TOL(nl) * NL_SUMSQ(nl))) {
+	    if (NL_OLDSQ(nl) >= NL_SUMSQ(nl)) {
+		call amovd (TRY(NL_TRY(nl)), PARAM(NL_PARAM(nl)),
+		    NL_NPARAMS(nl))
+		NL_LAMBDA(nl) = double (0.10) * NL_LAMBDA(nl)
+		break
+	    } else
+	        NL_LAMBDA(nl) = min (double(LAMBDAMAX),
+		    double (10.0) * NL_LAMBDA(nl))
+
+	} until (NL_LAMBDA(nl) <= double(0.0) ||
+	    NL_LAMBDA(nl) >= double(LAMBDAMAX))
+end
diff --git a/math/nlfit/nliterr.x b/math/nlfit/nliterr.x
new file mode 100644
index 00000000..cbc5a80a
--- /dev/null
+++ b/math/nlfit/nliterr.x
@@ -0,0 +1,55 @@
+include <mach.h>
+include	<math/nlfit.h>
+include "nlfitdefr.h"
+
+# NLITER - Routine to compute one iteration of the fitting process
+
+procedure nliterr (nl, x, z, w, npts, nvars, ier)
+
+pointer	nl		# pointer to nl fitting structure
+real	x[ARB]		# independent variables (npts * nvars)
+real	z[ARB]		# function values (npts)
+real	w[ARB]		# weights (npts)
+int	npts		# number of points
+int	nvars		# number of independent variables
+int	ier		# error code
+
+int	i, index
+real	nlacptsr(), nlresidr()
+
+begin
+	# Do the initial accumulation.
+   	NL_OLDSQ(nl) = nlacptsr (nl, x, z, w, npts, nvars)
+
+	# Set up temporary parameter array.
+	call amovr (PARAM(NL_PARAM(nl)), TRY(NL_TRY(nl)), NL_NPARAMS(nl))
+
+	repeat {
+
+	    # Solve the matrix.
+	    call nlsolver (nl, ier)
+	    if (ier == NO_DEG_FREEDOM)
+	        return
+
+	    # Increment the parameters.
+	    do i = 1, NL_NFPARAMS(nl) {
+	        index = PLIST(NL_PLIST(nl)+i-1)
+	        TRY(NL_TRY(nl)+index-1) = TRY(NL_TRY(nl)+index-1) +
+		    DPARAM(NL_DPARAM(nl)+i-1)
+	    }
+
+	    # Reaccumulate the residuals and increment lambda.
+	    NL_SUMSQ(nl) = nlresidr (nl, x, z, w, npts, nvars)
+	    #if (NL_OLDSQ(nl) > (NL_SUMSQ(nl) + NL_TOL(nl) * NL_SUMSQ(nl))) {
+	    if (NL_OLDSQ(nl) >= NL_SUMSQ(nl)) {
+		call amovr (TRY(NL_TRY(nl)), PARAM(NL_PARAM(nl)),
+		    NL_NPARAMS(nl))
+		NL_LAMBDA(nl) = real (0.10) * NL_LAMBDA(nl)
+		break
+	    } else
+	        NL_LAMBDA(nl) = min (real(LAMBDAMAX),
+		    real (10.0) * NL_LAMBDA(nl))
+
+	} until (NL_LAMBDA(nl) <= real(0.0) ||
+	    NL_LAMBDA(nl) >= real(LAMBDAMAX))
+end
diff --git a/math/nlfit/nllist.x b/math/nlfit/nllist.x
new file mode 100644
index 00000000..d96351c9
--- /dev/null
+++ b/math/nlfit/nllist.x
@@ -0,0 +1,30 @@
+# NL_LIST -- Procedure to order the list used when the NLFIT structure
+# is initialized.
+
+procedure nl_list (list, nlist, nfit)
+
+int	list[ARB]		# list
+int	nlist			# number of elements in the list
+int	nfit			# number of active list elments
+
+int	i, j, nfitp1, ifound
+
+begin
+	nfitp1 = nfit + 1
+
+	do i = 1, nlist {
+	    ifound = 0
+	    do j = 1, nfit {
+		if (list[j] == i)
+		    ifound = ifound + 1
+	    }
+	    if (ifound == 0) {
+		list[nfitp1] = i
+		nfitp1 = nfitp1 + 1
+	    } else if (ifound > 1)
+	        call error (0, "Incorrect parameter ordering in plist")
+	}
+
+	if (nfitp1 != (nlist + 1))
+	    call error (0, "Incorrect parameter ordering in plist")
+end
diff --git a/math/nlfit/nlpget.gx b/math/nlfit/nlpget.gx
new file mode 100644
index 00000000..3a74fa6f
--- /dev/null
+++ b/math/nlfit/nlpget.gx
@@ -0,0 +1,18 @@
+$if (datatype == r)
+include "nlfitdefr.h"
+$else
+include "nlfitdefd.h"
+$endif
+
+# NLPGET - Retreive parameter values
+
+procedure nlpget$t (nl, params, nparams)
+
+pointer	nl			# pointer to the nlfit structure
+PIXEL	params[ARB]		# parameter array
+int	nparams			# the number of the parameters
+
+begin
+	nparams = NL_NPARAMS(nl)
+	call amov$t (PARAM(NL_PARAM(nl)), params, nparams)
+end
diff --git a/math/nlfit/nlpgetd.x b/math/nlfit/nlpgetd.x
new file mode 100644
index 00000000..419bde09
--- /dev/null
+++ b/math/nlfit/nlpgetd.x
@@ -0,0 +1,14 @@
+include "nlfitdefd.h"
+
+# NLPGET - Retreive parameter values
+
+procedure nlpgetd (nl, params, nparams)
+
+pointer	nl			# pointer to the nlfit structure
+double	params[ARB]		# parameter array
+int	nparams			# the number of the parameters
+
+begin
+	nparams = NL_NPARAMS(nl)
+	call amovd (PARAM(NL_PARAM(nl)), params, nparams)
+end
diff --git a/math/nlfit/nlpgetr.x b/math/nlfit/nlpgetr.x
new file mode 100644
index 00000000..73fd7cbb
--- /dev/null
+++ b/math/nlfit/nlpgetr.x
@@ -0,0 +1,14 @@
+include "nlfitdefr.h"
+
+# NLPGET - Retreive parameter values
+
+procedure nlpgetr (nl, params, nparams)
+
+pointer	nl			# pointer to the nlfit structure
+real	params[ARB]		# parameter array
+int	nparams			# the number of the parameters
+
+begin
+	nparams = NL_NPARAMS(nl)
+	call amovr (PARAM(NL_PARAM(nl)), params, nparams)
+end
diff --git a/math/nlfit/nlsolve.gx b/math/nlfit/nlsolve.gx
new file mode 100644
index 00000000..a4055520
--- /dev/null
+++ b/math/nlfit/nlsolve.gx
@@ -0,0 +1,41 @@
+include <math/nlfit.h>
+$if (datatype == r)
+include "nlfitdefr.h"
+$else
+include "nlfitdefd.h"
+$endif
+
+# NLSOLVE -- Procedure to solve nonlinear system
+
+procedure nlsolve$t (nl, ier)
+
+pointer	nl	# pointer to the nlfit structure
+int	ier	# error code
+
+int	nfree
+
+begin
+	# Make temporary arrays.
+	call amov$t (ALPHA(NL_ALPHA(nl)), COVAR(NL_COVAR(nl)),
+	    NL_NFPARAMS(nl) ** 2)
+	call amov$t (BETA(NL_BETA(nl)), DPARAM(NL_DPARAM(nl)), NL_NFPARAMS(nl))
+
+	# Add the lambda damping factor.
+	call nl_damp$t (COVAR(NL_COVAR(nl)), COVAR(NL_COVAR(nl)),
+	    (1.0 + NL_LAMBDA(nl)), NL_NFPARAMS(nl), NL_NFPARAMS(nl))
+
+	ier = OK
+	nfree = NL_NPTS(nl) - NL_NFPARAMS(nl)
+	if (nfree < 0) {
+	    ier = NO_DEG_FREEDOM
+	    return
+	}
+
+	# Compute the factorization of the matrix.
+	call nl_chfac$t (COVAR(NL_COVAR(nl)), NL_NFPARAMS(nl), NL_NFPARAMS(nl),
+	    CHOFAC(NL_CHOFAC(nl)), ier)
+
+	# Solve the equations for the parameter increments.
+	call nl_chslv$t (CHOFAC(NL_CHOFAC(nl)), NL_NFPARAMS(nl), NL_NFPARAMS(nl),
+	    DPARAM(NL_DPARAM(nl)), DPARAM(NL_DPARAM(nl)))
+end
diff --git a/math/nlfit/nlsolved.x b/math/nlfit/nlsolved.x
new file mode 100644
index 00000000..77a43337
--- /dev/null
+++ b/math/nlfit/nlsolved.x
@@ -0,0 +1,37 @@
+include <math/nlfit.h>
+include "nlfitdefd.h"
+
+# NLSOLVE -- Procedure to solve nonlinear system
+
+procedure nlsolved (nl, ier)
+
+pointer	nl	# pointer to the nlfit structure
+int	ier	# error code
+
+int	nfree
+
+begin
+	# Make temporary arrays.
+	call amovd (ALPHA(NL_ALPHA(nl)), COVAR(NL_COVAR(nl)),
+	    NL_NFPARAMS(nl) ** 2)
+	call amovd (BETA(NL_BETA(nl)), DPARAM(NL_DPARAM(nl)), NL_NFPARAMS(nl))
+
+	# Add the lambda damping factor.
+	call nl_dampd (COVAR(NL_COVAR(nl)), COVAR(NL_COVAR(nl)),
+	    (1.0 + NL_LAMBDA(nl)), NL_NFPARAMS(nl), NL_NFPARAMS(nl))
+
+	ier = OK
+	nfree = NL_NPTS(nl) - NL_NFPARAMS(nl)
+	if (nfree < 0) {
+	    ier = NO_DEG_FREEDOM
+	    return
+	}
+
+	# Compute the factorization of the matrix.
+	call nl_chfacd (COVAR(NL_COVAR(nl)), NL_NFPARAMS(nl), NL_NFPARAMS(nl),
+	    CHOFAC(NL_CHOFAC(nl)), ier)
+
+	# Solve the equations for the parameter increments.
+	call nl_chslvd (CHOFAC(NL_CHOFAC(nl)), NL_NFPARAMS(nl), NL_NFPARAMS(nl),
+	    DPARAM(NL_DPARAM(nl)), DPARAM(NL_DPARAM(nl)))
+end
diff --git a/math/nlfit/nlsolver.x b/math/nlfit/nlsolver.x
new file mode 100644
index 00000000..85d238d6
--- /dev/null
+++ b/math/nlfit/nlsolver.x
@@ -0,0 +1,37 @@
+include <math/nlfit.h>
+include "nlfitdefr.h"
+
+# NLSOLVE -- Procedure to solve nonlinear system
+
+procedure nlsolver (nl, ier)
+
+pointer	nl	# pointer to the nlfit structure
+int	ier	# error code
+
+int	nfree
+
+begin
+	# Make temporary arrays.
+	call amovr (ALPHA(NL_ALPHA(nl)), COVAR(NL_COVAR(nl)),
+	    NL_NFPARAMS(nl) ** 2)
+	call amovr (BETA(NL_BETA(nl)), DPARAM(NL_DPARAM(nl)), NL_NFPARAMS(nl))
+
+	# Add the lambda damping factor.
+	call nl_dampr (COVAR(NL_COVAR(nl)), COVAR(NL_COVAR(nl)),
+	    (1.0 + NL_LAMBDA(nl)), NL_NFPARAMS(nl), NL_NFPARAMS(nl))
+
+	ier = OK
+	nfree = NL_NPTS(nl) - NL_NFPARAMS(nl)
+	if (nfree < 0) {
+	    ier = NO_DEG_FREEDOM
+	    return
+	}
+
+	# Compute the factorization of the matrix.
+	call nl_chfacr (COVAR(NL_COVAR(nl)), NL_NFPARAMS(nl), NL_NFPARAMS(nl),
+	    CHOFAC(NL_CHOFAC(nl)), ier)
+
+	# Solve the equations for the parameter increments.
+	call nl_chslvr (CHOFAC(NL_CHOFAC(nl)), NL_NFPARAMS(nl), NL_NFPARAMS(nl),
+	    DPARAM(NL_DPARAM(nl)), DPARAM(NL_DPARAM(nl)))
+end
diff --git a/math/nlfit/nlstat.gx b/math/nlfit/nlstat.gx
new file mode 100644
index 00000000..c17d5940
--- /dev/null
+++ b/math/nlfit/nlstat.gx
@@ -0,0 +1,30 @@
+include <math/nlfit.h>
+$if (datatype == r)
+include "nlfitdefr.h"
+$else
+include "nlfitdefd.h"
+$endif
+
+# NLSTAT[RD] - Fetch an NLFIT real/double parameter
+
+PIXEL procedure nlstat$t (nl, param)
+
+pointer	nl		# pointer to NLFIT structure
+int	param		# parameter to be fetched
+
+begin
+	switch (param) {
+	case NLSUMSQ:
+	    return (NL_SUMSQ(nl))
+	case NLOLDSQ:
+	    return (NL_OLDSQ(nl))
+	case NLTOL:
+	    return (NL_TOL(nl))
+	case NLLAMBDA:
+	    return (NL_LAMBDA(nl))
+	case NLSCATTER:
+	    return (NL_SCATTER(nl))
+	default:
+	    return ($INDEF$T)
+	}
+end
diff --git a/math/nlfit/nlstatd.x b/math/nlfit/nlstatd.x
new file mode 100644
index 00000000..73661ce8
--- /dev/null
+++ b/math/nlfit/nlstatd.x
@@ -0,0 +1,26 @@
+include <math/nlfit.h>
+include "nlfitdefd.h"
+
+# NLSTAT[RD] - Fetch an NLFIT real/double parameter
+
+double procedure nlstatd (nl, param)
+
+pointer	nl		# pointer to NLFIT structure
+int	param		# parameter to be fetched
+
+begin
+	switch (param) {
+	case NLSUMSQ:
+	    return (NL_SUMSQ(nl))
+	case NLOLDSQ:
+	    return (NL_OLDSQ(nl))
+	case NLTOL:
+	    return (NL_TOL(nl))
+	case NLLAMBDA:
+	    return (NL_LAMBDA(nl))
+	case NLSCATTER:
+	    return (NL_SCATTER(nl))
+	default:
+	    return (INDEFD)
+	}
+end
diff --git a/math/nlfit/nlstati.x b/math/nlfit/nlstati.x
new file mode 100644
index 00000000..04c8cf27
--- /dev/null
+++ b/math/nlfit/nlstati.x
@@ -0,0 +1,26 @@
+include "nlfitdefr.h"
+include <math/nlfit.h>
+
+# NLSTATI - Fetch a NLFIT integer parameter
+
+int procedure nlstati (nl, param)
+
+pointer	nl		# pointer to NLFIT structure
+int	param		# parameter to be fetched
+
+begin
+	switch (param) {
+	case NLNPARAMS:
+	    return (NL_NPARAMS(nl))
+	case NLNFPARAMS:
+	    return (NL_NFPARAMS(nl))
+	case NLITMAX:
+	    return (NL_ITMAX(nl))
+	case NLITER:
+	    return (NL_ITER(nl))
+	case NLNPTS:
+	    return (NL_NPTS(nl))
+	default:
+	    return (INDEFI)
+	}
+end
diff --git a/math/nlfit/nlstatr.x b/math/nlfit/nlstatr.x
new file mode 100644
index 00000000..c6360017
--- /dev/null
+++ b/math/nlfit/nlstatr.x
@@ -0,0 +1,26 @@
+include <math/nlfit.h>
+include "nlfitdefr.h"
+
+# NLSTAT[RD] - Fetch an NLFIT real/double parameter
+
+real procedure nlstatr (nl, param)
+
+pointer	nl		# pointer to NLFIT structure
+int	param		# parameter to be fetched
+
+begin
+	switch (param) {
+	case NLSUMSQ:
+	    return (NL_SUMSQ(nl))
+	case NLOLDSQ:
+	    return (NL_OLDSQ(nl))
+	case NLTOL:
+	    return (NL_TOL(nl))
+	case NLLAMBDA:
+	    return (NL_LAMBDA(nl))
+	case NLSCATTER:
+	    return (NL_SCATTER(nl))
+	default:
+	    return (INDEFR)
+	}
+end
diff --git a/math/nlfit/nlvector.gx b/math/nlfit/nlvector.gx
new file mode 100644
index 00000000..f33d063c
--- /dev/null
+++ b/math/nlfit/nlvector.gx
@@ -0,0 +1,27 @@
+include	<math/nlfit.h>
+$if (datatype == r)
+include "nlfitdefr.h"
+$else
+include "nlfitdefd.h"
+$endif
+
+define	VAR	(($1 - 1) * $2 + 1)
+
+# NLVECTOR -- Evaluate the fit for a series of data points.
+
+procedure nlvector$t (nl, x, zfit, npts, nvars)
+
+pointer	nl		# pointer to nl fitting structure
+PIXEL	x[ARB]		# independent variables (npts * nvars)
+PIXEL	zfit[ARB]	# function values (npts)
+int	npts		# number of points
+int	nvars		# number of independent variables
+
+int	i
+
+begin
+	# Compute the fitted function.
+	do i = 1, npts
+	    call zcall5 (NL_FUNC(nl), x[VAR (i, nvars)], nvars,
+		PARAM(NL_PARAM(nl)), NL_NPARAMS(nl), zfit[i])
+end
diff --git a/math/nlfit/nlvectord.x b/math/nlfit/nlvectord.x
new file mode 100644
index 00000000..34bce9a5
--- /dev/null
+++ b/math/nlfit/nlvectord.x
@@ -0,0 +1,23 @@
+include	<math/nlfit.h>
+include "nlfitdefd.h"
+
+define	VAR	(($1 - 1) * $2 + 1)
+
+# NLVECTOR -- Evaluate the fit for a series of data points.
+
+procedure nlvectord (nl, x, zfit, npts, nvars)
+
+pointer	nl		# pointer to nl fitting structure
+double	x[ARB]		# independent variables (npts * nvars)
+double	zfit[ARB]	# function values (npts)
+int	npts		# number of points
+int	nvars		# number of independent variables
+
+int	i
+
+begin
+	# Compute the fitted function.
+	do i = 1, npts
+	    call zcall5 (NL_FUNC(nl), x[VAR (i, nvars)], nvars,
+		PARAM(NL_PARAM(nl)), NL_NPARAMS(nl), zfit[i])
+end
diff --git a/math/nlfit/nlvectorr.x b/math/nlfit/nlvectorr.x
new file mode 100644
index 00000000..43450695
--- /dev/null
+++ b/math/nlfit/nlvectorr.x
@@ -0,0 +1,23 @@
+include	<math/nlfit.h>
+include "nlfitdefr.h"
+
+define	VAR	(($1 - 1) * $2 + 1)
+
+# NLVECTOR -- Evaluate the fit for a series of data points.
+
+procedure nlvectorr (nl, x, zfit, npts, nvars)
+
+pointer	nl		# pointer to nl fitting structure
+real	x[ARB]		# independent variables (npts * nvars)
+real	zfit[ARB]	# function values (npts)
+int	npts		# number of points
+int	nvars		# number of independent variables
+
+int	i
+
+begin
+	# Compute the fitted function.
+	do i = 1, npts
+	    call zcall5 (NL_FUNC(nl), x[VAR (i, nvars)], nvars,
+		PARAM(NL_PARAM(nl)), NL_NPARAMS(nl), zfit[i])
+end
diff --git a/math/nlfit/nlzero.gx b/math/nlfit/nlzero.gx
new file mode 100644
index 00000000..ffd1458b
--- /dev/null
+++ b/math/nlfit/nlzero.gx
@@ -0,0 +1,38 @@
+$if (datatype == r)
+include "nlfitdefr.h"
+$else
+include "nlfitdefd.h"
+$endif
+
+
+# NLZERO --  Zero the accumulators and reset the fitting parameter values to
+# their original values set by nlinit().
+
+procedure nlzero$t (nl)
+
+pointer	nl		# pointer to nl fitting structure
+
+int	nparams		# number of parameters
+int	nfparams	# number of fitted parameters
+
+begin
+	# Get number of parameters and fitting parameters.
+	nparams  = NL_NPARAMS(nl)
+	nfparams = NL_NFPARAMS(nl)
+
+	# Clear temporary array space.
+	call aclr$t (ALPHA(NL_ALPHA(nl)), nfparams * nfparams)
+	call aclr$t (COVAR(NL_COVAR(nl)), nfparams * nfparams)
+	call aclr$t (CHOFAC(NL_CHOFAC(nl)), nfparams * nfparams)
+	call aclr$t (BETA(NL_BETA(nl)), nfparams)
+
+	# Clear space for derivatives and trial parameter vectors.
+	call aclr$t (DERIV(NL_DERIV(nl)), nparams)
+	call aclr$t (TRY(NL_TRY(nl)), nparams)
+
+	# Reset parameters.
+	call amov$t (OPARAM(NL_OPARAM(nl)), PARAM(NL_PARAM(nl)), nparams)
+	call aclr$t (DPARAM(NL_DPARAM(nl)), nparams)
+
+	NL_SCATTER(nl) = PIXEL(0.0)
+end
diff --git a/math/nlfit/nlzerod.x b/math/nlfit/nlzerod.x
new file mode 100644
index 00000000..2ee86562
--- /dev/null
+++ b/math/nlfit/nlzerod.x
@@ -0,0 +1,34 @@
+include "nlfitdefd.h"
+
+
+# NLZERO --  Zero the accumulators and reset the fitting parameter values to
+# their original values set by nlinit().
+
+procedure nlzerod (nl)
+
+pointer	nl		# pointer to nl fitting structure
+
+int	nparams		# number of parameters
+int	nfparams	# number of fitted parameters
+
+begin
+	# Get number of parameters and fitting parameters.
+	nparams  = NL_NPARAMS(nl)
+	nfparams = NL_NFPARAMS(nl)
+
+	# Clear temporary array space.
+	call aclrd (ALPHA(NL_ALPHA(nl)), nfparams * nfparams)
+	call aclrd (COVAR(NL_COVAR(nl)), nfparams * nfparams)
+	call aclrd (CHOFAC(NL_CHOFAC(nl)), nfparams * nfparams)
+	call aclrd (BETA(NL_BETA(nl)), nfparams)
+
+	# Clear space for derivatives and trial parameter vectors.
+	call aclrd (DERIV(NL_DERIV(nl)), nparams)
+	call aclrd (TRY(NL_TRY(nl)), nparams)
+
+	# Reset parameters.
+	call amovd (OPARAM(NL_OPARAM(nl)), PARAM(NL_PARAM(nl)), nparams)
+	call aclrd (DPARAM(NL_DPARAM(nl)), nparams)
+
+	NL_SCATTER(nl) = double(0.0)
+end
diff --git a/math/nlfit/nlzeror.x b/math/nlfit/nlzeror.x
new file mode 100644
index 00000000..1cceba31
--- /dev/null
+++ b/math/nlfit/nlzeror.x
@@ -0,0 +1,34 @@
+include "nlfitdefr.h"
+
+
+# NLZERO --  Zero the accumulators and reset the fitting parameter values to
+# their original values set by nlinit().
+
+procedure nlzeror (nl)
+
+pointer	nl		# pointer to nl fitting structure
+
+int	nparams		# number of parameters
+int	nfparams	# number of fitted parameters
+
+begin
+	# Get number of parameters and fitting parameters.
+	nparams  = NL_NPARAMS(nl)
+	nfparams = NL_NFPARAMS(nl)
+
+	# Clear temporary array space.
+	call aclrr (ALPHA(NL_ALPHA(nl)), nfparams * nfparams)
+	call aclrr (COVAR(NL_COVAR(nl)), nfparams * nfparams)
+	call aclrr (CHOFAC(NL_CHOFAC(nl)), nfparams * nfparams)
+	call aclrr (BETA(NL_BETA(nl)), nfparams)
+
+	# Clear space for derivatives and trial parameter vectors.
+	call aclrr (DERIV(NL_DERIV(nl)), nparams)
+	call aclrr (TRY(NL_TRY(nl)), nparams)
+
+	# Reset parameters.
+	call amovr (OPARAM(NL_OPARAM(nl)), PARAM(NL_PARAM(nl)), nparams)
+	call aclrr (DPARAM(NL_DPARAM(nl)), nparams)
+
+	NL_SCATTER(nl) = real(0.0)
+end