Repatch (from linux) of OSX IRAF

1 files changed, 107 insertions, 0 deletions

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