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