Initial commit

1 files changed, 111 insertions, 0 deletions

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