Repatch (from linux) of OSX IRAF

1 files changed, 88 insertions, 0 deletions

diff --git a/noao/digiphot/apphot/center/apmctr1d.x b/noao/digiphot/apphot/center/apmctr1d.x
new file mode 100644
index 00000000..5dc47f60
--- /dev/null
+++ b/noao/digiphot/apphot/center/apmctr1d.x
@@ -0,0 +1,88 @@
+include <mach.h>
+include "../lib/center.h"
+
+# AP_MCTR1D -- Procedure to compute the center from the 1D marginal
+# distributions.
+
+int procedure ap_mctr1d (ctrpix, nx, ny, norm, xc, yc, xerr, yerr)
+
+real	ctrpix[nx, ny]		# object to be centered
+int	nx, ny			# dimensions of subarray
+real	norm			# the normalization factor
+real	xc, yc			# computed centers
+real	xerr, yerr		# estimate of centering error
+
+pointer	sp, xm, ym
+
+begin
+	call smark (sp)
+	call salloc (xm, nx, TY_REAL)
+	call salloc (ym, ny, TY_REAL)
+
+	# Compute the marginal distributions.
+	call ap_mkmarg (ctrpix, Memr[xm], Memr[ym], nx, ny)
+	call adivkr (Memr[xm], real (ny), Memr[xm], nx)
+	call adivkr (Memr[ym], real (nx), Memr[ym], ny)
+
+	# Get the centers and errors.
+	call ap_cmmarg (Memr[xm], nx, norm, xc, xerr)
+	call ap_cmmarg (Memr[ym], ny, norm, yc, yerr)
+
+	call sfree (sp)
+	return (AP_OK)
+end
+
+
+# AP_CMMARG -- Compute the center and estimate its error given the
+# marginal distribution and the number of points. The center is estimated
+# using a centroid which is an optimal solution for a gaussian star with
+# no background noise. The centering error is computing using a variance
+# estimate based on a gaussian star with no background noise.
+
+procedure ap_cmmarg (a, npts, norm, xc, err)
+
+real	a[npts]		# array
+int	npts		# number of points
+real	norm		# the normalization factor
+real	xc		# center value
+real	err		# error
+
+int	i, npos
+real	sumi, sumix, sumix2, mean, val
+real	asumr()
+
+begin
+	# Initialize.
+	mean = asumr (a, npts) / npts
+	npos = 0
+	sumi = 0.0
+	sumix = 0.0
+	sumix2 = 0.0
+
+	# Accumulate the sums.
+	do i = 1, npts {
+	    val = (a[i] - mean)
+	    if (val > 0.0) {
+		npos = npos + 1
+	        sumi = sumi + val
+	        sumix = sumix + val * i
+	        sumix2 = sumix2 + val * i ** 2
+	    }
+	}
+
+	# Compute the position and the error.
+	if (npos <= 0) {
+	    xc =  (1.0 + npts) / 2.0
+	    err = INDEFR
+	} else {
+	    xc = sumix / sumi
+	    err = (sumix2 / sumi - xc * xc)
+	    if (err <= 0.0) {
+		err = 0.0
+	    } else {
+	        err = sqrt (err / ((sumi + mean * npts) * norm))
+		if (err > real (npts))
+		    err = INDEFR
+	    }
+	}
+end