1 files changed, 194 insertions, 0 deletions
diff --git a/noao/digiphot/daophot/psf/dpspstars.x b/noao/digiphot/daophot/psf/dpspstars.x
new file mode 100644
index 00000000..69b14933
--- /dev/null
+++ b/noao/digiphot/daophot/psf/dpspstars.x
@@ -0,0 +1,194 @@
+include <mach.h>
+include "../lib/daophotdef.h"
+include "../lib/apseldef.h"
+include "../lib/psfdef.h"
+
+define	NCOLUMN	5
+
+# DP_GPSTARS -- Select the psf stars.
+
+procedure dp_gpstars (dao, im, tp_in, tp_out, text_file, maxnpsf, gd, mgd, id,
+	mkstars, interactive, use_cmdfile)
+
+pointer	dao			# pointer to the daophot structure
+pointer	im			# pointer to the input image
+int	tp_in			# the input file descriptor
+int	tp_out			# the output file descriptor
+bool	text_file		# text or table file
+int	maxnpsf			# the maximum number of psf stars
+pointer	gd			# pointer to the graphics stream
+pointer	mgd			# pointer to the plot metacode file
+pointer	id			# pointer to the image display stream
+bool	mkstars			# mark the accepted and deleted stars
+bool	interactive		# interactive mode
+bool	use_cmdfile		# cursor command file mode
+
+int	ier, npsf
+pointer	apsel, sp, ocolpoint, index
+real	radius
+int	dp_pfstars(), dp_ipfstars()
+
+begin
+	# Get some daophot pointers.
+	apsel = DP_APSEL(dao)
+
+	# Define the radius for determining proximity.
+	radius = DP_PSFRAD(dao) + DP_FITRAD(dao) + 2.0
+
+	# Allocate some working memory.
+	call smark (sp)
+	call salloc (ocolpoint, NCOLUMN, TY_POINTER)
+	call salloc (index, DP_APNUM(apsel), TY_INT)
+
+	# Initialize the output file.
+	if (text_file) {
+	    call seek (tp_in, BOF)
+	    call dp_apheader (tp_in, tp_out)
+	    call dp_xpselpars (tp_out, DP_INIMAGE(dao), maxnpsf, DP_SCALE(dao),
+	        DP_SPSFRAD(dao), DP_SFITRAD(dao))
+	    call dp_xpbanner (tp_out)
+	} else {
+	    call dp_tpdefcol (tp_out, Memi[ocolpoint])
+	    call tbhcal (tp_in, tp_out)
+	    call dp_tpselpars (tp_out, DP_INIMAGE(dao), maxnpsf, DP_SCALE(dao),
+	        DP_SPSFRAD(dao), DP_SFITRAD(dao))
+	}
+
+	# Change all the INDEFS to a large negative number.
+	call dp_chindef (Memr[DP_APMAG(apsel)], Memr[DP_APMAG(apsel)],
+	    DP_APNUM(apsel), -MAX_REAL)
+
+	# Sort the stars.
+	call quick (Memr[DP_APMAG(apsel)], DP_APNUM(apsel), Memi[index], ier)
+	call dp_irectify (Memi[DP_APID(apsel)], Memi[index], DP_APNUM(apsel)) 
+	call dp_rectify (Memr[DP_APXCEN(apsel)], Memi[index], DP_APNUM(apsel)) 
+	call dp_rectify (Memr[DP_APYCEN(apsel)], Memi[index], DP_APNUM(apsel)) 
+	call dp_rectify (Memr[DP_APMSKY(apsel)], Memi[index], DP_APNUM(apsel)) 
+
+	# Select the stars and write them to the output file.
+	if (use_cmdfile)
+	    npsf = dp_ipfstars (dao, im, maxnpsf, -MAX_REAL, radius, gd,
+		mgd, id, mkstars, false, false)
+	else if (interactive)
+	    npsf = dp_ipfstars (dao, im, maxnpsf, -MAX_REAL, radius, gd,
+		mgd, id, mkstars, true, true)
+	else
+	    npsf = dp_pfstars (dao, im, Memi[DP_APID(apsel)],
+	        Memr[DP_APXCEN(apsel)], Memr[DP_APYCEN(apsel)],
+		Memr[DP_APMAG(apsel)], Memr[DP_APMSKY(apsel)],
+		DP_APNUM(apsel), maxnpsf, -MAX_REAL, radius, mgd)
+
+	if (DP_VERBOSE(dao) == YES) {
+	    call printf ("\nTotal of %d PSF stars selected\n")
+		call pargi (npsf)
+	}
+
+	# Write out the stars.
+	call dp_wout (dao, im, Memr[DP_APXCEN(apsel)], Memr[DP_APYCEN(apsel)],
+	    Memr[DP_APXCEN(apsel)], Memr[DP_APYCEN(apsel)], npsf)
+	call dp_wpstars (tp_out, Memi[ocolpoint], text_file,
+	    Memi[DP_APID(apsel)], Memr[DP_APXCEN(apsel)],
+	    Memr[DP_APYCEN(apsel)], Memr[DP_APMAG(apsel)],
+	    Memr[DP_APMSKY(apsel)], npsf)
+
+	# Free memory.
+	call sfree (sp)
+end
+
+
+# DP_CHINDEF -- Change all the INDEFS to a legal number which is a lower limit.
+
+procedure dp_chindef (a, b, npix, lolimit)
+
+real	a[ARB]			# the input array.
+real	b[ARB]			# the output array.
+int	npix			# number of points
+real	lolimit			# the lower limit
+
+int	i
+
+begin
+	do i = 1, npix {
+	    if (IS_INDEFR(a[i]))
+		b[i] = lolimit
+	    else
+		b[i] = a[i]
+	}
+end
+
+
+# DP_PFSTARS -- Select the psf stars.
+
+int procedure dp_pfstars (dao, im, ids, xcen, ycen, mag, sky, nstar, maxnpsf,
+	lolimit, radius, mgd) 
+
+pointer	dao			# pointer to the daophot structure
+pointer	im			# pointer to the input image
+int	ids[ARB]		# array of star ids
+real	xcen[ARB]		# array of x coordinates
+real	ycen[ARB]		# array of y coordinates
+real	mag[ARB]		# array of magnitudes
+real	sky[ARB]		# array of sky values
+int	nstar			# the number of stars
+int	maxnpsf			# the maximum number of psf stars
+real	lolimit			# lower data limit
+real	radius			# minimum separation
+pointer	mgd			# pointer to the plot metacode file
+
+bool	omit
+int	istar, jstar, npsf
+real	radsq, dy2, dr2
+int	dp_addstar()
+
+begin
+	# Initialize the list.
+	call dp_pseti (dao, PNUM, 0)
+
+	# Set the radius.
+	radsq = radius * radius
+
+	# Get the first star.
+	if ((mag[1] > lolimit) && (dp_addstar (dao, im, xcen[1], ycen[1],
+	    INDEFR, ids[1], NULL, mgd, false)) == OK) {
+	    npsf = 1
+	} else
+	    npsf = 0
+
+	# Loop over the candidate psf stars.
+	do istar = 2, nstar {
+
+	    # Test for the maximum number of psf stars.
+	    if (npsf >= maxnpsf)
+		break
+
+	    # Test that the candidate psf stars are not saturated and
+	    # are sufficiently far from the edge of the frame.
+
+	    if (mag[istar] <= lolimit)
+		next
+
+	    # Text that there are no brighter stars with a distance squared
+	    # of radsq.
+	    omit = false
+	    do jstar = 1, istar - 1 {
+		dy2 = abs (ycen[jstar] - ycen[istar])
+		if (dy2 >= radius)
+		    next
+		dr2 = (xcen[jstar] - xcen[istar]) ** 2 + dy2 ** 2
+		if (dr2 >= radsq)
+		    next
+		omit = true
+		break
+	    }
+
+	    if (omit)
+		next
+
+	    if (dp_addstar (dao, im, xcen[istar], ycen[istar], INDEFR,
+	        ids[istar], NULL, mgd, false) == ERR)
+		next
+	    npsf = npsf + 1
+	}
+
+	return (npsf)
+end