Initial commit

1 files changed, 277 insertions, 0 deletions

diff --git a/noao/digiphot/daophot/peak/dppeakphot.x b/noao/digiphot/daophot/peak/dppeakphot.x
new file mode 100644
index 00000000..182b10d8
--- /dev/null
+++ b/noao/digiphot/daophot/peak/dppeakphot.x
@@ -0,0 +1,277 @@
+include	<imhdr.h>
+include <tbset.h>
+include "../lib/daophotdef.h"
+include "../lib/apseldef.h"
+include "../lib/peakdef.h"
+
+
+# DP_PEAKPHOT -- Fit the PSF to a single star.
+
+procedure dp_peakphot (dao, im, tp, tpout, tprej, ap_text)
+
+pointer	dao			# pointer to the daophot structure
+pointer	im			# the input image descriptor
+int	tp			# the input photometry file descriptor
+int	tpout			# the ouput photometry file descriptor
+int	tprej			# the rejections file descriptor
+bool	ap_text			# which style of photometry
+
+real	rel_bright, xold, yold, x, y, dx, dy, mag, sky, errmag
+real	chi, sharp, radius, itx, ity, otx, oty
+pointer	psffit, key, sp, subim, colpoint, indices, fields, perror
+int	id, in_nrow, instar, lowx, lowy, nxpix, nypix, niter, out_nrow
+int	rout_nrow, nterm, ier, plen
+
+int	tbpsta(), dp_rrphot(), dp_pkfit(), dp_gpkerr()
+pointer	dp_gsubrast()
+
+begin
+	# Get the daophot pointers.
+	psffit = DP_PSFFIT (dao)
+
+	# Store the original fitting radius.
+	radius = DP_FITRAD(dao)
+
+	# Check that the fitting radius is less than the psf radius.
+	DP_FITRAD(dao) = min (DP_FITRAD(dao), DP_PSFRAD(dao))
+	DP_SFITRAD(dao) = DP_FITRAD(dao) * DP_SCALE(dao)
+
+	# Allocate working space.
+	call smark (sp)
+	call salloc (indices, NAPPAR, TY_INT)
+	call salloc (fields, SZ_LINE, TY_CHAR)
+	call salloc (colpoint, PK_NOUTCOL, TY_INT)
+	call salloc (perror, SZ_FNAME, TY_CHAR)
+
+	# Initialze the output table.
+	if (DP_TEXT(dao) == YES) {
+	    call dp_xnewpeak (dao, tpout)
+	    if (tprej != NULL)
+	        call dp_xnewpeak (dao, tprej)
+	} else {
+	    call dp_tnewpeak (dao, tpout, Memi[colpoint])
+	    if (tprej != NULL)
+	        call dp_tnewpeak (dao, tprej, Memi[colpoint])
+	}
+
+	# Intialize the input table.
+	if (ap_text) {
+	    call pt_kyinit (key)
+	    Memi[indices] = DP_PAPID
+	    Memi[indices+1] = DP_PAPXCEN
+	    Memi[indices+2] = DP_PAPYCEN
+	    Memi[indices+3] = DP_PAPMAG1
+	    Memi[indices+4] = DP_PAPSKY
+	    call dp_gappsf (Memi[indices], Memc[fields], NAPRESULT)
+	    in_nrow = 0
+	} else {
+	    call dp_tpkinit (tp, Memi[indices])
+	    in_nrow = tbpsta (tp, TBL_NROWS)
+	}
+
+	# Initialize the photometry file reading code.
+	instar = 0
+
+	# Initialize the fitting code.
+	if (DP_RECENTER(dao) == YES)
+	    nterm = 3
+	else
+	    nterm = 1
+	if (DP_FITSKY(dao) == YES)
+	    nterm = nterm + 1
+	call dp_mempk (dao, nterm)
+
+	out_nrow = 0
+	rout_nrow = 0
+	repeat {
+
+	    # Read in the photometry for a single star.
+	    if (dp_rrphot (tp, key, Memc[fields], Memi[indices], id, itx,
+	        ity, sky, mag, instar, in_nrow) == EOF)
+		break
+
+	    # Convert to and from logical coordinates.
+	    call dp_win (dao, im, itx, ity, x, y, 1)
+	    call dp_wout (dao, im, x, y, otx, oty, 1)
+
+	    if (DP_VERBOSE(dao) == YES) {
+	        call printf (
+		    "Star: %5d X: %8.2f Y: %8.2f Mag: %8.2f Sky: %8.2f\n")
+		    call pargi (id)
+		    call pargr (otx)
+		    call pargr (oty)
+		    call pargr (mag)
+		    call pargr (sky)
+	    }
+
+	    # Check that the center is defined.
+	    if (IS_INDEFR(x) || IS_INDEFR(y)) {
+		if (DP_VERBOSE(dao) == YES) {
+		    call printf (
+		        "\tWarning: X and/or Y for star %d are undefined\n")
+		        call pargi (id)
+		}
+		next
+	    }
+
+    	    # Read in the subraster.
+	    subim = dp_gsubrast (im, x, y, DP_FITRAD(dao), lowx, lowy,
+	        nxpix, nypix)
+	    if (subim == NULL) {
+		if (DP_VERBOSE(dao) == YES) {
+		    call printf (
+		        "\tWarning: Cannot read in image data for star %d\n")
+		        call pargi (id)
+		}
+		next
+	    }
+
+	    # Save the old x and y values for use with the variable psf
+	    # option.
+	    xold = x
+	    yold = y
+	    call dp_wpsf (dao, im, xold, yold, xold, yold, 1)
+
+	    # Compute the relative centers and the relative brightness and
+	    # fit the star.
+	    if (IS_INDEFR(sky)) {
+
+		ier = PKERR_INDEFSKY
+
+	    } else {
+
+	        x = x - lowx + 1.0
+	        y = y - lowy + 1.0
+	        dx = (xold - 1.0) / DP_PSFX(psffit) - 1.0
+	        dy = (yold - 1.0) / DP_PSFY(psffit) - 1.0
+	        if (IS_INDEFR(mag))
+	            mag = DP_PSFMAG (psffit) + DELTA_MAG
+	        rel_bright = DAO_RELBRIGHT (psffit, mag)
+	        ier = dp_pkfit (dao, Memr[subim], nxpix, nypix, DP_FITRAD(dao),
+		    x, y, dx, dy, rel_bright, sky, errmag, chi, sharp, niter)
+	        x = x + lowx - 1.0
+	        y = y + lowy - 1.0
+	    }
+
+	    call dp_wout (dao, im, x, y, otx, oty, 1)
+
+	    if (ier != PKERR_OK) {
+
+		# Set fitting parameters to INDEF.
+		mag = INDEFR
+		niter = 0
+		errmag = INDEFR
+		chi = INDEFR
+		sharp = INDEFR
+
+		if (DP_VERBOSE(dao) == YES) {
+		    switch (ier) {
+		    case PKERR_INDEFSKY:
+		        call printf (
+		        "\tWarning: The sky value for star %d is undefined\n")
+		            call pargi (id)
+		    case PKERR_NOPIX:
+		        call printf (
+		            "\tWarning: Too few pixels to fit star %d\n")
+		            call pargi (id)
+		    case PKERR_SINGULAR:
+		        call printf (
+		            "\tWarning: Singular matrix computed for star %d\n")
+		            call pargi (id)
+		    case PKERR_FAINT:
+		        call printf ("\tWarning: Star %d too faint\n")
+		            call pargi (id)
+		    default:
+		        call printf (
+		           "\tWarning: Unknown error detected for star %d\n")
+		            call pargi (id)
+		    }
+		}
+
+	    } else {
+
+	        # Compute the results.
+	        mag = DP_PSFMAG (psffit) - 2.5 * log10 (rel_bright)
+	        errmag = 1.085736 * errmag / rel_bright
+		if (errmag >= 2.0)
+		    errmag = INDEFR
+
+	    	if (DP_VERBOSE (dao) == YES) {
+             	        call printf (
+		  "\tFIT:  Star: %5d X: %8.2f Y: %8.2f Mag: %8.2f Sky =%8.2f\n")
+	    	            call pargi (id)
+		    	    call pargr (otx)
+		    	    call pargr (oty)
+		   	    call pargr (mag)
+		    	    call pargr (sky)
+		}
+	    }
+
+	    # Get the error code.
+	    plen = dp_gpkerr (ier, Memc[perror], SZ_FNAME)
+
+	    # Now write the results to the output photometry or rejections
+	    # file.
+	    if (DP_TEXT(dao) == YES) {
+		if ((tprej != NULL) && (ier != PKERR_OK))
+		    call dp_xpkwrite (tprej, id, otx, oty, mag, errmag, sky,
+		        niter, chi, sharp, ier, Memc[perror], plen)
+		else
+		    call dp_xpkwrite (tpout, id, otx, oty, mag, errmag, sky,
+		        niter, chi, sharp, ier, Memc[perror], plen)
+	    } else {
+		if ((tprej != NULL) && (ier != PKERR_OK)) {
+	            rout_nrow = rout_nrow + 1
+		    call dp_tpkwrite (tprej, Memi[colpoint], id, otx, oty, mag,
+		        errmag, sky, niter, chi, sharp, ier,
+			Memc[perror], plen, rout_nrow)
+		} else {
+	            out_nrow = out_nrow + 1
+		    call dp_tpkwrite (tpout, Memi[colpoint], id, otx, oty, mag,
+		        errmag, sky, niter, chi, sharp, ier, Memc[perror],
+			plen, out_nrow)
+		}
+	    }
+	}
+
+	# Free the text file descriptor.
+	if (ap_text)
+	    call pt_kyfree (key)
+
+	# Restore the original fitting radius.
+	DP_FITRAD(dao) = radius
+	DP_SFITRAD(dao) = DP_FITRAD(dao) * DP_SCALE(dao)
+
+	call sfree (sp)
+end
+
+
+# DP_GPKERR -- Set the PEAK task error code string.
+
+int procedure dp_gpkerr (ier, perror, maxch)
+
+int	ier		# the integer error code
+char	perror		# the output error code string
+int	maxch		# the maximum size of the error code string
+
+int	plen
+int	gstrcpy()
+
+begin
+	switch (ier) {
+	case PKERR_OK:
+	    plen = gstrcpy ("No_error", perror, maxch)
+	case PKERR_INDEFSKY:
+	    plen = gstrcpy ("Bad_sky", perror, maxch)
+	case PKERR_NOPIX:
+	    plen = gstrcpy ("Npix_too_few", perror, maxch)
+	case PKERR_SINGULAR:
+	    plen = gstrcpy ("Singular", perror, maxch)
+	case PKERR_FAINT:
+	    plen = gstrcpy ("Too_faint", perror, maxch)
+	default:
+	    plen = gstrcpy ("No_error", perror, maxch)
+	}
+
+	return (plen)
+end