1 files changed, 113 insertions, 0 deletions
diff --git a/noao/digiphot/apphot/fitpsf/apsfrefit.x b/noao/digiphot/apphot/fitpsf/apsfrefit.x
new file mode 100644
index 00000000..a1b950fa
--- /dev/null
+++ b/noao/digiphot/apphot/fitpsf/apsfrefit.x
@@ -0,0 +1,113 @@
+include <mach.h>
+include "../lib/apphotdef.h"
+include "../lib/apphot.h"
+include "../lib/fitpsfdef.h"
+include "../lib/fitpsf.h"
+include "../lib/noisedef.h"
+
+# APSFREFIT -- Procedure to fit a function to the data already in the
+# data buffers.
+
+int procedure apsfrefit (ap, im)
+
+pointer	ap		# pointer to the apphot structure
+pointer	im		# the input image descriptor
+
+int	ier
+pointer	psf, nse
+int	apsfradgauss(), apsfelgauss(), apsfmoments()
+real	datamin, datamax, dmin, dmax, threshold
+
+begin
+	psf = AP_PPSF(ap)
+	nse = AP_NOISE(ap)
+	call amovkr (INDEFR, Memr[AP_PPARS(psf)], AP_MAXNPARS(psf))
+	call amovkr (INDEFR, Memr[AP_PPERRS(psf)], AP_MAXNPARS(psf))
+
+	# Compute the minimum and maximum good data value.
+	if (IS_INDEFR(AP_DATAMIN(ap)))
+	    datamin = -MAX_REAL
+	else
+	    datamin = AP_DATAMIN(ap)
+	if (IS_INDEFR(AP_DATAMAX(ap)))
+	    datamax = MAX_REAL
+	else
+	    datamax = AP_DATAMAX(ap)
+
+	switch (AP_PSFUNCTION(psf)) {
+
+	case AP_RADGAUSS:
+
+	    ier = apsfradgauss (Memr[AP_PSFPIX(psf)], AP_PNX(psf), AP_PNY(psf),
+		AP_POSITIVE(ap), AP_FWHMPSF(ap) * AP_SCALE(ap), datamin,
+		datamax, AP_NOISEFUNCTION(nse), AP_EPADU(nse),
+		AP_READNOISE(nse) / AP_EPADU(nse), AP_PMAXITER(psf),
+		AP_PK2(psf), AP_PNREJECT(psf), Memr[AP_PPARS(psf)],
+		Memr[AP_PPERRS(psf)], AP_PSFNPARS(psf))
+
+	    Memr[AP_PPARS(psf)+1] = Memr[AP_PPARS(psf)+1] + AP_PFXCUR(psf) -
+		AP_PXC(psf)
+	    Memr[AP_PPARS(psf)+2] = Memr[AP_PPARS(psf)+2] + AP_PFYCUR(psf) -
+		AP_PYC(psf)
+
+	case AP_ELLGAUSS:
+
+	    ier = apsfelgauss (Memr[AP_PSFPIX(psf)], AP_PNX(psf), AP_PNY(psf),
+		AP_POSITIVE(ap), AP_FWHMPSF(ap) * AP_SCALE(ap), datamin,
+		datamax, AP_NOISEFUNCTION(nse), AP_EPADU(nse),
+		AP_READNOISE(nse) / AP_EPADU(nse), AP_PMAXITER(psf),
+		AP_PK2(psf), AP_PNREJECT(psf), Memr[AP_PPARS(psf)],
+		Memr[AP_PPERRS(psf)], AP_PSFNPARS(psf))
+
+	    Memr[AP_PPARS(psf)+1] = Memr[AP_PPARS(psf)+1] + AP_PFXCUR(psf) -
+		AP_PXC(psf)
+	    Memr[AP_PPARS(psf)+2] = Memr[AP_PPARS(psf)+2] + AP_PFYCUR(psf) -
+		AP_PYC(psf)
+
+	case AP_MOMENTS:
+
+	    call alimr (Memr[AP_PSFPIX(psf)], AP_PNX(psf) * AP_PNY(psf),
+		dmin, dmax)
+	    dmin = max (dmin, datamin)
+	    dmax = min (dmax, datamax)
+	    threshold = 0.0
+
+	    if (AP_POSITIVE(ap) == YES)
+	        ier = apsfmoments (Memr[AP_PSFPIX(psf)], AP_PNX(psf),
+		    AP_PNY(psf), dmin + threshold, dmax,
+		    AP_POSITIVE(ap), Memr[AP_PPARS(psf)],
+		    Memr[AP_PPERRS(psf)], AP_PSFNPARS(psf))
+	    else 
+	        ier = apsfmoments (Memr[AP_PSFPIX(psf)], AP_PNX(psf),
+		    AP_PNY(psf), dmax + threshold, dmin,
+		    AP_POSITIVE(ap), Memr[AP_PPARS(psf)],
+		    Memr[AP_PPERRS(psf)], AP_PSFNPARS(psf))
+
+	    Memr[AP_PPARS(psf)+1] = Memr[AP_PPARS(psf)+1] + AP_PFXCUR(psf) -
+		AP_PXC(psf)
+	    Memr[AP_PPARS(psf)+2] = Memr[AP_PPARS(psf)+2] + AP_PFYCUR(psf) -
+		AP_PYC(psf)
+
+	default:
+
+	    # do nothing gracefully
+        }
+
+        switch (AP_WCSOUT(ap)) {
+        case WCS_WORLD, WCS_PHYSICAL:
+            call ap_ltoo (ap, Memr[AP_PPARS(psf)+1], Memr[AP_PPARS(psf)+2],
+                Memr[AP_PPARS(psf)+1], Memr[AP_PPARS(psf)+2], 1)
+        case WCS_TV:
+            call ap_ltov (im, Memr[AP_PPARS(psf)+1], Memr[AP_PPARS(psf)+2],
+                Memr[AP_PPARS(psf)+1], Memr[AP_PPARS(psf)+2], 1)
+        default:
+            ;
+        }
+
+	if (ier == AP_NPSF_TOO_SMALL) {
+	    call amovkr (INDEFR, Memr[AP_PPARS(psf)], AP_PSFNPARS(psf))
+	    call amovkr (INDEFR, Memr[AP_PPERRS(psf)], AP_PSFNPARS(psf))
+	}
+
+	return (ier)
+end