Initial commit

8 files changed, 1569 insertions, 0 deletions

diff --git a/noao/digiphot/daophot/peak/dpmempk.x b/noao/digiphot/daophot/peak/dpmempk.x
new file mode 100644
index 00000000..97c7ff48
--- /dev/null
+++ b/noao/digiphot/daophot/peak/dpmempk.x
@@ -0,0 +1,72 @@
+include "../lib/daophotdef.h"
+include "../lib/peakdef.h"
+
+
+# DP_PKSETUP -- Initialize the PEAK fitting structure.
+
+procedure dp_pksetup (dao)
+
+pointer	dao			# pointer to the daophot structure
+
+pointer	peak
+
+begin
+	call malloc (DP_PEAK(dao), LEN_PKSTRUCT, TY_STRUCT)
+	peak = DP_PEAK(dao)
+
+	DP_PKNTERM(peak) = 0
+	DP_PKCLAMP(peak) = NULL
+	DP_PKNORMAL(peak) = NULL
+	DP_PKRESID(peak) = NULL
+	DP_PKDERIV(peak) = NULL
+	DP_PKRESULT(peak) = NULL
+	DP_PKOLDRESULT(peak) = NULL
+end
+
+
+# DP_MEMPK -- Allocate memory for the PEAK fitting arrays.
+
+procedure dp_mempk (dao, nterm)
+
+pointer	dao			# pointer to the daophot structure
+int	nterm			# the number of terms to be fit
+
+pointer	peak
+
+begin
+	peak = DP_PEAK(dao)
+
+	call malloc (DP_PKCLAMP(peak), nterm, TY_REAL)
+	call malloc (DP_PKNORMAL(peak), nterm * nterm, TY_REAL)
+	call malloc (DP_PKRESID(peak), nterm * nterm, TY_REAL)
+	call malloc (DP_PKDERIV(peak), nterm * nterm, TY_REAL)
+	call malloc (DP_PKRESULT(peak), nterm * nterm, TY_REAL)
+	call malloc (DP_PKOLDRESULT(peak), nterm * nterm, TY_REAL)
+	DP_PKNTERM(peak) = nterm
+end
+
+
+# DP_PKCLOSE -- Free the PEAK fitting structure.
+
+procedure dp_pkclose (dao)
+
+pointer	dao			# pointer to the daophot structure
+
+pointer	peak
+
+begin
+	peak = DP_PEAK(dao)
+	if (DP_PKCLAMP(peak) != NULL)
+	    call mfree (DP_PKCLAMP(peak), TY_REAL)
+	if (DP_PKNORMAL(peak) != NULL)
+	    call mfree (DP_PKNORMAL(peak), TY_REAL)
+	if (DP_PKRESID(peak) != NULL)
+	    call mfree (DP_PKRESID(peak), TY_REAL)
+	if (DP_PKDERIV(peak) != NULL)
+	    call mfree (DP_PKDERIV(peak), TY_REAL)
+	if (DP_PKRESULT(peak) != NULL)
+	    call mfree (DP_PKRESULT(peak), TY_REAL)
+	if (DP_PKOLDRESULT(peak) != NULL)
+	    call mfree (DP_PKOLDRESULT(peak), TY_REAL)
+	call mfree (peak, TY_STRUCT)
+end
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
diff --git a/noao/digiphot/daophot/peak/dppkconfirm.x b/noao/digiphot/daophot/peak/dppkconfirm.x
new file mode 100644
index 00000000..c7dc6333
--- /dev/null
+++ b/noao/digiphot/daophot/peak/dppkconfirm.x
@@ -0,0 +1,25 @@
+# DP_PKCONFIRM -- Procedure to confirm the critical peak parameters.
+
+procedure dp_pkconfirm (dao)
+
+pointer	dao		# pointer to the daophot structure
+
+begin
+	call printf ("\n")
+
+	# Confirm recentering and sky fitting.
+	call dp_vrecenter (dao)
+	call dp_vfitsky (dao)
+
+	# Confirm the psf radius.
+	call dp_vpsfrad (dao)
+
+	# Confirm the fitting radius.
+	call dp_vfitrad (dao)
+
+	# Confirm the minimum and maximum good data values.
+	call dp_vdatamin (dao)
+	call dp_vdatamax (dao)
+
+	call printf ("\n")
+end
diff --git a/noao/digiphot/daophot/peak/dppkfit.x b/noao/digiphot/daophot/peak/dppkfit.x
new file mode 100644
index 00000000..7818100c
--- /dev/null
+++ b/noao/digiphot/daophot/peak/dppkfit.x
@@ -0,0 +1,411 @@
+include <mach.h>
+include	"../lib/daophotdef.h"
+include	"../lib/peakdef.h"
+
+# DP_PKFIT -- Do the actual fit.
+
+int procedure dp_pkfit (dao, subim, nx, ny, radius, x, y, dx, dy, rel_bright,
+	sky, errmag, chi, sharp, iter)
+
+pointer	dao		# pointer to the DAOPHOT Structure
+real	subim[nx,ny]	# pointer to the image subraster
+int	nx, ny		# size of the image subraster
+real	radius		# the fitting radius
+real	x, y		# initial estimate of the stellar position
+real	dx, dy		# distance of star from the psf position
+real	rel_bright	# initial estimate of stellar brightness
+real	sky		# the sky value associated with this star
+real	errmag		# error estimate for this star
+real	chi		# estimated goodness of fit parameter
+real	sharp		# broadness of the profile compared to the PSF
+int	iter		# number of iterations needed for a fit.
+
+bool	clip, redo
+int	i, flag, npix
+pointer	psffit, peak
+real	ronoise, numer, denom, chiold, sum_weight, noise, wcrit
+int	dp_fitbuild()
+
+begin
+	# Get the pointer to the PSF structure.
+	psffit = DP_PSFFIT (dao)
+	peak = DP_PEAK(dao)
+
+	# Initialize the parameters which control the fit.
+
+	chiold = 1.0
+	sharp = INDEFR
+	clip = false
+	ronoise = (DP_READNOISE(dao) / DP_PHOTADU(dao)) ** 2
+
+	call amovkr (1.0, Memr[DP_PKCLAMP(peak)], DP_PKNTERM(peak))
+	call amovkr (0.0, Memr[DP_PKOLDRESULT(peak)], DP_PKNTERM(peak))
+
+	# Iterate until a solution is found.
+
+	for (iter = 1; iter <= DP_MAXITER(dao); iter = iter + 1) {
+
+	    # Initialize the matrices and vectors required by the fit.
+
+	    chi = 0.0
+	    numer = 0.0
+	    denom = 0.0
+	    sum_weight = 0.0
+	    call aclrr (Memr[DP_PKRESID(peak)], DP_PKNTERM(peak))
+	    call aclrr (Memr[DP_PKNORMAL(peak)], DP_PKNTERM(peak) *
+	        DP_PKNTERM(peak))
+
+	    # Set up the critical error limit.
+	    if (iter >= WCRIT_NMAX)
+		wcrit = WCRIT_MAX
+	    else if (iter >= WCRIT_NMED)
+		wcrit = WCRIT_MED
+	    else if (iter >= WCRIT_NMIN)
+		wcrit = WCRIT_MIN
+	    else
+		wcrit = MAX_REAL
+
+
+	    # Build up the vector of residuals and the normal matrix.
+
+	    npix = dp_fitbuild (dao, subim, nx, ny, radius, x, y, dx, dy,
+	        rel_bright, sky, chiold, chi, clip, iter,
+		Memr[DP_PKCLAMP(peak)], Memr[DP_PKNORMAL(peak)],
+		Memr[DP_PKRESID(peak)], Memr[DP_PKDERIV(peak)],
+		DP_PKNTERM[(peak)], numer, denom, sum_weight)
+
+	    # Return if the iteration was unsuccessful.
+
+	    if (npix < MIN_NPIX)
+		return (PKERR_NOPIX)
+
+	    # Invert the matrix. Return if inversion was unsuccessful or
+	    # if any of the diagonal elements are less or equal to 0.0.
+
+	    call invers (Memr[DP_PKNORMAL(peak)], DP_PKNTERM(peak),
+	        DP_PKNTERM(peak), flag)
+	    if (flag != 0)
+		return (PKERR_SINGULAR)
+	    do i = 1, DP_PKNTERM(peak) {
+		if (Memr[DP_PKNORMAL(peak)+(i-1)*DP_PKNTERM(peak)+i-1] <= 0.0)
+		    return (PKERR_SINGULAR)
+	    }
+	    # Solve the matrix.
+
+	    call mvmul (Memr[DP_PKNORMAL(peak)], DP_PKNTERM(peak),
+	        DP_PKNTERM(peak), Memr[DP_PKRESID(peak)],
+		Memr[DP_PKRESULT(peak)])
+
+	    # In the beginning the brightness of the star will not be
+	    # permitted to change by more than two magnitudes per
+	    # iteration (that is to say if the estimate is getting
+	    # brighter, it may not get brighter by more than 525%
+	    # per iteration, and if it is getting fainter, it may not
+	    # get fainter by more than 84% per iteration). The x and y
+	    # coordinates of the centroid will be allowed to change by
+	    # no more than one-half pixel per iteration. Any time
+	    # that a parameter correction changes sign, the maximum
+	    # permissible change in that parameter will be reduced
+	    # by a factor of 2.
+
+	    # Perform the sign check.
+
+	    do i = 1, DP_PKNTERM(peak) {
+		if ((Memr[DP_PKOLDRESULT(peak)+i-1] *
+		    Memr[DP_PKRESULT(peak)+i-1]) < 0.0) 
+		    Memr[DP_PKCLAMP(peak)+i-1] = 0.5 *
+		        Memr[DP_PKCLAMP(peak)+i-1]
+		else
+		    Memr[DP_PKCLAMP(peak)+i-1] = min (1.0, 1.1 *
+		        Memr[DP_PKCLAMP(peak)+i-1])
+		Memr[DP_PKOLDRESULT(peak)+i-1] = Memr[DP_PKRESULT(peak)+i-1]
+	    }
+
+	    # Compute the new x, y, sky, and magnitude.
+	    rel_bright = rel_bright + Memr[DP_PKRESULT(peak)] /
+	        (1.0 + max (Memr[DP_PKRESULT(peak)] /
+	        (MAX_DELTA_BRIGHTER * rel_bright), -Memr[DP_PKRESULT(peak)] /
+		(MAX_DELTA_FAINTER * rel_bright)) / Memr[DP_PKCLAMP(peak)])
+
+	    # Return if the star becomes too faint)
+	    if (rel_bright < MIN_REL_BRIGHT)
+		return (PKERR_FAINT)
+
+	    if (DP_RECENTER(dao) == YES) {
+	        x = x + Memr[DP_PKRESULT(peak)+1] / (1.0 +
+	            abs (Memr[DP_PKRESULT(peak)+1]) / (MAX_DELTA_PIX *
+		    Memr[DP_PKCLAMP(peak)+1]))
+	        y = y + Memr[DP_PKRESULT(peak)+2] / (1.0 +
+	            abs (Memr[DP_PKRESULT(peak)+2]) / (MAX_DELTA_PIX *
+		    Memr[DP_PKCLAMP(peak)+2]))
+	    }
+
+	    if (DP_FITSKY(dao) == YES) {
+		noise = sqrt ((abs (sky / DP_PHOTADU(dao)) + ronoise)) 
+	        sky = sky + max (-3.0 * noise, min (Memr[DP_PKRESULT(peak)+
+		    DP_PKNTERM(peak)-1], 3.0 * noise)) 
+	    }
+
+	    # Force at least one iteration before checking for convergence.
+	    if (iter <= 1)
+		next
+
+	    # Start on the assumption the fit has converged.
+
+	    redo = false
+
+	    # Check for convergence. If the most recent computed correction 
+	    # to the brightness is larger than 0.1% of the brightness or
+	    # 0.05 * sigma of the brightness whichever is larger, convergence
+	    # has not been achieved.
+
+	    errmag = chiold * sqrt (Memr[DP_PKNORMAL(peak)])
+	    if (clip) {
+	        if (abs (Memr[DP_PKRESULT(peak)]) > max ((MAX_NEW_ERRMAG *
+		    errmag), (MAX_NEW_RELBRIGHT1 * rel_bright))) {
+		    redo = true
+	        } else {
+		    if (DP_RECENTER(dao) == YES) {
+		        if (max (abs (Memr[DP_PKRESULT(peak)+1]),
+	                    abs (Memr[DP_PKRESULT(peak)+2])) > MAX_PIXERR1)
+		            redo = true
+		    }
+		    if (DP_FITSKY(dao) == YES) {
+		        if (abs (Memr[DP_PKRESULT(peak)+DP_PKNTERM(peak)-1]) >
+			    1.0e-4 * sky)
+			    redo = true
+		    }
+		}
+	    } else {
+	        if (abs (Memr[DP_PKRESULT(peak)]) > max (errmag,
+		    (MAX_NEW_RELBRIGHT2 * rel_bright))) {
+		    redo = true
+	        } else {
+		    if (DP_RECENTER(dao) == YES) {
+		        if (max (abs (Memr[DP_PKRESULT(peak)+1]),
+	                    abs (Memr[DP_PKRESULT(peak)+2])) > MAX_PIXERR2)
+		            redo = true
+		    }
+		    if (DP_FITSKY(dao) == YES) {
+		        if (abs (Memr[DP_PKRESULT(peak)+DP_PKNTERM(peak)-1]) >
+			    1.0e-4 * sky)
+			    redo = true
+		    }
+		}
+	    }
+	    if (redo)
+		next
+
+	    if ((iter < DP_MAXITER(dao)) && (! clip)) {
+		if (DP_CLIPEXP(dao) > 0)
+		    clip = true
+		call aclrr (Memr[DP_PKOLDRESULT(peak)], DP_PKNTERM(peak))
+		call amaxkr (Memr[DP_PKCLAMP(peak)], MAX_CLAMPFACTOR,
+		    Memr[DP_PKCLAMP(peak)], DP_PKNTERM(peak))
+	    } else {
+		sharp = 1.4427 * Memr[DP_PSFPARS(psffit)] *
+	    	    Memr[DP_PSFPARS(psffit)+1] * numer / (DP_PSFHEIGHT(psffit) *
+	    	    rel_bright * denom)
+		if (sharp <= MIN_SHARP || sharp >= MAX_SHARP)
+	    	    sharp = INDEFR
+		break
+	    }
+	}
+
+	if (iter > DP_MAXITER(dao))
+	    iter = iter - 1
+	if ((errmag / rel_bright) >= wcrit)
+	    return (PKERR_FAINT)
+	else
+	    return (PKERR_OK)
+end
+
+
+# DP_FITBUILD -- Build the normal and vector of residuals for the fit.
+
+int procedure dp_fitbuild (dao, subim, nx, ny, radius, x, y, xfrom_psf,
+	yfrom_psf, rel_bright, sky, chiold, chi, clip, iter, clamp, normal,
+	resid, deriv, nterm, numer, denom, sum_weight)
+
+pointer	dao			# pointer to the DAOPHOT Structure
+real	subim[nx,ny]		# subimage containing star
+int	nx, ny			# size of the image subraster
+real	radius			# the fitting radius
+real	x, y			# initial estimate of the position
+real	xfrom_psf, yfrom_psf	# distance from the psf star
+real	rel_bright		# initial estimate of stellar brightness
+real	sky			# the sky value associated with this star
+real	chiold			# previous estimate of gof
+real	chi			# estimated goodness of fit parameter
+bool	clip			# clip the weights ?
+int	iter			# current iteration number
+real	clamp[ARB]		# clamp array
+real	normal[nterm,ARB]	# normal matrix
+real	resid[ARB]		# residual matrix
+real	deriv[ARB]		# derivative matrix
+int	nterm			# the number of terms to be fit
+real	numer, denom		# used in sharpness calculation
+real	sum_weight		# sum of the weights
+
+int	i, j, ix, iy, lowx, lowy, highx, highy, npix
+pointer	psffit
+real	fitradsq, pererr, peakerr, datamin, datamax, read_noise
+real	dx, dy, dxsq, dysq, radsq, dvdx, dvdy, d_pixval
+real	pred_pixval, sigma, sigmasq, relerr, weight
+real 	rhosq, pixval, dfdsig
+
+real	dp_usepsf()
+
+begin
+	# Get the pointer to the PSF structure.
+	psffit = DP_PSFFIT (dao)
+
+	# Set up some constants.
+	fitradsq = radius * radius
+	read_noise = (DP_READNOISE(dao) / DP_PHOTADU(dao)) ** 2
+	pererr = 0.01 * DP_FLATERR(dao)
+	peakerr = 0.01 * DP_PROFERR(dao) /
+	    (Memr[DP_PSFPARS(psffit)] * Memr[DP_PSFPARS(psffit)+1])
+	if (IS_INDEFR (DP_MINGDATA(dao)))
+	    datamin = -MAX_REAL
+	else
+	    datamin = DP_MINGDATA(dao)
+	if (IS_INDEFR (DP_MAXGDATA(dao)))
+	    datamax = MAX_REAL
+	else
+	    datamax = DP_MAXGDATA(dao)
+
+	# Define the size of the subraster to be used in the fit.
+	lowx = max (1, int (x - radius))
+	lowy = max (1, int (y - radius))
+	highx = min (nx, int (x + radius) + 1)	
+	highy = min (ny, int (y + radius) + 1)	
+
+	npix = 0
+	do iy = lowy, highy {
+
+	    dy = real (iy) - y
+	    dysq = dy * dy
+
+	    do ix = lowx, highx {
+
+		# Is this pixel in the good data range.
+		pixval = subim[ix,iy]
+		if (pixval < datamin || pixval > datamax)
+		    next
+
+		# Is this pixel inside the fitting radius?
+		dx = real(ix) - x
+		dxsq = dx * dx
+		radsq = (dxsq + dysq) / fitradsq
+		if ((1.0 - radsq)  <= PEAK_EPSILONR)
+		    next
+
+		# We have a good pixel within the fitting radius.
+		deriv[1] = dp_usepsf (DP_PSFUNCTION(psffit), dx, dy,
+		    DP_PSFHEIGHT(psffit), Memr[DP_PSFPARS(psffit)],
+		    Memr[DP_PSFLUT(psffit)], DP_PSFSIZE(psffit),
+		    DP_NVLTABLE(psffit), DP_NFEXTABLE(psffit), xfrom_psf,
+		    yfrom_psf, dvdx, dvdy)
+		if (((rel_bright * deriv[1] + sky) > datamax) && (iter >= 3))
+		    next
+		if (DP_RECENTER(dao) == YES) {
+		    deriv[2] = rel_bright * dvdx
+		    deriv[3] = rel_bright * dvdy
+		}
+		if (DP_FITSKY(dao) == YES)
+		    deriv[nterm] = 1.0
+
+	        # Get the residual from the PSF fit and the pixel
+		# intensity as predicted by the fit
+		d_pixval = (pixval - sky) - rel_bright * deriv[1]
+		pred_pixval = max (0.0, pixval - d_pixval)
+
+		# Calculate the anticipated error in the intensity of
+		# in this pixel including READOUT noise, photon statistics
+		# and the error of interpolating within the PSF.
+
+		sigmasq = pred_pixval / DP_PHOTADU (dao) +  read_noise +
+		    (pererr * pred_pixval) ** 2 + (peakerr *
+		    (pred_pixval - sky)) ** 2
+		if (sigmasq <= 0.0)
+		    next
+		sigma = sqrt (sigmasq)
+		relerr = abs (d_pixval / sigma)
+
+		# Compute the radial wweighting function.
+		weight = 5.0 / (5.0 + radsq / (1.0 - radsq))
+
+		# Now add this pixel into the quantities which go to make
+		# up the sharpness index.
+		rhosq = dxsq / Memr[DP_PSFPARS(psffit)] ** 2 + dysq /
+		        Memr[DP_PSFPARS(psffit)+1] ** 2
+
+		# Include in the sharpness calculation only those pixels
+		# which are within NCORE_SIGMASQ core-sigmasq of the
+		# centroid. This saves time and floating underflows
+		# bu excluding pixels which contribute less than one
+		# part in a million to the fit.
+
+		if (rhosq <= NCORE_SIGMASQ) {
+		    rhosq = 0.6931472 * rhosq
+		    dfdsig = exp (-rhosq) * (rhosq - 1.0)
+		    #pred_pixval = max (0.0, pixval - sky) + sky
+		    numer = numer + dfdsig * d_pixval / sigmasq
+		    denom = denom + (dfdsig ** 2) / sigmasq
+		}
+
+
+		# Derive the weight of this pixel. First of all the weight
+		# depends on the distance of the pixel from the centroid of
+		# the star --- it is determined from a function which is very
+		# nearly unity for radii much smaller than the fitting radius,
+		# and which goes to zero for radii very near the fitting
+		# radius. Then reject any pixels with 10 sigma residuals
+		# after the first iteration. 
+
+		chi = chi + weight * relerr
+		sum_weight = sum_weight + weight
+
+		# Now the weight is scaled to the inverse square of the
+		# expected mean error.
+
+		weight = weight / sigmasq
+
+		# Reduce the weight of a bad pixel. A pixel having a residual
+		# of 2.5 sigma gets reduced to half weight; a pixel having
+		# a residual of of 5.0 sigma gets weight 1/257.
+
+		if (clip) 
+		    weight = weight / (1.0 + (relerr / (DP_CLIPRANGE(dao) *
+		        chiold)) ** DP_CLIPEXP(dao))
+
+		# Now add this pixel into the vector of residuals
+		# and the normal matrix.
+		do i = 1, nterm {
+		    resid[i] = resid[i] + d_pixval * deriv[i] * weight
+		    do j = 1, nterm {
+			normal[i,j] = normal[i,j] + deriv[i] * deriv[j] *
+			    weight
+		    }
+		}
+
+		npix = npix + 1
+	    }
+	}
+
+	# Compute the goodness of fit index CHI. CHI is pulled toward its
+	# expected value of unity before being stored in CHIOLD to keep the
+	# statistics of a small number of pixels from dominating the error
+	# analysis.
+
+	if (sum_weight > nterm) {
+	    chi = CHI_NORM * chi / sqrt (sum_weight * (sum_weight - nterm))
+	    chiold = ((sum_weight - MIN_SUMWT) * chi + MIN_SUMWT) / sum_weight
+	} else {
+	    chi = 1.0
+	    chiold = 1.0
+	}
+
+	return (npix)
+end
diff --git a/noao/digiphot/daophot/peak/dppkwrite.x b/noao/digiphot/daophot/peak/dppkwrite.x
new file mode 100644
index 00000000..97f770bf
--- /dev/null
+++ b/noao/digiphot/daophot/peak/dppkwrite.x
@@ -0,0 +1,365 @@
+include	<tbset.h>
+include <time.h>
+include "../lib/daophotdef.h"
+include "../lib/apseldef.h"
+include "../lib/peakdef.h"
+
+# DP_TNEWPEAK -- Create a new output PEAK ST table.
+
+procedure dp_tnewpeak (dao, tp, colpoint)
+
+pointer	dao		# pointer to the daophot strucuture
+pointer	tp		# pointer to the output table
+int	colpoint[ARB]	# array of column pointers
+
+int	i
+pointer	sp, colnames, colunits, colformat, col_dtype, col_len
+
+begin
+	# Allocate space for the table definition.
+	call smark (sp)
+	call salloc (colnames, PK_NOUTCOL * (SZ_COLNAME + 1), TY_CHAR)
+	call salloc (colunits, PK_NOUTCOL * (SZ_COLUNITS + 1), TY_CHAR)
+	call salloc (colformat, PK_NOUTCOL * (SZ_COLFMT + 1), TY_CHAR)
+	call salloc (col_dtype, PK_NOUTCOL, TY_INT)
+	call salloc (col_len, PK_NOUTCOL, TY_INT)
+
+	# Set up the column definitions.
+	call strcpy (ID, Memc[colnames], SZ_COLNAME)
+	call strcpy (XCENTER, Memc[colnames+SZ_COLNAME+1], SZ_COLNAME)
+	call strcpy (YCENTER, Memc[colnames+2*SZ_COLNAME+2], SZ_COLNAME)
+	call strcpy (MAG, Memc[colnames+3*SZ_COLNAME+3], SZ_COLNAME)
+	call strcpy (MAGERR, Memc[colnames+4*SZ_COLNAME+4], SZ_COLNAME)
+	call strcpy (SKY, Memc[colnames+5*SZ_COLNAME+5], SZ_COLNAME)
+	call strcpy (NITER, Memc[colnames+6*SZ_COLNAME+6], SZ_COLNAME)
+	call strcpy (SHARP, Memc[colnames+7*SZ_COLNAME+7], SZ_COLNAME)
+	call strcpy (CHI, Memc[colnames+8*SZ_COLNAME+8], SZ_COLNAME)
+	call strcpy (PIER, Memc[colnames+9*SZ_COLNAME+9], SZ_COLNAME)
+	call strcpy (PERROR, Memc[colnames+10*SZ_COLNAME+10], SZ_COLNAME)
+
+	# Define the column formats.
+	call strcpy ("%6d", Memc[colformat], SZ_COLFMT)
+	call strcpy ("%10.3f", Memc[colformat+SZ_COLFMT+1], SZ_COLFMT)
+	call strcpy ("%10.3f", Memc[colformat+2*SZ_COLFMT+2], SZ_COLFMT)
+	call strcpy ("%12.3f", Memc[colformat+3*SZ_COLFMT+3], SZ_COLFMT)
+	call strcpy ("%14.3f", Memc[colformat+4*SZ_COLFMT+4], SZ_COLFMT)
+	call strcpy ("%15.7g", Memc[colformat+5*SZ_COLFMT+5], SZ_COLFMT)
+	call strcpy ("%12.3f", Memc[colformat+6*SZ_COLFMT+6], SZ_COLFMT)
+	call strcpy ("%12.3f", Memc[colformat+7*SZ_COLFMT+7], SZ_COLFMT)
+	call strcpy ("%12.3f", Memc[colformat+8*SZ_COLFMT+8], SZ_COLFMT)
+	call strcpy ("%6d", Memc[colformat+9*SZ_COLFMT+9], SZ_COLFMT)
+	call strcpy ("%13s", Memc[colformat+10*SZ_COLFMT+10], SZ_COLFMT)
+
+	# Define the column units.
+	call strcpy ("NUMBER", Memc[colunits], SZ_COLUNITS)
+	call strcpy ("PIXELS", Memc[colunits+SZ_COLUNITS+1], SZ_COLUNITS)
+	call strcpy ("PIXELS", Memc[colunits+2*SZ_COLUNITS+2], SZ_COLUNITS)
+	call strcpy ("MAGNITUDES", Memc[colunits+3*SZ_COLUNITS+3], SZ_COLUNITS)
+	call strcpy ("MAGNITUDES", Memc[colunits+4*SZ_COLUNITS+4], SZ_COLUNITS)
+	call strcpy ("ADU", Memc[colunits+5*SZ_COLUNITS+5], SZ_COLUNITS)
+	call strcpy ("NUMBER", Memc[colunits+6*SZ_COLUNITS+6], SZ_COLUNITS)
+	call strcpy ("NUMBER", Memc[colunits+7*SZ_COLUNITS+7], SZ_COLUNITS)
+	call strcpy ("NUMBER", Memc[colunits+8*SZ_COLUNITS+8], SZ_COLUNITS)
+	call strcpy ("NUMBER", Memc[colunits+9*SZ_COLUNITS+9], SZ_COLUNITS)
+	call strcpy ("PERRORS", Memc[colunits+10*SZ_COLUNITS+10], SZ_COLUNITS)
+
+	# Define the column data types.
+	Memi[col_dtype] = TY_INT
+	Memi[col_dtype+1] = TY_REAL
+	Memi[col_dtype+2] = TY_REAL
+	Memi[col_dtype+3] = TY_REAL
+	Memi[col_dtype+4] = TY_REAL
+	Memi[col_dtype+5] = TY_REAL
+	Memi[col_dtype+6] = TY_INT
+	Memi[col_dtype+7] = TY_REAL
+	Memi[col_dtype+8] = TY_REAL
+	Memi[col_dtype+9] = TY_INT
+	Memi[col_dtype+10] = -13
+
+	# Define the column lengths.
+	do i = 1, PK_NOUTCOL 
+	    Memi[col_len+i-1] = 1
+	
+	# Define and create the table.
+	call tbcdef (tp, colpoint, Memc[colnames], Memc[colunits],
+	    Memc[colformat], Memi[col_dtype], Memi[col_len], PK_NOUTCOL)
+	call tbtcre (tp)
+
+	# Write out some header parameters.
+	call dp_tpeakpars (dao, tp)
+
+	call sfree (sp)
+end
+
+
+define PK_NAME1STR "#N%4tID%10tXCENTER%20tYCENTER%30tMAG%42tMERR%56tMSKY%71t\
+NITER%80t\\\n"
+define PK_UNIT1STR "#U%4t##%10tpixels%20tpixels%30tmagnitudes%42t\
+magnitudes%56tcounts%71t##%80t\\\n"
+define PK_FORMAT1STR "#F%4t%%-9d%10t%%-10.3f%20t%%-10.3f%30t%%-12.3f%42t\
+%%-14.3f%56t%%-15.7g%71t%%-6d%80t \n"
+
+define PK_NAME2STR "#N%12tSHARPNESS%24tCHI%36tPIER%42tPERROR%80t\\\n"
+define PK_UNIT2STR "#U%12t##%24t##%36t##%42tperrors%80t\\\n"
+define PK_FORMAT2STR "#F%12t%%-23.3f%24t%%-12.3f%36t%%-6d%42t%%-13s%80t \n"
+
+
+# DP_XNEWPEAK -- Create a new PEAK output text file.
+
+procedure dp_xnewpeak (dao, tp)
+
+pointer	dao		# pointer to the daophot structure
+int	tp		# the output file descriptor
+
+begin
+	# Write out some header parameters.
+	call dp_xpeakpars (dao, tp)
+
+	# Print out the banner file.
+	call fprintf (tp, "#\n")
+	call fprintf (tp, PK_NAME1STR)
+	call fprintf (tp, PK_UNIT1STR)
+	call fprintf (tp, PK_FORMAT1STR)
+	call fprintf (tp, "#\n")
+	call fprintf (tp, PK_NAME2STR)
+	call fprintf (tp, PK_UNIT2STR)
+	call fprintf (tp, PK_FORMAT2STR)
+	call fprintf (tp, "#\n")
+end
+
+
+# DP_TPEAKPARS -- Add parameters to the header of the PEAK ST table.
+
+procedure dp_tpeakpars (dao, tp)
+
+pointer	dao			# pointer to the daophot structure
+pointer	tp			# pointer to the output table
+
+pointer	sp, psffit, outstr, date, time
+bool	itob()
+int	envfind()
+
+begin
+	# Define some daophot pointers.
+	psffit = DP_PSFFIT(dao)
+
+	# Allocate working space.
+	call smark (sp)
+	call salloc (outstr, SZ_LINE, TY_CHAR)
+	call salloc (date, SZ_DATE, TY_CHAR)
+	call salloc (time, SZ_DATE, TY_CHAR)
+
+	# Write the id.
+	if (envfind ("version", Memc[outstr], SZ_LINE) <= 0)
+	    call strcpy ("NOAO/IRAF", Memc[outstr], SZ_LINE)
+	call dp_rmwhite (Memc[outstr], Memc[outstr], SZ_LINE)
+	call tbhadt (tp, "IRAF", Memc[outstr])
+	if (envfind ("userid", Memc[outstr], SZ_LINE) <= 0)
+	    call tbhadt (tp, "USER", Memc[outstr])
+	call gethost (Memc[outstr], SZ_LINE)
+	call tbhadt (tp, "HOST", Memc[outstr])
+	call dp_date (Memc[date], Memc[time], SZ_DATE)
+	call tbhadt (tp, "DATE", Memc[date])
+	call tbhadt (tp, "TIME", Memc[time])
+	call tbhadt (tp, "PACKAGE", "daophot")
+	call tbhadt (tp, "TASK", "peak")
+
+	# Write out the files names.
+	call dp_imroot (DP_INIMAGE(dao), Memc[outstr], SZ_LINE)
+	call tbhadt (tp, "IMAGE", Memc[outstr])
+	call dp_froot (DP_INPHOTFILE(dao), Memc[outstr], SZ_LINE)
+	call tbhadt (tp, "PHOTFILE", Memc[outstr])
+	call dp_imroot (DP_PSFIMAGE(dao), Memc[outstr], SZ_LINE)
+	call tbhadt (tp, "PSFIMAGE", Memc[outstr])
+	call dp_froot (DP_OUTPHOTFILE(dao), Memc[outstr], SZ_LINE)
+	call tbhadt (tp, "PEAKFILE", Memc[outstr])
+	if (DP_OUTREJFILE(dao) == EOS)
+	    call tbhadt (tp, "REJFILE", "\"\"")
+	else {
+	    call dp_froot (DP_OUTREJFILE(dao), Memc[outstr], SZ_LINE)
+	    call tbhadt (tp, "REJFILE", Memc[outstr])
+	}
+
+	# Define the data characteristics.
+	call tbhadr (tp, "SCALE", DP_SCALE(dao))
+	call tbhadr (tp, "DATAMIN", DP_MINGDATA(dao))
+	call tbhadr (tp, "DATAMAX", DP_MAXGDATA(dao))
+	call tbhadr (tp, "GAIN", DP_PHOTADU(dao))
+	call tbhadr (tp, "READNOISE", DP_READNOISE(dao))
+
+	# Define the observing parameters.
+	call tbhadt (tp, "OTIME", DP_OTIME(dao))
+	call tbhadr (tp, "XAIRMASS", DP_XAIRMASS(dao))
+	call tbhadt (tp, "IFILTER", DP_IFILTER(dao))
+
+	# Define the fitting parameters.
+	call tbhadb (tp, "RECENTER", itob (DP_RECENTER(dao)))
+	call tbhadb (tp, "FITSKY", itob (DP_FITSKY(dao)))
+	call tbhadr (tp, "PSFRAD", DP_SPSFRAD(dao))
+	call tbhadr (tp, "FITRAD", DP_SFITRAD(dao))
+	call tbhadr (tp, "PSFMAG", DP_PSFMAG(psffit))
+	call tbhadi (tp, "MAXITER", DP_MAXITER(dao))
+	call tbhadr (tp, "FLATERROR", DP_FLATERR(dao))
+	call tbhadr (tp, "PROFERROR", DP_PROFERR(dao))
+	call tbhadi (tp, "CLIPEXP", DP_CLIPEXP(dao))
+	call tbhadr (tp, "CLIPRANGE", DP_CLIPRANGE(dao))
+	#call tbhadi (tp, "MAXNSTAR", DP_MAXNSTAR(dao))
+
+	call sfree(sp)
+end
+
+
+# DP_XPEAKPARS -- Add parameters to the header of the output PEAK text file.
+
+procedure dp_xpeakpars (dao, tp)
+
+pointer	dao			# pointer to the daophot structure
+int	tp			# the output file descriptor
+
+pointer	sp, psffit, outstr, date, time
+bool	itob()
+int	envfind()
+
+begin
+	# Define some daophot pointers.
+	psffit = DP_PSFFIT(dao)
+
+	# Allocate working space.
+	call smark (sp)
+	call salloc (outstr, SZ_LINE, TY_CHAR)
+	call salloc (date, SZ_DATE, TY_CHAR)
+	call salloc (time, SZ_DATE, TY_CHAR)
+
+	# Write the id.
+	if (envfind ("version", Memc[outstr], SZ_LINE) <= 0)
+	    call strcpy ("NOAO/IRAF", Memc[outstr], SZ_LINE)
+	call dp_rmwhite (Memc[outstr], Memc[outstr], SZ_LINE)
+	call dp_sparam (tp, "IRAF", Memc[outstr], "version", "")
+	if (envfind ("userid", Memc[outstr], SZ_LINE) <= 0)
+	    call dp_sparam (tp, "USER", Memc[outstr], "name", "")
+	call gethost (Memc[outstr], SZ_LINE)
+	call dp_sparam (tp, "HOST", Memc[outstr], "computer", "")
+	call dp_date (Memc[date], Memc[time], SZ_DATE)
+	call dp_sparam (tp, "DATE", Memc[date], "yyyy-mm-dd", "")
+	call dp_sparam (tp, "TIME", Memc[time], "hh:mm:ss", "")
+	call dp_sparam (tp, "PACKAGE", "daophot", "name", "")
+	call dp_sparam (tp, "TASK", "peak", "name", "")
+
+	# Write out the files names.
+	call dp_imroot (DP_INIMAGE(dao), Memc[outstr], SZ_LINE)
+	call dp_sparam (tp, "IMAGE", Memc[outstr], "imagename", "")
+	call dp_froot (DP_INPHOTFILE(dao), Memc[outstr], SZ_LINE)
+	call dp_sparam (tp, "PHOTFILE", Memc[outstr], "filename", "")
+	call dp_imroot (DP_PSFIMAGE(dao), Memc[outstr], SZ_LINE)
+	call dp_sparam (tp, "PSFIMAGE", Memc[outstr], "imagename", "")
+	call dp_froot (DP_OUTPHOTFILE(dao), Memc[outstr], SZ_LINE)
+	call dp_sparam (tp, "PEAKFILE", Memc[outstr], "filename", "")
+
+	if (DP_OUTREJFILE(dao) == EOS)
+	    call dp_sparam (tp, "REJFILE", "\"\"", "filename", "")
+	else {
+	    call dp_froot (DP_OUTREJFILE(dao), Memc[outstr], SZ_LINE)
+	    call dp_sparam (tp, "REJFILE", Memc[outstr], "filename", "")
+	}
+
+	# Define the data characteristics.
+	call dp_rparam (tp, "SCALE", DP_SCALE(dao), "units/pix", "")
+	call dp_rparam (tp, "DATAMIN", DP_MINGDATA(dao), "counts", "")
+	call dp_rparam (tp, "DATAMAX", DP_MAXGDATA(dao), "counts", "")
+	call dp_rparam (tp, "GAIN", DP_PHOTADU(dao), "e-/adu", "")
+	call dp_rparam (tp, "READNOISE", DP_READNOISE(dao), "electrons", "")
+
+	# Define the observing parameters.
+	call dp_sparam (tp, "OTIME", DP_OTIME(dao), "timeunit", "")
+	call dp_rparam (tp, "XAIRMASS", DP_XAIRMASS(dao), "number", "")
+	call dp_sparam (tp, "IFILTER", DP_IFILTER(dao), "filter", "")
+
+	# Define the fitting parameters.
+	call dp_bparam (tp, "RECENTER", itob (DP_RECENTER(dao)), "switch", "")
+	call dp_bparam (tp, "FITSKY", itob (DP_FITSKY(dao)), "switch", "")
+	call dp_rparam (tp, "PSFRAD", DP_SPSFRAD(dao), "scaleunit", "")
+	call dp_rparam (tp, "FITRAD", DP_SFITRAD(dao), "scaleunit", "")
+	call dp_rparam (tp, "PSFMAG", DP_PSFMAG(psffit), "magnitudes", "")
+	call dp_iparam (tp, "MAXITER", DP_MAXITER(dao), "number", "")
+	call dp_rparam (tp, "FLATERROR", DP_FLATERR(dao), "percentage", "")
+	call dp_rparam (tp, "PROFERROR", DP_PROFERR(dao), "percentage", "")
+	call dp_iparam (tp, "CLIPEXP", DP_CLIPEXP(dao), "number", "")
+	call dp_rparam (tp, "CLIPRANGE", DP_CLIPRANGE(dao), "sigma", "")
+	#call dp_iparam (tp, "MAXNSTAR", DP_MAXNSTAR(dao), "number", "")
+
+	call sfree(sp)
+end
+
+
+# DP_TPKWRITE -- Write the output photometry record to a PEAK ST table.
+
+procedure dp_tpkwrite (tpout, colpoint, id, x, y, mag, errmag, sky, niter, chi,
+	sharp, pier, perror, plen, star)
+
+pointer	tpout		# pointer to the output table
+int	colpoint[ARB]	# array of column pointers
+int	id		# id of the star
+real	x, y		# position of the star
+real	mag		# magnitude of the star
+real	errmag		# error magnitude of the star
+real	sky		# value of sky
+int	niter		# number of iterations
+real	chi		# chi squared value
+real	sharp		# sharpness characteristic
+int	pier		# the error code
+char	perror[ARB]	# the error string
+int	plen		# the length of the error code string
+int	star		# row number
+
+begin
+	call tbrpti (tpout, colpoint[1], id, 1, star)
+	call tbrptr (tpout, colpoint[2], x, 1, star)
+	call tbrptr (tpout, colpoint[3], y, 1, star)
+	call tbrptr (tpout, colpoint[4], mag, 1, star)
+	call tbrptr (tpout, colpoint[5], errmag, 1, star)
+	call tbrptr (tpout, colpoint[6], sky, 1, star)
+	call tbrpti (tpout, colpoint[7], niter, 1, star)
+	call tbrptr (tpout, colpoint[8], sharp, 1, star)
+	call tbrptr (tpout, colpoint[9], chi, 1, star)
+	call tbrpti (tpout, colpoint[10], pier, 1, star)
+	call tbrptt (tpout, colpoint[11], perror, plen, 1, star)
+end
+
+
+define PK_DATA1STR "%-9d%10t%-10.3f%20t%-10.3f%30t%-12.3f%42t%-14.3f%56t\
+%-15.7g%71t%-6d%80t\\\n"
+define PK_DATA2STR "%12t%-12.3f%24t%-12.3f%36t%-6d%42t%-13.13s%80t \n"
+
+# DP_XPKWRITE -- Write the output photometry record to a PEAK text file.
+
+procedure dp_xpkwrite (tpout, id, x, y, mag, errmag, sky, niter, chi, sharp,
+	pier, perror, plen)
+
+int 	tpout		# the output file descriptor
+int	id		# id of the star
+real	x, y		# position of the star
+real	mag		# magnitude of the star
+real	errmag		# error magnitude of the star
+real	sky		# value of sky
+int	niter		# number of iterations
+real	chi		# chi squared value
+real	sharp		# sharpness characteristic
+int	pier		# the error code
+char	perror[ARB]	# the error string
+int	plen		# the length of the error code string
+
+begin
+	call fprintf (tpout, PK_DATA1STR)
+	    call pargi (id)
+	    call pargr (x)
+	    call pargr (y)
+	    call pargr (mag)
+	    call pargr (errmag)
+	    call pargr (sky)
+	    call pargi (niter)
+	call fprintf (tpout, PK_DATA2STR)
+	    call pargr (sharp)
+	    call pargr (chi)
+	    call pargi (pier)
+	    call pargstr (perror)
+end
diff --git a/noao/digiphot/daophot/peak/dprrphot.x b/noao/digiphot/daophot/peak/dprrphot.x
new file mode 100644
index 00000000..36aed88f
--- /dev/null
+++ b/noao/digiphot/daophot/peak/dprrphot.x
@@ -0,0 +1,98 @@
+include	"../lib/apseldef.h"
+
+# DP_TPKINIT -- Procedure to initialize for reading the "standard" fields from
+# a photometry table. The "standard" fields being ID, X, Y, MAG, ERR, and SKY
+
+procedure dp_tpkinit (tp, colpoint)
+
+pointer	tp			# the table descriptor
+pointer	colpoint[ARB]		# the column descriptor
+
+begin
+	# Get the results one by one
+	# First the ID
+	call tbcfnd (tp, ID, colpoint[1], 1)
+	if (colpoint[1] == NULL) {
+	    call eprintf ("Column %s not found\n")
+	        call pargstr (ID)
+	}
+
+	# Then the position
+	call tbcfnd (tp, XCENTER, colpoint[2], 1)
+	if (colpoint[2] == NULL) {
+	    call eprintf ("Column %s not found\n")
+	        call pargstr (XCENTER)
+	}
+
+	call tbcfnd (tp, YCENTER, colpoint[3], 1)
+	if (colpoint[3] == NULL) {
+	    call eprintf ("Column %s not found\n")
+	        call pargstr (YCENTER)
+	}
+
+	# Now the Magnitude
+	call tbcfnd (tp, MAG, colpoint[4], 1)
+	if (colpoint[4] == NULL)		# No column
+	    call tbcfnd (tp, APMAG, colpoint[4], 1)
+	if (colpoint[4] == NULL) {
+	    call eprintf ("Column %s not found\n")
+	    call pargstr (APMAG)
+	}
+
+	# The sky
+	call tbcfnd (tp, SKY, colpoint[5], 1)
+	if (colpoint[5] == NULL)
+	    call tbcfnd (tp, SKY, colpoint[5], 1)
+	if (colpoint[5] == NULL) {
+	    call eprintf ("Column %s not found\n")
+	        call pargstr (SKY)
+	}
+end
+
+
+# DP_RRPHOT -- Fetch the photometry for a single star from either a
+# table or a text photometry file.
+
+int procedure dp_rrphot (tp, key, fields, indices, id, x, y, sky, mag,
+	instar, nrow)
+
+int	tp		# the input file descriptor
+pointer	key		# pointer to text apphot structure
+char	fields[ARB]	# character fields
+int	indices[ARB]	# columns pointers
+int	id		# star id
+real	x		# x center value
+real	y		# y center value
+real	sky		# sky value
+real	mag		# magnitude
+int	instar		# current record
+int	nrow		# maximum number of rows for ST table
+
+bool	nullflag
+int	nrec
+int	dp_apsel()
+
+begin
+	# If nrow is 0 the file file is a text file otherwise it is a table.
+
+	if (nrow == 0) {
+
+	    nrec = dp_apsel (key, tp, fields, indices, id, x, y, sky, mag)
+	    if (nrec != EOF)
+		instar = instar + 1
+
+	} else if ((instar + 1) <= nrow) {
+
+	    instar = instar + 1
+	    call tbrgti (tp, indices[1], id, nullflag, 1, instar)
+	    call tbrgtr (tp, indices[2], x, nullflag, 1, instar)
+	    call tbrgtr (tp, indices[3], y, nullflag, 1, instar)
+	    call tbrgtr (tp, indices[4], mag, nullflag, 1, instar)
+	    call tbrgtr (tp, indices[5], sky, nullflag, 1, instar)
+	    nrec = instar
+
+	} else
+	    nrec = EOF
+
+	return (nrec)
+end
diff --git a/noao/digiphot/daophot/peak/mkpkg b/noao/digiphot/daophot/peak/mkpkg
new file mode 100644
index 00000000..f3570728
--- /dev/null
+++ b/noao/digiphot/daophot/peak/mkpkg
@@ -0,0 +1,22 @@
+# PEAK task
+
+$checkout libpkg.a ".."
+$update libpkg.a
+$checkin libpkg.a ".."
+$exit
+
+libpkg.a:
+	dpmempk.x	../lib/daophotdef.h	../lib/peakdef.h
+	dppeakphot.x	<imhdr.h>		<tbset.h>                   \
+			../lib/daophotdef.h	../lib/apseldef.h           \
+			../lib/peakdef.h
+	dppkconfirm.x
+	dppkfit.x	<mach.h>                ../lib/daophotdef.h         \
+			../lib/peakdef.h
+	dppkwrite.x	<tbset.h>		<time.h>		    \
+			../lib/daophotdef.h     ../lib/apseldef.h           \
+			../lib/peakdef.h
+	dprrphot.x      ../lib/apseldef.h
+	t_peak.x	<fset.h>		<imhdr.h>                   \
+	                ../lib/daophotdef.h
+	;
diff --git a/noao/digiphot/daophot/peak/t_peak.x b/noao/digiphot/daophot/peak/t_peak.x
new file mode 100644
index 00000000..2022fe16
--- /dev/null
+++ b/noao/digiphot/daophot/peak/t_peak.x
@@ -0,0 +1,299 @@
+include	<fset.h>
+include <imhdr.h>
+include	"../lib/daophotdef.h"
+
+# T_PEAK  -- Procedure to fit the PSF to single stars.
+
+procedure t_peak ()
+
+pointer	image				# name of the image
+pointer	photfile			# input photometry file
+pointer	psfimage			# the PSF image
+pointer	peakfile			# output PEAK photometry file
+pointer	rejfile				# output PEAK rejections file
+
+pointer	sp, im, psfim, outfname, dao, str
+int	apd, root, verbose, verify, cache, update, pkfd, rejfd
+int	imlist, limlist, alist, lalist, pimlist, lpimlist, olist, lolist
+int	rlist, lrlist, wcs, req_size, old_size, buf_size, memstat
+bool	ap_text
+
+pointer	immap(), tbtopn()
+int	open(), fnldir(), strlen(), strncmp(), fstati(), btoi()
+int	access(), imtopen(), imtlen(), imtgetim(), fntopnb(), fntlenb()
+int	fntgfnb(), clgwrd(), sizeof(), dp_memstat()
+bool	clgetb(), itob()
+
+begin
+	# Set the standard output to flush on newline.
+	if (fstati (STDOUT, F_REDIR) == NO)
+	    call fseti (STDOUT, F_FLUSHNL, YES)
+
+	# Get some working memory.
+	call smark (sp)
+	call salloc (image, SZ_FNAME, TY_CHAR)
+	call salloc (photfile, SZ_FNAME, TY_CHAR)
+	call salloc (psfimage, SZ_FNAME, TY_CHAR)
+	call salloc (peakfile, SZ_FNAME, TY_CHAR)
+	call salloc (rejfile, SZ_FNAME, TY_CHAR)
+	call salloc (outfname, SZ_FNAME, TY_CHAR)
+	call salloc (str, SZ_FNAME, TY_CHAR)
+
+	# Get the various task parameters
+	call clgstr ("image", Memc[image], SZ_FNAME)
+	call clgstr ("photfile", Memc[photfile], SZ_FNAME)
+	call clgstr ("psfimage", Memc[psfimage], SZ_FNAME)
+	call clgstr ("peakfile", Memc[peakfile], SZ_FNAME)
+	call clgstr ("rejfile", Memc[rejfile], SZ_FNAME)
+	verbose = btoi (clgetb ("verbose"))
+	verify = btoi (clgetb ("verify"))
+	update = btoi (clgetb ("update"))
+	cache = btoi (clgetb ("cache"))
+
+	# Get the lists.
+	imlist = imtopen (Memc[image])
+	limlist = imtlen (imlist)
+	alist = fntopnb (Memc[photfile], NO)
+	lalist = fntlenb (alist)
+	pimlist = imtopen (Memc[psfimage])
+	lpimlist = imtlen (pimlist)
+	olist = fntopnb (Memc[peakfile], NO)
+	lolist =  fntlenb (olist)
+	rlist = fntopnb (Memc[rejfile], NO)
+	lrlist =  fntlenb (rlist)
+
+	# Test that the lengths of the photometry file, psf image, and
+	# output file lists are the same as the length of the input image
+	# list.
+
+	if ((limlist != lalist) && (strncmp (Memc[photfile],
+	    DEF_DEFNAME, DEF_LENDEFNAME) != 0)) {
+	    call imtclose (imlist)
+	    call fntclsb (alist)
+	    call imtclose (pimlist)
+	    call fntclsb (olist)
+	    call fntclsb (rlist)
+	    call sfree (sp)
+	    call error (0,
+	        "Incompatable image and input photometry file list lengths")
+	}
+
+	if ((limlist != lpimlist) && (strncmp (Memc[psfimage],
+	    DEF_DEFNAME, DEF_LENDEFNAME) != 0)) {
+	    call imtclose (imlist)
+	    call fntclsb (alist)
+	    call imtclose (pimlist)
+	    call fntclsb (olist)
+	    call fntclsb (rlist)
+	    call sfree (sp)
+	    call error (0,
+	        "Incompatable image and psf image list lengths")
+	}
+
+	if ((limlist != lolist) && (strncmp (Memc[peakfile],
+	    DEF_DEFNAME, DEF_LENDEFNAME) != 0)) {
+	    call imtclose (imlist)
+	    call fntclsb (alist)
+	    call imtclose (pimlist)
+	    call fntclsb (olist)
+	    call fntclsb (rlist)
+	    call sfree (sp)
+	    call error (0,
+	        "Incompatable image and output photometry file list lengths")
+	}
+
+	if ((lrlist != 0) && (limlist != lolist) && (strncmp (Memc[rejfile],
+	    DEF_DEFNAME, DEF_LENDEFNAME) != 0)) {
+	    call imtclose (imlist)
+	    call fntclsb (alist)
+	    call imtclose (pimlist)
+	    call fntclsb (olist)
+	    call fntclsb (rlist)
+	    call sfree (sp)
+	    call error (0,
+	        "Incompatable image and rejections file list lengths")
+	}
+
+	# Initialize the DAOPHOT structure, and get the pset parameters.
+	call dp_gppars (dao)	
+	call dp_seti (dao, VERBOSE, verbose)
+
+	# Verify the standard algorithm parameters.
+	if (verify == YES) {
+	    call dp_pkconfirm (dao)
+	    if (update == YES)
+		call dp_pppars (dao)
+	}
+
+        # Get the wcs information.
+        wcs = clgwrd ("wcsin", Memc[str], SZ_FNAME, WCSINSTR)
+        if (wcs <= 0) {
+            call eprintf (
+                "Warning: Setting the input coordinate system to logical\n")
+            wcs = WCS_LOGICAL
+        }
+        call dp_seti (dao, WCSIN, wcs)
+        wcs = clgwrd ("wcsout", Memc[str], SZ_FNAME, WCSOUTSTR)
+        if (wcs <= 0) {
+            call eprintf (
+                "Warning: Setting the output coordinate system to logical\n")
+            wcs = WCS_LOGICAL
+        }
+        call dp_seti (dao, WCSOUT, wcs)
+        wcs = clgwrd ("wcspsf", Memc[str], SZ_FNAME, WCSPSFSTR)
+        if (wcs <= 0) {
+            call eprintf (
+                "Warning: Setting the psf coordinate system to logical\n")
+            wcs = WCS_LOGICAL
+        }
+        call dp_seti (dao, WCSPSF, wcs)
+
+	# Initialize the PSF structure.
+	call dp_fitsetup (dao)
+
+	# Initialize the PEAK structure.
+	call dp_pksetup (dao)
+
+	# Loop over the list of images.
+	while (imtgetim (imlist, Memc[image], SZ_FNAME) != EOF) {
+
+	    # Open the image.
+	    im = immap (Memc[image], READ_ONLY, 0)		
+	    call dp_imkeys (dao, im)
+	    call dp_sets (dao, INIMAGE, Memc[image])
+
+            # Cache the input image pixels.
+            req_size = MEMFUDGE * IM_LEN(im,1) * IM_LEN(im,2) *
+                sizeof (IM_PIXTYPE(im))
+            memstat = dp_memstat (cache, req_size, old_size)
+            if (memstat == YES)
+                call dp_pcache (im, INDEFI, buf_size)
+
+	    # Open the input photometry file.
+	    if (fntgfnb (alist, Memc[photfile], SZ_FNAME) == EOF)
+		call strcpy (DEF_DEFNAME, Memc[photfile], SZ_FNAME)
+	    root = fnldir (Memc[photfile], Memc[outfname], SZ_FNAME)
+	    if (strncmp (DEF_DEFNAME, Memc[photfile+root],
+	        DEF_LENDEFNAME) == 0 || root == strlen (Memc[photfile]))
+	        call dp_inname (Memc[image], Memc[outfname], "mag",
+		    Memc[outfname], SZ_FNAME)
+	    else
+	        call strcpy (Memc[photfile], Memc[outfname], SZ_FNAME)
+	    ap_text =  itob (access (Memc[outfname], 0, TEXT_FILE))
+	    if (ap_text)
+	        apd = open (Memc[outfname], READ_ONLY, TEXT_FILE)
+	    else
+	        apd = tbtopn (Memc[outfname], READ_ONLY, 0)
+	    call dp_sets (dao, INPHOTFILE, Memc[outfname])
+
+	    # Read in the PSF function.
+	    if (imtgetim (pimlist, Memc[psfimage], SZ_FNAME) == EOF)
+		call strcpy (DEF_DEFNAME, Memc[psfimage], SZ_FNAME)
+	    root = fnldir (Memc[psfimage], Memc[outfname], SZ_FNAME)
+	    if (strncmp (DEF_DEFNAME, Memc[psfimage+root],
+	        DEF_LENDEFNAME) == 0 || root == strlen (Memc[psfimage]))
+	        call dp_iimname (Memc[image], Memc[outfname], "psf",
+		    Memc[outfname], SZ_FNAME)
+	    else
+	        call strcpy (Memc[psfimage], Memc[outfname], SZ_FNAME)
+	    psfim = immap (Memc[outfname], READ_ONLY, 0)
+	    call dp_readpsf (dao, psfim)
+	    call dp_sets (dao, PSFIMAGE, Memc[outfname])
+	
+	    # Open the output photometry file. If the output is "default",
+	    # dir$default or a directory specification then the extension .pk
+	    # is added to the image name and a suitable version number is
+	    # appended to the output name.
+
+	    if (fntgfnb (olist, Memc[peakfile], SZ_FNAME) == EOF)
+		call strcpy (DEF_DEFNAME, Memc[peakfile], SZ_FNAME)
+	    root = fnldir (Memc[peakfile], Memc[outfname], SZ_FNAME)
+	    if (strncmp (DEF_DEFNAME, Memc[peakfile + root],
+	        DEF_LENDEFNAME) == 0 || root == strlen (Memc[peakfile]))
+	        call dp_outname (Memc[image], Memc[outfname], "pk",
+		    Memc[outfname], SZ_FNAME)
+	    else
+	        call strcpy (Memc[peakfile], Memc[outfname], SZ_FNAME)
+	    if (DP_TEXT(dao) == YES)
+	        pkfd = open (Memc[outfname], NEW_FILE, TEXT_FILE)
+	    else
+	        pkfd = tbtopn (Memc[outfname], NEW_FILE, 0)
+	    call dp_sets (dao, OUTPHOTFILE, Memc[outfname])
+
+	    # Open the output rejections file if any are defined. If the
+	    # output is "default", dir$default or a directory specification
+	    # then the extension .prj is added to the image name and a
+	    # suitable version number is appended to the output file name.
+
+	    if (lrlist <= 0) {
+		rejfd = NULL
+		Memc[outfname] = EOS
+	    } else {
+	        if (fntgfnb (rlist, Memc[rejfile], SZ_FNAME) == EOF)
+		    call strcpy (DEF_DEFNAME, Memc[rejfile], SZ_FNAME)
+	        root = fnldir (Memc[rejfile], Memc[outfname], SZ_FNAME)
+	        if (strncmp (DEF_DEFNAME, Memc[rejfile+root],
+		    DEF_LENDEFNAME) == 0 || root == strlen (Memc[rejfile]))
+	            call dp_outname (Memc[image], Memc[outfname], "prj",
+		        Memc[outfname], SZ_FNAME)
+	        else
+	            call strcpy (Memc[rejfile], Memc[outfname], SZ_FNAME)
+	        if (DP_TEXT(dao) == YES)
+	            rejfd = open (Memc[outfname], NEW_FILE, TEXT_FILE)
+	        else
+	            rejfd = tbtopn (Memc[outfname], NEW_FILE, 0)
+	    }
+	    call dp_sets (dao, OUTREJFILE, Memc[outfname])
+
+	    # Now go and do the PSF fitting.
+	    call dp_peakphot (dao, im, apd, pkfd, rejfd, ap_text)
+
+	    # Close the input image.
+	    call imunmap (im)
+
+	    # Close the input photometry file. 
+	    if (ap_text)
+	        call close (apd)
+	    else
+	        call tbtclo (apd)
+
+	    # Close the PSF image.
+	    call imunmap (psfim)
+
+	    # Close the output photometry file.
+	    if (DP_TEXT(dao) == YES)
+	        call close (pkfd)
+	    else
+	        call tbtclo (pkfd)
+
+	    # Close the output rejections file.
+	    if (rejfd != NULL) {
+	        if (DP_TEXT(dao) == YES)
+	            call close (rejfd)
+	        else
+	            call tbtclo (rejfd)
+	    }
+
+            # Uncache memory.
+            call fixmem (old_size)
+
+	}
+
+	# Close the image/file lists.
+	call imtclose (imlist)
+	call fntclsb (alist)
+	call imtclose (pimlist)
+	call fntclsb (olist)
+	call fntclsb (rlist)
+
+	# Free the PEAK structure.
+	call dp_pkclose (dao)
+
+	# Free the PSF structure.
+	call dp_fitclose (dao)
+	
+	# Free the daophot structure.
+	call dp_free (dao)
+
+	call sfree(sp)
+end