Repatch (from linux) of OSX IRAF

1 files changed, 171 insertions, 0 deletions

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