Initial commit

1 files changed, 327 insertions, 0 deletions

diff --git a/noao/digiphot/daophot/allstar/dpastar.x b/noao/digiphot/daophot/allstar/dpastar.x
new file mode 100644
index 00000000..c8004f6d
--- /dev/null
+++ b/noao/digiphot/daophot/allstar/dpastar.x
@@ -0,0 +1,327 @@
+include <mach.h>
+include <imhdr.h>
+include "../lib/daophotdef.h"
+include "../lib/apseldef.h"
+include	"../lib/allstardef.h"
+
+
+# DP_ASTAR -- Begin doing the photometry.
+
+procedure dp_astar (dao, im, subim, allfd, rejfd, cache, savesub, version)
+
+pointer	dao		# pointer to the daophot structure
+pointer	im		# pointer to the input image
+pointer	subim		# pointer to the subtracted image
+int	allfd		# file descriptor for the output photometry file
+int	rejfd		# file descriptor for rejections files
+int	cache 		# cache the data ?
+int	savesub		# save the subtracted image ?
+int	version		# version number
+
+bool	clip
+int	nstar, row1_number, row2_number, niter, itsky, x1, x2, y1, y2, istar
+int	lstar, ldummy, newsky
+pointer	apsel, psffit, allstar, sp, indices, colpoints
+real	fradius, sepcrt, sepmin, clmpmx, wcrit, radius, xmin, xmax, ymin, ymax
+real	csharp, rdummy
+
+int	dp_alphot(), dp_alnstar()
+pointer	dp_gst(), dp_gwt(), dp_gdc()
+real	dp_usepsf()
+
+begin
+	# Define some daophot structure pointers.
+	apsel = DP_APSEL(dao)
+	psffit = DP_PSFFIT (dao)
+	allstar = DP_ALLSTAR (dao)
+
+	# Store the original fitting radius
+	fradius = DP_FITRAD(dao)
+	DP_FITRAD(dao) = min (DP_FITRAD(dao), DP_PSFRAD(dao))
+	DP_SFITRAD(dao) = DP_FITRAD(dao) * DP_SCALE(dao)
+
+	# Get some working memory.
+	call smark (sp)
+	call salloc (indices, NAPPAR, TY_INT)
+	call salloc (colpoints, ALL_NOUTCOLUMN, TY_INT)
+
+	# Define some daophot constants.
+	nstar = DP_APNUM(apsel)
+	csharp = 1.4427 * Memr[DP_PSFPARS(psffit)] *
+	    Memr[DP_PSFPARS(psffit)+1] / dp_usepsf (DP_PSFUNCTION(psffit), 0.0,
+	    0.0, DP_PSFHEIGHT(psffit), Memr[DP_PSFPARS(psffit)],
+	    Memr[DP_PSFLUT(psffit)], DP_PSFSIZE(psffit), DP_NVLTABLE(psffit),
+	    DP_NFEXTABLE(psffit), 0.0, 0.0, rdummy, rdummy)
+	if (IS_INDEFR(DP_MERGERAD(dao)))
+	    sepcrt = 2.0 * (Memr[DP_PSFPARS(psffit)] ** 2 +
+	        Memr[DP_PSFPARS(psffit)+1] ** 2)
+	else
+	    sepcrt = DP_MERGERAD(dao) ** 2
+	sepmin = min (1.0, FRACTION_MINSEP * sepcrt) 
+	itsky = 3
+
+	# Initialize the output photometry file for writing.
+	row1_number = 0
+	row2_number = 0
+	if (DP_TEXT(dao) == YES) {
+	    call dp_xnewal (dao, allfd)
+	    if (rejfd != NULL)
+	        call dp_xnewal (dao, rejfd)
+	} else {
+	    call dp_tnewal (dao, allfd, Memi[colpoints])
+	    if (rejfd != NULL)
+	        call dp_tnewal (dao, rejfd, Memi[colpoints])
+	}
+
+	# Allocate memory for the input, output and fitting arrays.
+	call dp_gnindices (Memi[indices])
+	call dp_rmemapsel (dao, Memi[indices], NAPPAR, nstar)
+	call dp_almemstar (dao, nstar, DP_MAXGROUP(dao))
+	call dp_cache (dao, im, subim, cache)
+
+	# Read in the input data and assign the initial weights.
+	call dp_setwt (dao, im)
+
+	# Convert the magnitudes to relative brightnesses.
+	call dp_alzero (dao, Memr[DP_APMAG(apsel)], nstar)
+
+	# Initialize all the fit/nofit indices and the error codes.
+	call amovki (NO, Memi[DP_ASKIP(allstar)], nstar)
+	call amovki (ALLERR_OK, Memi[DP_AIER(allstar)], nstar)
+
+	# Remove stars that have INDEF centers, are off the image altogether,
+	# or are too close to another star before beginning to fit.
+	call dp_strip (Memi[DP_APID(apsel)], Memr[DP_APXCEN(apsel)],
+	    Memr[DP_APYCEN(apsel)], Memr[DP_APMAG(apsel)],
+	    Memr[DP_APMSKY(apsel)], Memi[DP_ASKIP(allstar)],
+	    Memi[DP_AIER(allstar)], nstar, sepmin, int(IM_LEN(im,1)),
+	    int(IM_LEN(im,2)), DP_FITRAD(dao), DP_VERBOSE(dao))
+
+	# Write out results for the rejected stars.
+	call dp_wout (dao, im, Memr[DP_APXCEN(apsel)+nstar],
+	    Memr[DP_APYCEN(apsel)+nstar], Memr[DP_APXCEN(apsel)+nstar],
+	    Memr[DP_APYCEN(apsel)+nstar], DP_APNUM(apsel) - nstar)
+	if (DP_TEXT(dao) == YES) {
+	    call dp_xalwrite (allfd, rejfd, Memi[DP_APID(apsel)],
+		Memr[DP_APXCEN(apsel)], Memr[DP_APYCEN(apsel)],
+		Memr[DP_APMAG(apsel)], Memr[DP_APERR(apsel)],
+		Memr[DP_APMSKY(apsel)], Memr[DP_APCHI(apsel)],
+		Memr[DP_ANUMER(allstar)], Memr[DP_ADENOM(allstar)],
+		Memi[DP_ASKIP(allstar)], Memi[DP_AIER(allstar)], niter,
+		nstar + 1, DP_APNUM(apsel), DP_PSFMAG(psffit), csharp)
+	} else {
+	    call dp_talwrite (allfd, rejfd, Memi[colpoints],
+	        Memi[DP_APID(apsel)], Memr[DP_APXCEN(apsel)],
+	        Memr[DP_APYCEN(apsel)], Memr[DP_APMAG(apsel)],
+	        Memr[DP_APERR(apsel)], Memr[DP_APMSKY(apsel)],
+	        Memr[DP_APCHI(apsel)], Memr[DP_ANUMER(allstar)],
+	        Memr[DP_ADENOM(allstar)], Memi[DP_ASKIP(allstar)],
+		Memi[DP_AIER(allstar)], niter, nstar + 1, DP_APNUM(apsel),
+		row1_number, row2_number, DP_PSFMAG(psffit), csharp)
+	}
+
+	# Do some initialization for the fit.
+	clmpmx = CLAMP_FRACTION * DP_FITRAD(dao)
+	call aclrr (Memr[DP_AXOLD(allstar)], nstar)
+	call aclrr (Memr[DP_AYOLD(allstar)], nstar)
+	call amovkr (clmpmx, Memr[DP_AXCLAMP(allstar)], nstar)
+	call amovkr (clmpmx, Memr[DP_AYCLAMP(allstar)], nstar)
+	call amovkr (1.0, Memr[DP_ASUMWT(allstar)], nstar)
+
+	# Begin iterating.
+	for (niter = 1; (nstar > 0) && (niter <= DP_MAXITER (dao));
+	    niter = niter + 1) {
+
+	    if (DP_VERBOSE(dao) == YES) {
+	        call printf ("NITER = %d\n")
+		    call pargi (niter)
+	    }
+
+	    if ((DP_CLIPEXP(dao) != 0) && (niter >= 4))
+		clip = true
+	    else
+		clip = false
+
+	    # Define the critical errors.
+	    if (niter >= 15)
+	        wcrit = WCRIT15
+	    else if (niter >= 10)
+	        wcrit = WCRIT10
+	    else if (niter >= 5)
+	        wcrit = WCRIT5
+	    else
+	        wcrit = WCRIT0 
+
+	    # Sort the data on y.
+	    call dp_alsort (Memi[DP_APID(apsel)], Memr[DP_APXCEN(apsel)], 
+		Memr[DP_APYCEN(apsel)], Memr[DP_APMAG(apsel)], 
+		Memr[DP_APMSKY(apsel)], Memr[DP_ASUMWT(allstar)],
+		Memr[DP_AXOLD(allstar)], Memr[DP_AYOLD(allstar)],
+		Memr[DP_AXCLAMP(allstar)], Memr[DP_AYCLAMP(allstar)], nstar)
+
+	    # Compute the working radius. This is either the fitting radius
+	    # or the outer sky radius if this is an iteration during which
+	    # sky is to be determined whichever is larger.
+	    if ((DP_FITSKY(dao) == YES) && (mod (niter, itsky) == 0)) {
+		newsky = YES
+	        if ((DP_ANNULUS(dao) + DP_DANNULUS(dao)) > DP_FITRAD(dao))
+		    radius = DP_ANNULUS(dao) + DP_DANNULUS(dao)
+		else
+		    radius = DP_FITRAD(dao)
+	    } else {
+		newsky = NO
+		radius = DP_FITRAD(dao)
+	    }
+
+	    # Define region of the image used to fit the remaining stars.
+	    call alimr (Memr[DP_APXCEN(apsel)], nstar, xmin, xmax)
+	    call alimr (Memr[DP_APYCEN(apsel)], nstar, ymin, ymax)
+	    x1 = max (1, min (IM_LEN(im,1), int (xmin-radius)+1))
+	    x2 = max (1, min (IM_LEN(im,1), int (xmax+radius)))
+	    y1 = max (1, min (IM_LEN(im,2), int (ymin-radius)+1))
+	    y2 = max (1, min (IM_LEN (im,2), int (ymax+radius)))
+
+	    # Reinitialize the weight and scratch arrays / images .
+	    call dp_wstinit (dao, im, Memr[DP_APXCEN(apsel)],
+	        Memr[DP_APYCEN(apsel)], Memr[DP_APMAG(apsel)],
+		nstar, radius, x1, x2, y1, y2)
+
+	    # Recompute the initial sky estimates if that switch is enabled.
+	    if (newsky == YES)
+	        call dp_alsky (dao, im, Memr[DP_APXCEN(apsel)],
+		    Memr[DP_APYCEN(apsel)], Memr[DP_APMSKY(apsel)], nstar,
+		    x1, x2, y1, y2, DP_ANNULUS(dao), DP_ANNULUS(dao) +
+		    DP_DANNULUS(dao), 100, -MAX_REAL)
+
+	    # Group the remaining stars.
+	    call dp_regroup (Memi[DP_APID(apsel)], Memr[DP_APXCEN(apsel)], 
+		Memr[DP_APYCEN(apsel)], Memr[DP_APMAG(apsel)], 
+		Memr[DP_APMSKY(apsel)], Memr[DP_ASUMWT(allstar)],
+		Memr[DP_AXOLD(allstar)], Memr[DP_AYOLD(allstar)],
+		Memr[DP_AXCLAMP(allstar)], Memr[DP_AYCLAMP(allstar)], nstar,
+		DP_FITRAD(dao), Memi[DP_ALAST(allstar)])
+
+	    # Reset the error codes.
+	    call amovki (ALLERR_OK, Memi[DP_AIER(allstar)], nstar)
+
+	    # Do the serious fitting one group at a time.
+	    for (istar = 1; istar <= nstar; istar = lstar + 1) {
+
+		# Do the fitting a group at a time.
+	        lstar = dp_alphot (dao, im, nstar, istar, niter, sepcrt,
+		    sepmin, wcrit, clip, clmpmx, version) 
+
+		# Write the results.
+	        if (DP_TEXT(dao) == YES) {
+		    call dp_xalwrite (allfd, rejfd, Memi[DP_APID(apsel)],
+		        Memr[DP_APXCEN(apsel)], Memr[DP_APYCEN(apsel)],
+		        Memr[DP_APMAG(apsel)], Memr[DP_APERR(apsel)],
+		        Memr[DP_APMSKY(apsel)], Memr[DP_APCHI(apsel)],
+		        Memr[DP_ANUMER(allstar)], Memr[DP_ADENOM(allstar)],
+		        Memi[DP_ASKIP(allstar)], Memi[DP_AIER(allstar)],
+			niter, istar, lstar, DP_PSFMAG(psffit), csharp)
+	        } else {
+		    call dp_talwrite (allfd, rejfd, Memi[colpoints],
+		    Memi[DP_APID(apsel)], Memr[DP_APXCEN(apsel)],
+		    Memr[DP_APYCEN(apsel)], Memr[DP_APMAG(apsel)],
+		    Memr[DP_APERR(apsel)], Memr[DP_APMSKY(apsel)],
+		    Memr[DP_APCHI(apsel)], Memr[DP_ANUMER(allstar)],
+		    Memr[DP_ADENOM(allstar)], Memi[DP_ASKIP(allstar)],
+		    Memi[DP_AIER(allstar)], niter, istar, lstar, row1_number,
+		    row2_number, DP_PSFMAG(psffit), csharp)
+		}
+
+	    }
+
+	    # Find the last star in the list that still needs more work.
+	    nstar = dp_alnstar (Memi[DP_ASKIP(allstar)], nstar)
+
+	    # Remove the fitted stars from the list.
+	    for (istar = 1;  (istar < nstar) && nstar > 0; istar = istar + 1) {
+		call dp_alswitch (Memi[DP_APID(apsel)],
+		    Memr[DP_APXCEN(apsel)], Memr[DP_APYCEN(apsel)],
+		    Memr[DP_APMAG(apsel)], Memr[DP_APERR(apsel)],
+		    Memr[DP_APMSKY(apsel)], Memr[DP_ASUMWT(allstar)],
+		    Memr[DP_AXOLD(allstar)], Memr[DP_AYOLD(allstar)],
+		    Memr[DP_AXCLAMP(allstar)], Memr[DP_AYCLAMP(allstar)],
+		    Memi[DP_ASKIP(allstar)], istar, nstar)
+	    }
+
+	    # Flush the output buffers.
+	    if (dp_gwt (dao, im, ldummy, ldummy, READ_WRITE, YES) == NULL)
+		;
+	    if (dp_gst (dao, im, ldummy, ldummy, READ_ONLY, YES) == NULL)
+		;
+	    if (dp_gdc (dao, im, ldummy, ldummy, READ_WRITE, YES) == NULL)
+		;
+	}
+
+	# Release cached memory.
+	call dp_uncache (dao, subim, savesub)
+
+	call sfree (sp)
+
+	# Restore the original fitting radius
+	DP_FITRAD(dao) = fradius
+	DP_SFITRAD(dao) = DP_FITRAD(dao) * DP_SCALE(dao)
+end
+
+
+# DP_ALNSTAR -- Compute the number of stars left
+
+int procedure dp_alnstar (skip, nstar)
+
+int	skip[ARB]		# skip array
+int	nstar			# number of stars
+
+begin
+	while (nstar > 0) {
+	    if (skip[nstar] == NO)
+		break
+	    nstar = nstar - 1
+	}
+
+	return (nstar)
+end
+
+
+# DP_ALSWITCH -- Switch array elements.
+
+procedure dp_alswitch (id, xcen, ycen, mag, magerr, sky, chi, xold, yold,
+	xclamp, yclamp, skip, istar, nstar)
+
+int	id[ARB]			# array of ids
+real	xcen[ARB]		# array of x centers
+real	ycen[ARB]		# array of y centers
+real	mag[ARB]		# array of magnitudes
+real	magerr[ARB]		# array of magnitude errors
+real	sky[ARB]		# array of sky values
+real	chi[ARB]		# array of chi values
+real	xold[ARB]		# old x array
+real	yold[ARB]		# old y array
+real	xclamp[ARB]		# array of x clamps
+real	yclamp[ARB]		# array of y clamps
+int	skip[ARB]		# skip array
+int	istar			# next star
+int	nstar			# total number of stars
+
+int	dp_alnstar()
+
+begin
+	if (skip[istar] == YES) {
+	    id[istar] = id[nstar]
+	    xcen[istar] = xcen[nstar]
+	    ycen[istar] = ycen[nstar]
+	    mag[istar] = mag[nstar]
+	    magerr[istar] = magerr[nstar]
+	    sky[istar] = sky[nstar]
+	    chi[istar] = chi[nstar]
+	    xold[istar] = xold[nstar]
+	    yold[istar] = yold[nstar]
+	    xclamp[istar] = xclamp[nstar]
+	    yclamp[istar] = yclamp[nstar]
+	    skip[istar] = NO
+	    nstar = nstar - 1
+	    nstar = dp_alnstar (skip, nstar)
+	}
+end