1 files changed, 87 insertions, 0 deletions
diff --git a/noao/digiphot/apphot/center/aplgctr1d.x b/noao/digiphot/apphot/center/aplgctr1d.x
new file mode 100644
index 00000000..ce6d7c53
--- /dev/null
+++ b/noao/digiphot/apphot/center/aplgctr1d.x
@@ -0,0 +1,87 @@
+include <math.h>
+include "../lib/center.h"
+
+define	TOL	0.001		# tolerance for fitting algorithm
+
+# AP_LGCTR1D -- Procedure to compute the center from the 1D marginal
+# distributions using a simplified version of the optimal filtering
+# technique and addopting a Gaussian model for the fit. The method
+# is streamlined by replacing the Gaussian with a simple triangle
+# following L.Goad.
+
+int procedure ap_lgctr1d (ctrpix, nx, ny, cx, cy, sigma, maxiter, norm,
+	skysigma, xc, yc, xerr, yerr)
+
+real	ctrpix[nx, ny]		# object to be centered
+int	nx, ny			# dimensions of subarray
+real	cx, cy			# center in subraster coordinates
+real	sigma			# sigma of PSF
+int	maxiter			# maximum number of iterations
+real	norm			# the normalization factor
+real	skysigma		# standard deviation of the pixels
+real	xc, yc			# computed centers
+real	xerr, yerr		# first guess at errors
+
+int	nxiter, nyiter
+pointer	sp, xm, ym
+real	ratio, constant
+
+int	aptopt()
+real	asumr()
+
+begin
+	# Allocate working space.
+	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)
+	xerr = asumr (Memr[xm], nx)
+	yerr = asumr (Memr[ym], ny)
+	call adivkr (Memr[xm], real (nx), Memr[xm], nx)
+	call adivkr (Memr[ym], real (ny), Memr[ym], ny)
+
+	# Compute the x center and error.
+	xc = cx
+	nxiter = aptopt (Memr[xm], nx, xc, sigma, TOL, maxiter, YES)
+	if (xerr <= 0.0)
+	    xerr = INDEFR
+	else {
+	    if (IS_INDEFR(skysigma))
+	        constant = 0.0
+	    else
+	        constant = 4.0 * SQRTOFPI * sigma * skysigma ** 2 
+	    ratio = constant / xerr
+	    xerr = sigma ** 2 / (xerr * norm)
+	    xerr = sqrt (max (xerr, ratio * xerr))
+	    if (xerr > real (nx))
+		xerr = INDEFR
+	}
+
+	# Compute the y center and error.
+	yc = cy
+	nyiter = aptopt (Memr[ym], ny, yc, sigma, TOL, maxiter, YES)
+	if (yerr <= 0.0)
+	    yerr = INDEFR
+	else {
+	    if (IS_INDEFR(skysigma))
+	        constant = 0.0
+	    else
+	        constant = 4.0 * SQRTOFPI * sigma * skysigma ** 2 
+	    ratio = constant / yerr
+	    yerr = sigma ** 2 / (yerr * norm)
+	    yerr = sqrt (max (yerr, ratio * yerr))
+	    if (yerr > real (ny))
+		yerr = INDEFR
+	}
+
+	# Return appropriate error code.
+	call sfree (sp)
+	if (nxiter < 0 || nyiter < 0)
+	    return (AP_CTR_SINGULAR)
+	else if (nxiter > maxiter || nyiter > maxiter)
+	    return (AP_CTR_NOCONVERGE)
+	else
+	    return (AP_OK)
+end