Initial commit

13 files changed, 4641 insertions, 0 deletions

diff --git a/noao/digiphot/daophot/allstar/dpabuf.x b/noao/digiphot/daophot/allstar/dpabuf.x
new file mode 100644
index 00000000..13aac109
--- /dev/null
+++ b/noao/digiphot/daophot/allstar/dpabuf.x
@@ -0,0 +1,495 @@
+include	<imhdr.h>
+include "../lib/daophotdef.h"
+include "../lib/allstardef.h"
+
+# DP_GWT -- Return a pointer to the desired weights pixels. The weights buffer
+# pointer is assumed to be NULL on the initial all. If the new line limits
+# are totally contained within the old line limits no action is taken and
+# the old line limits are returned. If the weights pixels are stored in a 
+# scratch image the pixels must be accessed  sequentially, otherwise an
+# error is returned.
+
+pointer procedure dp_gwt (dao, im, line1, line2, mode, flush)
+
+pointer	dao			# pointer to the daophot strucuture
+pointer	im			# pointer to the input image
+int	line1			# the lower line limit
+int	line2			# the upper line limit
+int	mode			# input / output mode
+int	flush 			# flush the output
+
+int	nx, ny, xoff, yoff, llast1, llast2
+pointer	allstar
+
+begin
+	allstar = DP_ALLSTAR(dao)
+	if (DP_WBUF(allstar) == NULL) {
+	    llast1 = 0
+	    llast2 = 0
+	} else if (flush == YES) {
+	    line1 = llast1
+	    line2 = llast2
+	} else if  ((line1 >= llast1) && (line2 <= llast2)) {
+	    line1 = llast1
+	    line2 = llast2
+	    return (DP_WBUF(allstar))
+	} else if (line1 < llast1)
+	    call error (0,
+	        "ERROR: Attempting to access the weights pixels randomly")
+
+	# All the data is cached.
+	if (DP_CACHE (allstar, A_WEIGHT) == YES) {
+
+	    nx = IM_LEN(im,1)
+	    ny = IM_LEN(im,2)
+	    xoff = 1
+	    yoff = 1
+	    DP_WBUF(allstar) = DP_WEIGHTS(allstar)
+
+	# Read in some new data.
+	} else if (mode == READ_ONLY) {
+
+            call dp_albufl2r (DP_WEIGHTS(allstar), DP_WBUF(allstar),
+	        llast1, llast2, line1, line2, flush)
+
+	    if (flush == NO) {
+		nx = IM_LEN(im,1)
+	        ny = line2 - line1 + 1
+		xoff = 1
+	        yoff = line1
+	    } else {
+		nx = 0
+		ny = 0
+		xoff = 0
+		yoff = 0
+	    }
+
+	# Write out the old data and read in some new data.
+	} else if (mode == READ_WRITE) {
+
+            call dp_albufl2rw (DP_WEIGHTS(allstar), DP_WEIGHTS(allstar),
+	        DP_WBUF(allstar), llast1, llast2, line1, line2, flush)
+
+	    if (flush == NO) {
+		nx = IM_LEN(im,1)
+	        ny = line2 - line1 + 1
+		xoff = 1
+	        yoff = line1
+	    } else {
+		nx = 0
+		ny = 0
+		xoff = 0
+		yoff = 0
+	    }
+	}
+
+	# Update the required buffer parameters.
+	DP_WNX(allstar) = nx
+	DP_WNY(allstar) = ny
+	DP_WXOFF(allstar) = xoff
+	DP_WYOFF(allstar) = yoff
+
+	# Update the buffer definition which is currently not used.
+	DP_WLX(allstar) = xoff
+	DP_WMX(allstar) = max (xoff + nx - 1, 0)
+	DP_WLY(allstar) = yoff
+	DP_WMY(allstar) = max (yoff + ny - 1, 0)
+
+	return (DP_WBUF(allstar))
+end
+
+
+# DP_GST -- Return a pointer to the scratch image pixels. The scratch image
+# pointer is assumed to be NULL on the initial all. If the new line limits
+# are totally contained within the old line limits no action is taken and
+# the old line limits are returned. If the scratch image pixels are stored in a 
+# scratch image the pixels must be accessed  sequentially, otherwise an
+# error is returned.
+
+pointer procedure dp_gst (dao, im, line1, line2, mode, flush)
+
+pointer	dao			# pointer to the daophot structure
+pointer	im			# pointer to the input image
+int	line1			# the lower line limit
+int	line2			# the upper line limit
+int	mode			# input / output mode
+int	flush			# flush the buffer
+
+int	nx, ny, xoff, yoff, llast1, llast2
+pointer	allstar
+
+begin
+	allstar = DP_ALLSTAR(dao)
+	if (DP_SBUF(allstar) == NULL) {
+	    llast1 = 0
+	    llast2 = 0
+	} else if (flush == YES) {
+	    line1 = llast1
+	    line2 = llast2
+	} else if ((line1 >= llast1) && (line2 <= llast2)) {
+	    line1 = llast1
+	    line2 = llast2
+	    return (DP_SBUF(allstar))
+	} else if (line1 < llast1)
+	    call error (0,
+		"ERROR: Attempting to access scratch image pixels randomly")
+
+	# All the data is cached.
+	if (DP_CACHE (allstar, A_SUBT) == YES) {
+
+	    nx = IM_LEN(im,1)
+	    ny = IM_LEN(im,2)
+	    xoff = 1
+	    yoff = 1
+	    DP_SBUF(allstar) = DP_SUBT(allstar)
+
+	# Read in some new data.
+	} else if (mode == READ_ONLY) {
+
+            call dp_albufl2r (DP_SUBT(allstar), DP_SBUF(allstar),
+	        llast1, llast2, line1, line2, flush)
+
+	    if (flush == NO) {
+		nx = IM_LEN(im,1)
+	        ny = line2 - line1 + 1
+		xoff = 1
+	        yoff = line1
+	    } else {
+		nx = 0
+		ny = 0
+		xoff = 0
+		yoff = 0
+	    }
+
+	# Write out the old data and read in some new data.
+	} else if (mode == READ_WRITE) {
+
+            call dp_albufl2rw (DP_SUBT(allstar), DP_SUBT(allstar),
+	        DP_SBUF(allstar), llast1, llast2, line1, line2, flush)
+
+	    if (flush == NO) {
+		nx = IM_LEN(im,1)
+	        ny = line2 - line1 + 1
+		xoff = 1
+	        yoff = line1
+	    } else {
+		nx = 0
+		ny = 0
+		xoff = 0
+		yoff = 0
+	    }
+	}
+
+	# Update the required buffer parameters.
+	DP_SNX(allstar) = nx
+	DP_SNY(allstar) = ny
+	DP_SXOFF(allstar) = xoff
+	DP_SYOFF(allstar) = yoff
+
+	# Update the buffer description which is not currently used.
+	DP_SLX(allstar) = xoff
+	DP_SMX(allstar) = max (xoff + nx - 1, 0)
+	DP_SLY(allstar) = yoff
+	DP_SMY(allstar) = max (yoff + ny - 1, 0)
+
+	return (DP_SBUF(allstar))
+end
+
+
+# DP_GDC -- Return a pointer to the subtracted image pixels. The subtracted
+# image pixels pointer is assumed to be NULL on the initial all. If the new
+# line limits are totally contained within the old line limits no action is
+# taken and the old line limits are returned. If the subtracted image pixels are
+# stored in a scratch image the pixels must be accessed  sequentially,
+# otherwise an error is returned.
+
+pointer procedure dp_gdc (dao, im, line1, line2, mode, flush)
+
+pointer	dao			# pointer to the daophot strucuture
+pointer	im			# pointer to the input image
+int	line1, line2		# the upper and lower line limits
+int	mode			# input / output mode
+int	flush			# flush the input data
+
+int	nx, ny, xoff, yoff, llast1, llast2
+pointer	allstar
+
+begin
+	allstar = DP_ALLSTAR(dao)
+	if (DP_DBUF(allstar) == NULL) {
+	    llast1 = 0
+	    llast2 = 0
+	} else if (flush == YES) {
+	    line1 = llast1
+	    line2 = llast2
+	} else if ((line1 >= llast1) && (line2 <= llast2)) {
+	    line1 = llast1
+	    line2 = llast2
+	    return (DP_DBUF(allstar))
+	} else if (line1 < llast1)
+	    call error (0,
+		"ERROR: Attempting to access subtracted image pixels randomly")
+
+	# All the data is cached.
+	if (DP_CACHE (allstar, A_DCOPY) == YES) {
+
+	    nx = IM_LEN(im,1)
+	    ny = IM_LEN(im,2)
+	    xoff = 1
+	    yoff = 1
+	    DP_DBUF(allstar) = DP_DATA(allstar)
+
+	# Read in some new data.
+	} else if (mode == READ_ONLY) {
+
+            call dp_albufl2r (DP_DATA(allstar), DP_DBUF(allstar),
+	        llast1, llast2, line1, line2, flush)
+
+	    if (flush == NO) {
+		nx = IM_LEN(im,1)
+	        ny = line2 - line1 + 1
+		xoff = 1
+	        yoff = line1
+	    } else {
+		nx = 0
+		ny = 0
+		xoff = 0
+		yoff = 0
+	    }
+
+	# Write out the old data and read in some new data.
+	} else if (mode == READ_WRITE) {
+
+            call dp_albufl2rw (DP_DATA(allstar), DP_DATA(allstar),
+	        DP_DBUF(allstar), llast1, llast2, line1, line2, flush)
+
+	    if (flush == NO) {
+		nx = IM_LEN(im,1)
+	        ny = line2 - line1 + 1
+		xoff = 1
+	        yoff = line1
+	    } else { 
+		nx = 0
+		ny = 0
+		xoff = 0
+		yoff = 0
+	    }
+	}
+
+	# Update the required buffer parameters.
+	DP_DNX(allstar) = nx
+	DP_DNY(allstar) = ny
+	DP_DXOFF(allstar) = xoff
+	DP_DYOFF(allstar) = yoff
+
+	# Update the buffer definition which is currently not used.
+	DP_DLX(allstar) = xoff
+	DP_DMX(allstar) = max (xoff + nx - 1, 0)
+	DP_DLY(allstar) = yoff
+	DP_DMY(allstar) = max (yoff + ny - 1, 0)
+
+	return (DP_DBUF(allstar))
+end
+
+
+# DP_ALBUFL2RW -- Maintain a buffer of image lines and which can be read from
+# an input image and written to an output image. The input and output images
+# may be the same but must have the same dimensions. A new buffer is created
+# when the buffer pointer is null and reallocated when the buffer changes size.
+
+procedure dp_albufl2rw (inim, outim, buf, llast1, llast2, line1, line2, flush)
+
+pointer	inim		# the input image pointer
+pointer	outim		# the output image pointer
+pointer	buf		# pointer to the internal buffer
+int	llast1, llast2	# the line limits of the previous buffer
+int	line1, line2	# the line limits of the requested buffer
+int	flush		# flush the output buffer
+
+int	i, ncols, nlines, nllast, lout, nmove
+pointer	buf1, buf2
+pointer	imgl2r(), impl2r()
+
+define	flush_ 11
+
+begin
+	# Write the data in the buffer to the output image and free
+	# the buffer.
+
+	if (flush == YES)
+	    go to flush_
+
+	# Define the size of the current buffer.
+	ncols = IM_LEN(inim,1)
+	nlines = line2 - line1 + 1
+
+	# If the buffer pointer is undefined then allocate memory for the
+	# buffer.
+
+	if (buf == NULL) {
+	    call malloc (buf, ncols * nlines, TY_REAL)
+	    llast1 = 0
+	    llast2 = 0
+	    nllast = 0
+	} else
+	    nllast = llast2 - llast1 + 1
+
+	# Write out the lines that are not needed any more.
+
+	if ((line1 > llast1) && (llast1 > 0)) {
+	    buf2 = buf
+	    lout = min (llast2, line1 - 1)
+	    do i = llast1, lout {
+		buf1 = impl2r (outim, i)
+		call amovr (Memr[buf2], Memr[buf1], ncols)
+		buf2 = buf2 + ncols
+	    }
+	}
+
+	# Now move the remaining lines to the beginning of the buffer.
+
+	nmove = 0
+	if (llast2 >= line1) {
+	    buf1 = buf
+	    buf2 = buf + ncols * (line1 - llast1)
+	    do i = line1, llast2 {
+		if (buf1 != buf2)
+	            call amovr (Memr[buf2], Memr[buf1], ncols)
+		nmove = nmove + 1
+	        buf2 = buf2 + ncols
+	        buf1 = buf1 + ncols
+	    }
+	}
+
+	# Resize the buffer if necessary.
+
+	if (nlines > nllast)
+	    call realloc (buf, ncols * nlines, TY_REAL)
+
+	# Read only the image lines with are different from the last buffer.
+
+	if (line2 > llast2) {
+	    buf1 = buf + ncols * nmove
+	    lout = max (line1, llast2 + 1)
+	    do i = lout, line2 {
+	        buf2 = imgl2r (inim, i)
+	        call amovr (Memr[buf2], Memr[buf1], ncols)
+	        buf1 = buf1 + ncols
+	    }
+	}
+
+	# Save the buffer parameters.
+	llast1 = line1
+	llast2 = line2
+
+	if (flush == NO)
+	    return
+flush_
+	# Flush the remaining portion of the buffer to the output image
+	# and free the buffer space.
+
+	if (buf != NULL) {
+
+	    buf2 = buf
+	    do i = llast1, llast2 {
+	        buf1 = impl2r (outim, i)
+	        call amovr (Memr[buf2], Memr[buf1], ncols)
+	        buf2 = buf2 + ncols
+	    }
+	    call imflush (outim)
+
+	    call mfree (buf, TY_REAL)
+	}
+
+	buf = NULL
+end
+
+
+# DP_ALBUFL2R -- Maintain a buffer of image lines which can be read
+# sequentially from an input image. The buffer pointer must be initialized
+# to NULL. A new buffer is created when the buffer pointer is null and
+# reallocated when the buffer increases in size.
+
+procedure dp_albufl2r (inim, buf, llast1, llast2, line1, line2, flush)
+
+pointer	inim		# pointer to the input image
+pointer	buf		# pointer to the internal buffer
+int	llast1, llast2	# the line limits of the previous buffer
+int	line1, line2	# the line limits of the requested buffer
+int	flush		# flush the output buffer
+
+int	i, ncols, nlines, nllast, lout, nmove
+pointer	buf1, buf2
+pointer	imgl2r()
+
+define	flush_ 11
+
+begin
+
+	# Write the data in the buffer to the output image and free
+	# the buffer.
+
+	if (flush == YES)
+	    go to flush_
+
+	# Define the size parameters for the new buffer.
+
+	ncols = IM_LEN(inim,1)
+	nlines = line2 - line1 + 1
+
+	# If the buffer pointer is undefined then allocate memory for the
+	# buffer.
+
+	if (buf == NULL) {
+	    call malloc (buf, ncols * nlines, TY_REAL)
+	    llast1 = 0
+	    llast2 = 0
+	    nllast = 0
+	} else
+	    nllast = llast2 - llast1 + 1
+	
+	# Now move the remaining lines to the beginning of the buffer.
+
+	nmove = 0
+	if (llast2 >= line1)  {
+	    buf1 = buf
+	    buf2 = buf + ncols * (line1 - llast1)
+	    do i = line1, llast2 {
+		if (buf1 != buf2)
+	            call amovr (Memr[buf2], Memr[buf1], ncols)
+		nmove = nmove + 1
+	        buf2 = buf2 + ncols
+	        buf1 = buf1 + ncols
+	    }
+	}
+
+	# Resize the buffer if necessary.
+
+	if (nlines >  nllast)
+	    call realloc (buf, ncols * nlines, TY_REAL)
+
+	# Read only the image lines with are different from the last buffer.
+
+	if (line2 > llast2) {
+	    buf1 = buf + ncols * nmove
+	    lout = max (line1, llast2 + 1)
+	    do i = lout, line2 {
+	        buf2 = imgl2r (inim, i)
+	        call amovr (Memr[buf2], Memr[buf1], ncols)
+	        buf1 = buf1 + ncols
+	    }
+	}
+
+	# Save the buffer parameters.
+	llast1 = line1
+	llast2 = line2
+
+	if (flush == NO)
+	    return
+flush_
+
+	# Free the buffer space.
+	if (buf != NULL)
+	    call mfree (buf, TY_REAL)
+	buf = NULL
+end
diff --git a/noao/digiphot/daophot/allstar/dpaconfirm.x b/noao/digiphot/daophot/allstar/dpaconfirm.x
new file mode 100644
index 00000000..76e24798
--- /dev/null
+++ b/noao/digiphot/daophot/allstar/dpaconfirm.x
@@ -0,0 +1,37 @@
+include "../lib/daophotdef.h"
+
+# DP_ACONFIRM -- Procedure to confirm the critical allstar parameters.
+
+procedure dp_aconfirm (dao)
+
+pointer	dao		# pointer to the group structure
+
+int	dp_stati()
+
+begin
+	call printf ("\n")
+
+	# Confirm the recentering and sky fitting parameters.
+	call dp_vrecenter (dao)
+	call dp_vgroupsky (dao)
+	call dp_vfitsky (dao)
+	if (dp_stati (dao, FITSKY) == YES) {
+	    call dp_vsannulus (dao)
+	    call dp_vwsannulus (dao)
+	}
+
+	# Confirm the psf radius.
+	call dp_vpsfrad (dao)
+
+	# Confirm the fitting radius.
+	call dp_vfitrad (dao)
+
+	# Confirm the maximum group size.
+	call dp_vmaxgroup (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/allstar/dpalinit.x b/noao/digiphot/daophot/allstar/dpalinit.x
new file mode 100644
index 00000000..4f1afb89
--- /dev/null
+++ b/noao/digiphot/daophot/allstar/dpalinit.x
@@ -0,0 +1,665 @@
+include <mach.h>
+include	<imhdr.h>
+include <math.h>
+include "../lib/daophotdef.h"
+include	"../lib/allstardef.h"
+
+# DP_SETWT -- Initialize the subtracted image and compute the initial weights
+# image from the input image.
+
+procedure dp_setwt (dao, im) 
+
+pointer	dao				# pointer to daophot structure
+pointer	im				# input image decriptor
+
+int	i, ncol, nline
+pointer	sp, v1, v2, v3, v4, line1, line2, line3, line4
+pointer	allstar, wtim, dataim, subtim
+real	rdnoise, mingdata, maxgdata
+int	imgnlr(), impnlr()
+
+begin
+	# Allocate some working space.
+	call smark (sp)
+	call salloc (v1, IM_MAXDIM, TY_LONG)
+	call salloc (v2, IM_MAXDIM, TY_LONG)
+	call salloc (v3, IM_MAXDIM, TY_LONG)
+	call salloc (v4, IM_MAXDIM, TY_LONG)
+
+	# Define the allstar pointer.
+	allstar = DP_ALLSTAR (dao)
+
+	# Define some useful constants.
+	rdnoise = (DP_READNOISE(dao) / DP_PHOTADU(dao)) ** 2
+	if (IS_INDEFR(DP_MINGDATA(dao)))
+	    mingdata = -MAX_REAL
+	else
+	    mingdata = DP_MINGDATA(dao)
+	if (IS_INDEFR(DP_MAXGDATA(dao)))
+	    maxgdata = MAX_REAL
+	else
+	    maxgdata = DP_MAXGDATA(dao)
+
+	wtim = DP_WEIGHTS (allstar)
+	dataim = DP_DATA(allstar)
+	subtim = DP_SUBT(allstar)
+	ncol = IM_LEN (im,1)
+	nline = IM_LEN(im,2)
+
+	call amovkl (long(1), Meml[v1], IM_MAXDIM)
+	call amovkl (long(1), Meml[v2], IM_MAXDIM)
+	call amovkl (long(1), Meml[v3], IM_MAXDIM)
+	call amovkl (long(1), Meml[v4], IM_MAXDIM)
+	do i = 1, nline {
+
+	    # Get the next line
+	    if (imgnlr (im, line1, Meml[v1]) == EOF) 
+		break
+
+	    # Initialize the subtracted image.
+	    if (DP_CACHE(allstar, A_DCOPY) == YES) {
+		call amovr (Memr[line1], Memr[dataim], ncol)
+		dataim = dataim + ncol
+	    } else if (impnlr (dataim, line2, Meml[v2]) != EOF) {
+		call amovr (Memr[line1], Memr[line2], ncol)
+	    }
+
+	    # Initialize the weights image.
+	    if (DP_CACHE(allstar, A_WEIGHT) == YES) {
+		call dp_wtvector (Memr[line1], Memr[wtim], ncol, mingdata,
+		    maxgdata, rdnoise)
+		wtim = wtim + ncol
+	    } else if (impnlr (wtim, line3, Meml[v3]) != EOF) {
+		call dp_wtvector (Memr[line1], Memr[line3], ncol, mingdata,
+		    maxgdata, rdnoise)
+	    }
+
+	    # Initilize the subtracted image.
+	    if (DP_CACHE(allstar, A_WEIGHT) == YES) {
+		;
+	    } else if (impnlr (subtim, line4, Meml[v4]) != EOF) {
+		call amovkr (0.0, Memr[line4], ncol)
+	    }
+
+	}
+
+	# Make sure all the changes are written to disk.
+	if (DP_CACHE(allstar, A_WEIGHT) == NO)
+	    call imflush (DP_WEIGHTS(allstar))
+	if (DP_CACHE(allstar, A_DCOPY) == NO)
+	    call imflush (DP_DATA(allstar))
+	if (DP_CACHE(allstar, A_SUBT) == NO)
+	    call imflush (DP_SUBT(allstar))
+
+	call sfree (sp)
+end
+
+
+# DP_WTVECTOR -- Compute the initial weights for a vector of input data.
+
+procedure dp_wtvector (a, b, ncols, mingdata, maxgdata, rnoisesq)
+
+real	a[ARB]			# input array
+real	b[ARB]			# output array
+int	ncols			# number of points
+real	mingdata		# minimum good data value
+real	maxgdata		# maximum good data value
+real	rnoisesq		# read noise squared in adu
+
+int	i
+
+begin
+	do i = 1, ncols {
+	    if (a[i] < mingdata || a[i] > maxgdata)
+		b[i] = -MAX_REAL
+	    else
+		b[i] = rnoisesq
+	}
+end
+
+
+define	DELTA_MAG		12.5
+define	INIT_REL_BRIGHT		0.003
+
+# DP_ALZERO --  Convert from magnitudes to relative brightnesses before
+# fitting the stars. If the star's magnitude is undefined or more than
+# 12.5 magnitudes fainter than the PSF magnitude, a default relative
+# brightness is defined.
+
+procedure dp_alzero (dao, mag, nstars)
+
+pointer	dao			# pointer to the daophot strucuture
+real	mag[ARB]		# the magnitude array
+int	nstars			# number of stars
+
+int	i
+pointer	psffit
+real	faint
+
+begin
+	psffit = DP_PSFFIT (dao)
+	faint = DP_PSFMAG(psffit) + DELTA_MAG
+
+	do i = 1, nstars {
+	    if (IS_INDEFR(mag[i])) {
+		mag[i] = INIT_REL_BRIGHT
+	    } else if (mag[i] >= faint) {
+		mag[i] = INIT_REL_BRIGHT
+	    } else {
+	        mag[i] = DAO_RELBRIGHT (psffit, mag[i])
+	    }
+	}
+end
+
+
+# DP_STRIP -- Remove stars with undefined centers, stars with centers that
+# are off the input image, and duplicate stars from the input photometry
+# list, where duplicate stars are defined as those within a distance of
+# radius pixels of each other. The duplicate stars are moved to the end of
+# the star list and the number of stars is recomputed.
+
+procedure dp_strip (id, x, y, mag, sky, skip, aier, nstar, sepradsq, ncol,
+	nline, fitradius, verbose)
+
+int	id[ARB]			# array of star ids
+real	x[ARB]			# array of star x values
+real	y[ARB]			# array of star y values
+real	mag[ARB]		# array of star magnitudes
+real	sky[ARB]		# array of star sky values
+int	skip[ARB]		# array of fit/nofit indicators
+int	aier[ARB]		# array of error codes
+int	nstar			# number of stars (may change)
+real	sepradsq		# separation radius squared
+int	ncol			# number of columns in the image
+int	nline			# number of columns in the image
+real	fitradius		# the fitting radius
+int	verbose			# print error messages ?
+
+int	i, j, ier, ahold
+pointer	sp, index, idhold
+real	seprad, dx, dy, xhold, yhold, maghold, skyhold
+
+begin
+	# Duplicate stars are impossible.
+	if (nstar <= 1)
+	    return
+
+	# Allocate some working space.
+	call smark (sp)
+	call salloc (index, nstar, TY_INT)
+
+	# Sort the data in y.
+	call quick (y, nstar, Memi[index], ier)
+
+	# Rectify the remaining arrays.
+	call dp_irectify (id, Memi[index], nstar)
+	call dp_rectify (x, Memi[index], nstar)
+	call dp_rectify (mag, Memi[index], nstar)
+	call dp_rectify (sky, Memi[index], nstar)
+
+	# Determine whether any stars are close to star i.
+	seprad = sqrt (sepradsq)
+	do i = 1, nstar - 1 {
+
+	    # This star was rejected on a previous loop.
+	    if (skip[i] == YES)
+		next
+
+	    # Reject if star has an INDEF valued position or is off image.
+	    if (IS_INDEFR(x[i]) || IS_INDEFR(y[i])) {
+		mag[i] = INDEFR
+		skip[i] = YES
+		aier[i] = ALLERR_OFFIMAGE
+		if (verbose == YES)
+		    call printf (
+		        "REJECTING: Star %d has an undefined x or y value\n")
+			call pargi (id[i])
+	    } else if ((int (x[i] - fitradius) + 1) > ncol) {
+		mag[i] = INDEFR
+		skip[i] = YES
+		aier[i] = ALLERR_OFFIMAGE
+		if (verbose == YES)
+		    call printf ("REJECTING: Star %d is outside the image\n")
+			call pargi (id[i])
+	    } else if (int (x[i] + fitradius) < 1) {
+		mag[i] = INDEFR
+		skip[i] = YES
+		aier[i] = ALLERR_OFFIMAGE
+		if (verbose == YES)
+		    call printf ("REJECTING: Star %d is outside the image\n")
+			call pargi (id[i])
+	    } else if ((int (y[i] - fitradius) + 1) > nline) {
+		mag[i] = INDEFR
+		skip[i] = YES
+		aier[i] = ALLERR_OFFIMAGE
+		if (verbose == YES)
+		    call printf ("REJECTING: Star %d is outside the image\n")
+			call pargi (id[i])
+	    } else if (int (y[i] + fitradius) < 1) {
+		mag[i] = INDEFR
+		skip[i] = YES
+		aier[i] = ALLERR_OFFIMAGE
+		if (verbose == YES)
+		    call printf ("REJECTING: Star %d is outside the image\n")
+			call pargi (id[i])
+	    }
+
+	    # This star was rejected on this loop. 
+	    if (skip[i] == YES)
+		next
+
+	    # Loop over the remaining stars.
+	    do j = i + 1, nstar {
+
+	        # Star has already been rejected on previous loop.
+	        if (skip[j] == YES)
+		    next
+
+		# Test for INDEF.
+	        if (IS_INDEFR(x[j]) || IS_INDEFR(y[j])) {
+		    mag[j] = INDEFR
+		    skip[j] = YES
+		    aier[j] = ALLERR_OFFIMAGE
+		    if (verbose == YES)
+		        call printf (
+		    	"REJECTING: Star %d has an undefined x or y value\n")
+			    call pargi (id[j])
+	        } else if ((int (x[j] - fitradius) + 1) > ncol) {
+		    mag[j] = INDEFR
+		    skip[j] = YES
+		    aier[j] = ALLERR_OFFIMAGE
+		    if (verbose == YES)
+		        call printf (
+			    "REJECTING: Star %d is outside the image\n")
+			    call pargi (id[j])
+	        } else if (int (x[j] + fitradius) < 1) {
+		    mag[j] = INDEFR
+		    skip[j] = YES
+		    aier[j] = ALLERR_OFFIMAGE
+		    if (verbose == YES)
+		        call printf (
+			    "REJECTING: Star %d is outside the image\n")
+			    call pargi (id[j])
+	        } else if ((int (y[j] - fitradius) + 1) > nline) {
+		    mag[j] = INDEFR
+		    skip[j] = YES
+		    aier[j] = ALLERR_OFFIMAGE
+		    if (verbose == YES)
+		        call printf (
+			    "REJECTING: Star %d is outside the image\n")
+			    call pargi (id[j])
+	        } else if (int (y[j] + fitradius) < 1) {
+		    mag[j] = INDEFR
+		    skip[j] = YES
+		    aier[j] = ALLERR_OFFIMAGE
+		    if (verbose == YES)
+		        call printf (
+			    "REJECTING: Star %d is outside the image\n")
+			    call pargi (id[j])
+	        }
+
+		# Star was rejected.
+		if (skip[j] == YES)
+		    next
+
+		# Test for proximity.
+		dy = y[j] - y[i]
+		if (dy > seprad)
+		    break
+		dx = x[j] - x[i]
+		if (abs (dx) > seprad)
+		    next
+		if ((dx * dx + dy * dy) > sepradsq)
+		    next
+
+		# Set the magnitude of the star to INDEF and skip.
+		if (mag[j] <= mag[i]) {
+		    mag[j] = INDEFR
+		    skip[j] = YES
+		    aier[j] = ALLERR_MERGE
+		    if (verbose == YES) {
+		        call printf (
+			    "REJECTING: Star %d has merged with star %d\n")
+			    call pargi (id[j])
+		            call pargi (id[i])
+		    }
+		} else {
+		    mag[i] = INDEFR
+		    skip[i] = YES
+		    aier[i] = ALLERR_MERGE
+		    if (verbose == YES) {
+		        call printf (
+			    "REJECTING: Star %d has merged with star %d\n")
+			    call pargi (id[i])
+			    call pargi (id[j])
+		    }
+		    break
+		}
+	    }
+	}
+
+	# Remove the duplicate stars.
+	for (i = 1; i <= nstar; i = i + 1) {
+
+	    # Redefine the number of stars by removing skipped stars from
+	    # the end of the star list.
+	    while (nstar >= 1) {
+		if (skip[nstar] == NO)
+		    break
+		nstar = nstar - 1
+	    }
+	    if (i > nstar)
+		break
+	    if (skip[i] == NO)
+		next
+
+	    # Switch the rejected star with the one at the end of the list.
+	    idhold = id[i]
+	    xhold = x[i]
+	    yhold = y[i]
+	    maghold = mag[i]
+	    skyhold = sky[i]
+	    ahold = aier[i]
+
+	    id[i] = id[nstar]
+	    x[i] = x[nstar]
+	    y[i] = y[nstar]
+	    mag[i] = mag[nstar]
+	    sky[i] = sky[nstar]
+	    skip[i] = NO
+	    aier[i] = aier[nstar]
+
+	    id[nstar] = idhold
+	    x[nstar] = xhold
+	    y[nstar] = yhold
+	    mag[nstar] = maghold
+	    sky[nstar] = skyhold
+	    skip[nstar] = YES
+	    aier[nstar] = ahold
+
+	    nstar = nstar - 1
+	}
+
+	call sfree (sp)
+end
+
+
+# DP_WSTINIT -- Initialize the weight and scratch arrays / images.
+
+procedure dp_wstinit (dao, im, xcen, ycen, mag, nstar, radius, x1, x2, y1, y2)
+
+pointer	dao				# pointer to the daophot structure
+pointer	im				# pointer to the input image
+real	xcen[ARB]			# the x centers array
+real	ycen[ARB]			# the y centers array
+real	mag[ARB]			# the magnitude array
+int	nstar				# the number of stars
+real	radius				# radius for subraction
+int	x1, x2				# column limits
+int	y1, y2				# line limits
+
+int	j, l, ncol, npix, nl, ns, lx, mx
+pointer psffit, allstar, databuf, wtbuf, owtbuf, subtbuf
+real	rsq, psfradius, psfradsq, x, y, dy, dysq, deltax, deltay
+pointer	imgs2r(), imps2r()
+
+begin
+	# Set up some pointers.
+	psffit = DP_PSFFIT(dao)
+	allstar = DP_ALLSTAR(dao)
+
+	# Set up some constants.
+	ncol = IM_LEN(im,1)
+	npix = x2 - x1 + 1
+	rsq = radius ** 2
+	if (DP_PSFSIZE(psffit) == 0)
+	    psfradius = DP_PSFRAD(dao)
+	else
+	    psfradius = (real (DP_PSFSIZE(psffit) - 1) / 2.0 - 1.0) / 2.0
+	psfradsq = psfradius * psfradius
+
+	# Begin the search of the groups at the first group.
+	ns = 0
+	nl = 0
+
+	ns = 0
+	do j = y1, y2 {
+
+	    # Get the data.
+	    if (DP_CACHE(allstar,A_DCOPY) == YES)
+		databuf = DP_DATA(allstar) + (j - 1) * ncol + x1 - 1
+	    else
+		databuf = imgs2r (DP_DATA(allstar), x1, x2, j, j)
+	    if (DP_CACHE(allstar,A_WEIGHT) == YES) {
+		wtbuf = DP_WEIGHTS(allstar) + (j - 1) * ncol + x1 - 1
+		owtbuf = wtbuf
+	    } else {
+		owtbuf = imps2r (DP_WEIGHTS(allstar), x1, x2, j, j)
+		wtbuf = imgs2r (DP_WEIGHTS(allstar), x1, x2, j, j)
+	    }
+	    if (DP_CACHE(allstar,A_SUBT) == YES)
+		subtbuf = DP_SUBT(allstar) + (j - 1) * ncol + x1 - 1
+	    else
+		subtbuf = imps2r (DP_SUBT(allstar), x1, x2, j, j)
+
+	    # Set all the weights in the working array negative.
+	    call aabsr (Memr[wtbuf], Memr[owtbuf], npix)
+	    call anegr (Memr[owtbuf], Memr[owtbuf], npix)
+	    call amovkr (0.0, Memr[subtbuf], npix)
+
+	    # Set the y coordinate.
+	    y = real (j)
+
+	    # Initialize the weight and scratch arrays.
+
+	    # Find all the points within one working radius of each star.
+	    # Set all the sigmas positive again and copy them from DATA
+	    # to SUBT.
+	    do l = ns + 1, nstar {
+		dy = y - ycen[l] 
+		if (dy > radius) {
+		    ns = l
+		    next
+		} else if (dy < -radius)
+		    break 
+		dysq = dy ** 2
+		lx = max (x1, min (x2, int (xcen[l] - radius) + 1))
+		mx = max (x1, min (x2, int (xcen[l] + radius)))
+		x = xcen[l] - lx + 1.0
+		lx = lx - x1 + 1
+		mx = mx - x1 + 1
+		call dp_wstvector (x, Memr[databuf+lx-1], Memr[owtbuf+lx-1],
+		    Memr[subtbuf+lx-1], mx - lx + 1, dysq, rsq, -MAX_REAL)
+	    }
+
+	    do l = nl + 1, nstar {
+		dy = y - ycen[l] 
+		if (dy > psfradius) {
+		    nl = l
+		    next
+		} else if (dy < -psfradius)
+		    break 
+		dysq = dy ** 2
+		lx = max (x1, min (x2, int (xcen[l] - psfradius) + 1))
+		mx = max (x1, min (x2, int (xcen[l] + psfradius)))
+		x = xcen[l] - lx + 1.0
+		lx = lx - x1 + 1
+		mx = mx - x1 + 1
+		call dp_wpsf (dao, im, xcen[l], ycen[l], deltax, deltay, 1)
+		deltax = (deltax - 1.0) / DP_PSFX(psffit) - 1.0
+		deltay = (deltay - 1.0) / DP_PSFY(psffit) - 1.0
+		call dp_alsubstar (psffit, x, dy, deltax, deltay, mag[l],
+		    Memr[subtbuf+lx-1], Memr[owtbuf+lx-1], mx - lx + 1, dysq,
+		    psfradsq)
+	    }
+	}
+
+	if (DP_CACHE(allstar,A_WEIGHT) == NO)
+	    call imflush (DP_WEIGHTS(allstar))
+	if (DP_CACHE(allstar,A_SUBT) == NO)
+	    call imflush (DP_SUBT(allstar))
+end
+
+
+# DP_WSTVECTOR -- Set the weight and scratch vector
+
+procedure dp_wstvector (xcen, data, weight, subt, npix, dysq, rsq, badwt)
+
+real	xcen			# x coordinate of center
+real	data[ARB]		# the data array
+real	weight[ARB]		# the weight array
+real	subt[ARB]		# the subtracted array array
+int	npix			# the number of pixels
+real	dysq			# the y distance squared
+real	rsq			# the inclusion radius squared
+real	badwt			# badwt value
+
+int	i
+real	dx
+
+begin
+	do i = 1, npix {
+	    dx = (real (i) - xcen)
+	    if ((dx ** 2 + dysq) <= rsq) {
+		if (weight[i] <= badwt)
+		    next
+		weight[i] = abs (weight[i])
+		subt[i] = data[i]
+	    } else if (dx >= 0.0)
+		break
+	}
+end
+
+
+# DP_ALSUBSTAR -- The subtraction vector.
+
+procedure dp_alsubstar (psffit, xcen, dy, deltax, deltay, mag, subt, weight,
+	npix, dysq, psfradsq)
+
+pointer	psffit			# pointer to the psf fitting structure
+real	xcen			# x coordinate of center
+real	dy			# y offset from psf center
+real	deltax			# x distance from psf position
+real	deltay			# y distance from psf position
+real	mag			# the magnitude of the star
+real	subt[ARB]		# the subtracted array array
+real	weight[ARB]		# the weight array
+int	npix			# the number of pixels
+real	dysq			# the y distance squared
+real	psfradsq		# the inclusion radius squared
+
+int	i
+real	dx, dvdxc, dvdyc
+real	dp_usepsf()
+
+begin
+	do i = 1, npix {
+	    if (weight[i] < 0.0)
+		next
+	    dx = real (i) - xcen
+	    if ((dx ** 2 + dysq) < psfradsq) {
+		subt[i] = subt[i] - mag * 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), deltax, deltay,
+		    dvdxc, dvdyc)
+	    } else if (dx > 0.0)
+		break
+	}
+end
+
+
+# DP_ALSKY -- Recompute the sky values.
+
+procedure dp_alsky (dao, im, xcen, ycen, sky, nstar, x1, x2, y1, y2, rinsky,
+	routsky, minnsky, badwt)
+
+pointer	dao			# pointer to the daophot structure
+int	im			# the input image pointer
+real	xcen[ARB]		# x coordinates
+real	ycen[ARB]		# y coordinates
+real	sky[ARB]		# sky values
+int	nstar			# number of stars
+int	x1, x2			# coordinate limits
+int	y1, y2			# coordinate limits
+real	rinsky			# inner radius of the sky annulus
+real	routsky			# outer radius of the sky annulus
+int	minnsky			# minimum number of sky pixels
+real	badwt			# the bad weight value
+
+int	istar, i, j, ncols, lenskybuf, lx, mx, ly, my, line1, line2, nsky, ier
+pointer	allstar, sub, psub, wgt, pwgt, skyvals, index
+real	risq, rosq, dannulus, dysq, dx, rsq
+pointer	dp_gst(), dp_gwt()
+
+begin
+	# Get the allstar pointer.
+	allstar = DP_ALLSTAR(dao)
+
+	# Define some constants
+	ncols = IM_LEN(im,1)
+	risq = rinsky ** 2 
+	rosq = routsky ** 2
+	dannulus = routsky - rinsky
+
+	# Allcoate some working memory.
+	lenskybuf = PI * (2.0 * rinsky + dannulus + 1.0) * (dannulus + 0.5)
+	call malloc (skyvals, lenskybuf, TY_REAL)
+	call malloc (index, lenskybuf, TY_INT)
+
+	# Accumulate the sky buffer.
+	do istar = 1, nstar {
+
+	    # Get the data.
+	    lx = max (x1, min (x2, int (xcen[istar] - routsky) + 1))
+	    mx = max (x1, min (x2, int (xcen[istar] + routsky)))
+	    ly = max (y1, min (y2, int (ycen[istar] - routsky) + 1))
+	    my = max (y1, min (y2, int (ycen[istar] + routsky)))
+	    line1 = ly
+	    line2 = my
+	    sub = dp_gst (dao, im, line1, line2, READ_ONLY, NO)
+	    wgt = dp_gwt (dao, im, line1, line2, READ_ONLY, NO)
+
+	    nsky = 0
+	    psub = sub + (ly - DP_SYOFF(allstar)) * ncols 
+	    pwgt = wgt + (ly - DP_WYOFF(allstar)) * ncols
+	    do j = ly, my {
+		dysq = (real (j) - ycen[istar]) ** 2
+		do i = lx, mx {
+		    dx = (real (i) - xcen[istar])
+		    rsq = dx ** 2 + dysq
+		    if (rsq > rosq) {
+			if (dx > 0.0)
+			    break
+			else 
+			    next
+		    }
+		    if (rsq < risq)
+			next
+		    if (Memr[pwgt+i-1] <= badwt)
+			next
+		    Memr[skyvals+nsky] = Memr[psub+i-1]
+		    nsky = nsky + 1
+		}
+		psub = psub + ncols
+		pwgt = pwgt + ncols
+	    }
+
+	    # Compute the new sky value.
+	    if (nsky > minnsky) {
+	        call quick (Memr[skyvals], nsky, Memi[index], ier)
+	        j = nint (0.2 * nsky)
+	        dx = 0.0
+	        do i = (nsky + 1) / 2 - j, (nsky / 2) + j + 1
+		    dx = dx + Memr[skyvals+i-1]
+	        sky[istar] = dx / real ((nsky / 2) + 2 * j + 2 -
+		    (nsky + 1) / 2) 
+	    }
+	}
+
+
+	call mfree (skyvals, TY_REAL)
+	call mfree (index, TY_INT)
+	sub = dp_gst (dao, im, line1, line2, READ_ONLY, YES)
+	wgt = dp_gwt (dao, im, line1, line2, READ_ONLY, YES)
+end
diff --git a/noao/digiphot/daophot/allstar/dpalmemstar.x b/noao/digiphot/daophot/allstar/dpalmemstar.x
new file mode 100644
index 00000000..fd9ee80e
--- /dev/null
+++ b/noao/digiphot/daophot/allstar/dpalmemstar.x
@@ -0,0 +1,182 @@
+include "../lib/daophotdef.h"
+include "../lib/allstardef.h"
+
+# DP_ALLSTARSETUP -- Procedure to set up the ALLSTAR parameters.
+
+procedure dp_allstarsetup (dp)
+
+pointer	dp		# pointer to daophot structure
+
+pointer	allstar
+
+begin
+	# Allocate memory.
+	call malloc (DP_ALLSTAR(dp), LEN_ALLSTARSTRUCT, TY_STRUCT)
+	allstar = DP_ALLSTAR(dp)
+
+	DP_DATA(allstar) = NULL
+	DP_SUBT(allstar) = NULL
+	DP_WEIGHTS(allstar) = NULL
+	DP_DBUF(allstar) = NULL
+	DP_SBUF(allstar) = NULL
+	DP_WBUF(allstar) = NULL
+
+	DP_ANUMER(allstar) = NULL
+	DP_ADENOM(allstar) = NULL
+	DP_ASUMWT(allstar) = NULL
+	DP_ARPIXSQ(allstar) = NULL
+	DP_ASKIP(allstar) = NULL
+	DP_AXCLAMP(allstar) = NULL
+	DP_AYCLAMP(allstar) = NULL
+	DP_AXOLD(allstar) = NULL
+	DP_AYOLD(allstar) = NULL
+	DP_AX(allstar) = NULL
+	DP_AV(allstar) = NULL
+	DP_AC(allstar) = NULL
+	DP_ALAST(allstar) = NULL
+	DP_ANPIX(allstar) = NULL
+	DP_AIER(allstar) = NULL
+end
+
+
+# DP_ALMEMSTAR -- Procedure to allocate sufficient memory for ALLSTAR.
+
+procedure dp_almemstar (dao, max_star, max_group)
+
+pointer	dao			# pointer to daophot structure
+int	max_star		# maximum number of stars
+int	max_group		# maximum group size
+
+pointer	allstar
+
+begin
+	allstar = DP_ALLSTAR(dao)
+
+	if (DP_ANUMER (allstar) != NULL)
+	    call mfree (DP_ANUMER (allstar), TY_REAL)
+	call malloc (DP_ANUMER (allstar), max_star, TY_REAL)
+
+	if (DP_ADENOM (allstar) != NULL)
+	    call mfree (DP_ADENOM (allstar), TY_REAL)
+	call malloc (DP_ADENOM (allstar), max_star, TY_REAL)
+
+	if (DP_ARPIXSQ (allstar) != NULL)
+	    call mfree (DP_ARPIXSQ (allstar), TY_REAL)
+	call malloc (DP_ARPIXSQ (allstar), max_star, TY_REAL)
+
+	if (DP_ASUMWT (allstar) != NULL)
+	    call mfree (DP_ASUMWT (allstar), TY_REAL)
+	call malloc (DP_ASUMWT (allstar), max_star, TY_REAL)
+
+	if (DP_ANPIX (allstar) != NULL)
+	    call mfree (DP_ANPIX (allstar), TY_INT)
+	call malloc (DP_ANPIX (allstar), max_star, TY_INT)
+
+	if (DP_AIER (allstar) != NULL)
+	    call mfree (DP_AIER (allstar), TY_INT)
+	call malloc (DP_AIER (allstar), max_star, TY_INT)
+
+	if (DP_ASKIP (allstar) != NULL)
+	    call mfree (DP_ASKIP (allstar), TY_INT)
+	call malloc (DP_ASKIP (allstar), max_star, TY_INT)
+
+	if (DP_ALAST (allstar) != NULL)
+	    call mfree (DP_ALAST (allstar), TY_INT)
+	call malloc (DP_ALAST (allstar), max_star, TY_INT)
+
+	if (DP_AXCLAMP (allstar) != NULL)
+	    call mfree (DP_AXCLAMP (allstar), TY_REAL)
+	call malloc (DP_AXCLAMP (allstar), max_star, TY_REAL)
+
+	if (DP_AYCLAMP (allstar) != NULL)
+	    call mfree (DP_AYCLAMP (allstar), TY_REAL)
+	call malloc (DP_AYCLAMP (allstar), max_star, TY_REAL)
+
+	if (DP_AXOLD (allstar) != NULL)
+	    call mfree (DP_AXOLD (allstar), TY_REAL)
+	call malloc (DP_AXOLD (allstar), max_star, TY_REAL)
+
+	if (DP_AYOLD (allstar) != NULL)
+	    call mfree (DP_AYOLD (allstar), TY_REAL)
+	call malloc (DP_AYOLD (allstar), max_star, TY_REAL)
+
+        # Allocate space for the fitting matrices. Note that nine
+        # times less space is required if recentering is turned
+        # off.
+        if (DP_RECENTER(dao) == YES) {
+
+	    if (DP_AX (allstar) != NULL)
+	        call mfree (DP_AX (allstar), TY_REAL)
+	    call malloc (DP_AX (allstar), 3 * max_group + 1, TY_REAL)
+
+	    if (DP_AV (allstar) != NULL)
+	        call mfree (DP_AV (allstar), TY_REAL)
+	    call malloc (DP_AV (allstar), 3 * max_group + 1, TY_REAL)
+
+	    if (DP_AC (allstar) != NULL)
+	        call mfree (DP_AC (allstar), TY_REAL)
+	    call malloc (DP_AC (allstar), (3 * max_group + 1) *
+	        (3 * max_group + 1), TY_REAL)
+
+        } else {
+
+            if (DP_AX (allstar) != NULL)
+                call mfree (DP_AX (allstar), TY_REAL)
+            call malloc (DP_AX (allstar), max_group + 1, TY_REAL)
+
+            if (DP_AV (allstar) != NULL)
+                call mfree (DP_AV (allstar), TY_REAL)
+            call malloc (DP_AV (allstar), max_group + 1, TY_REAL)
+
+            if (DP_AC (allstar) != NULL)
+                call mfree (DP_AC (allstar), TY_REAL)
+            call malloc (DP_AC (allstar), (max_group + 1) * (max_group + 1),
+                TY_REAL)
+        }
+end
+
+
+# DP_ALCLOSE -- Procedure to close up the ALLSTAR parameters.
+
+procedure dp_alclose (dp)
+
+pointer	dp		# pointer to daophot structure
+
+pointer	allstar
+
+begin
+	allstar = DP_ALLSTAR(dp)
+
+	if (DP_ANUMER (allstar) != NULL)
+	    call mfree (DP_ANUMER (allstar), TY_REAL)
+	if (DP_ADENOM (allstar) != NULL)
+	    call mfree (DP_ADENOM (allstar), TY_REAL)
+	if (DP_ARPIXSQ (allstar) != NULL)
+	    call mfree (DP_ARPIXSQ (allstar), TY_REAL)
+	if (DP_ASUMWT (allstar) != NULL)
+	    call mfree (DP_ASUMWT (allstar), TY_REAL)
+	if (DP_ASKIP (allstar) != NULL)
+	    call mfree (DP_ASKIP (allstar), TY_INT)
+	if (DP_ALAST (allstar) != NULL)
+	    call mfree (DP_ALAST (allstar), TY_INT)
+	if (DP_AXCLAMP (allstar) != NULL)
+	    call mfree (DP_AXCLAMP (allstar), TY_REAL)
+	if (DP_AYCLAMP (allstar) != NULL)
+	    call mfree (DP_AYCLAMP (allstar), TY_REAL)
+	if (DP_AXOLD (allstar) != NULL)
+	    call mfree (DP_AXOLD (allstar), TY_REAL)
+	if (DP_AYOLD (allstar) != NULL)
+	    call mfree (DP_AYOLD (allstar), TY_REAL)
+	if (DP_AX (allstar) != NULL)
+	    call mfree (DP_AX (allstar), TY_REAL)
+	if (DP_AV (allstar) != NULL)
+	    call mfree (DP_AV (allstar), TY_REAL)
+	if (DP_AC (allstar) != NULL)
+	    call mfree (DP_AC (allstar), TY_REAL)
+	if (DP_ANPIX (allstar) != NULL)
+	    call mfree (DP_ANPIX (allstar), TY_INT)
+	if (DP_AIER (allstar) != NULL)
+	    call mfree (DP_AIER (allstar), TY_INT)
+
+	call mfree (allstar, TY_STRUCT)
+end
diff --git a/noao/digiphot/daophot/allstar/dpalphot.x b/noao/digiphot/daophot/allstar/dpalphot.x
new file mode 100644
index 00000000..74741958
--- /dev/null
+++ b/noao/digiphot/daophot/allstar/dpalphot.x
@@ -0,0 +1,1438 @@
+include <mach.h>
+include	<imhdr.h>
+include "../lib/daophotdef.h"
+include "../lib/apseldef.h"
+include	"../lib/allstardef.h"
+
+# DP_ALPHOT -- Perform one iteration on each group of stars.
+
+int procedure dp_alphot (dao, im, ntot, istar, niter, sepcrit, sepmin, wcrit,
+	clip, clampmax, version)
+
+pointer	dao			# pointer to the daophot structure
+pointer	im			# pointer to the input image
+int	ntot			# the total number of stars left unfit
+int	istar			# first star in the current group
+int	niter			# the current iteration
+real	sepcrit			# critical separation for merging
+real	sepmin			# minimum separation for merging
+real	wcrit			# the critical error for merging
+bool	clip			# clip the data
+real	clampmax		# maximum clamping factor
+int	version			# version number
+
+bool	regroup
+int	cdimen, lstar, nstar, fstar, nterm, ixmin, ixmax, iymin, iymax
+int	fixmin, fixmax, fiymin, fiymax, nxpix, nypix, flag, i, j, k, ifaint
+pointer	apsel, psffit, allstar, data, subt, weights
+real	radius, totradius, mean_sky, pererr, peakerr, faint
+real	xmin, xmax, ymin, ymax, sumres, grpwt, chigrp
+
+bool	dp_alredo(), dp_checkc(), dp_almerge(), dp_alclamp(), dp_alfaint()
+int	dp_laststar(), dp_delfaintest()
+pointer	dp_gwt(), dp_gst(), dp_gdc()
+real	dp_almsky(), asumr()
+
+begin
+	# Get some daophot pointers.
+	psffit = DP_PSFFIT (dao)
+	apsel = DP_APSEL(dao)
+	allstar = DP_ALLSTAR (dao)
+
+	# Define some constants. At some point these should be stored
+	# in the allstar structure at task initialization. When the final
+	# rewrite gets done this will occur.
+
+	radius = DP_FITRAD(dao)
+	totradius = DP_FITRAD(dao) + DP_PSFRAD(dao) + 1
+	if (DP_RECENTER(dao) == YES)
+	    cdimen = 3 * DP_MAXGROUP(dao) + 1
+	else
+	    cdimen = DP_MAXGROUP(dao) + 1
+	pererr = 0.01 * DP_FLATERR(dao)
+	peakerr = 0.01 * DP_PROFERR(dao) / (Memr[DP_PSFPARS(psffit)] *
+	    Memr[DP_PSFPARS(psffit)+1])
+	regroup = false
+
+	repeat {
+
+	    # Given the last star arrary as computed by dp_regroup, the first
+	    # star in the current group, and the total number of stars,
+	    # find the last star in the current group and compute the
+	    # total number of stars in the group. 
+
+	    lstar = dp_laststar (Memi[DP_ALAST(allstar)], istar, ntot)
+	    nstar = lstar - istar + 1
+
+	    # Don't compute the subraster limits if no regroup has been
+	    # performed.
+	    if (! regroup) {
+	        call alimr (Memr[DP_APXCEN(apsel)+istar-1], nstar, xmin, xmax)	
+	        call alimr (Memr[DP_APYCEN(apsel)+istar-1], nstar, ymin, ymax)	
+	    }
+
+	    # Determine whether the group is too large to be fit.
+
+	    if (nstar > DP_MAXGROUP (dao)) {
+
+	        if (DP_VERBOSE(dao) == YES) {
+	            call printf (
+		    "Group too large: %5d  Fitting Radius: %5.2f  Regrouping\n")
+		        call pargi (nstar)
+		        call pargr (radius)
+	        }
+
+	        # If size of group too small go to next group.
+	        radius = RADIUS_FRACTION * radius
+	        if ((DP_RECENTER(dao) == YES) && (niter >= 2)) {
+	            if (radius < DENSE_RADIUS1) {
+			fstar = dp_delfaintest (Memr[DP_APMAG(apsel)], istar,
+			    lstar)
+		        Memr[DP_APMAG(apsel)+fstar-1] = INDEFR
+			Memi[DP_ASKIP(allstar)+fstar-1] = YES
+			Memi[DP_AIER(allstar)+fstar-1] = ALLERR_BIGGROUP
+		        if (DP_VERBOSE(dao) == YES) {
+		            call printf (
+		 "REJECTING: Star ID: %d Group too dense to reduce\n")
+				call pargi (Memi[DP_APID(apsel)+fstar-1])
+			}
+		        return (lstar)
+	            }
+	        } else {
+	            if (radius < DENSE_RADIUS2) {
+			fstar = dp_delfaintest (Memr[DP_APMAG(apsel)], istar,
+			    lstar)
+		        Memr[DP_APMAG(apsel)+fstar-1] = INDEFR
+			Memi[DP_ASKIP(allstar)+fstar-1] = YES
+			Memi[DP_AIER(allstar)+fstar-1] = ALLERR_BIGGROUP
+		        if (DP_VERBOSE(dao) == YES) {
+		            call printf (
+			"REJECTING: Star ID: %d Group too dense to reduce\n")
+				call pargi (Memi[DP_APID(apsel)+fstar-1])
+			}
+			return (lstar)
+	            }
+	        }
+
+	        # Regroup the stars.
+	        call dp_regroup (Memi[DP_APID(apsel)+istar-1],
+	            Memr[DP_APXCEN(apsel)+istar-1], Memr[DP_APYCEN(apsel)+
+		    istar-1], Memr[DP_APMAG(apsel)+istar-1],
+		    Memr[DP_APMSKY(apsel)+istar-1],
+		    Memr[DP_ASUMWT(allstar)+istar-1],
+		    Memr[DP_AXOLD(allstar)+ istar-1],
+		    Memr[DP_AYOLD(allstar)+istar-1],
+		    Memr[DP_AXCLAMP(allstar)+istar-1],
+		    Memr[DP_AYCLAMP(allstar)+istar-1], nstar, radius,
+		    Memi[DP_ALAST(allstar)+istar-1])
+		regroup = true
+	        next
+	    }
+
+	    # Compute the mean sky for the group. If the sky is undefined
+	    # reject the group.
+	    mean_sky = dp_almsky (Memr[DP_APMSKY(apsel)+istar-1], nstar)
+	    if (IS_INDEFR(mean_sky)) {
+		call amovkr (INDEFR, Memr[DP_APMAG(apsel)+istar-1], nstar)
+		call amovki (YES, Memi[DP_ASKIP(allstar)+istar-1], nstar)
+		call amovki (ALLERR_INDEFSKY, Memi[DP_AIER(allstar)+istar-1],
+		    nstar)
+		if (DP_VERBOSE(dao) == YES) {
+		    do i = istar, lstar {
+		        call printf (
+			"REJECTING: Star ID: %d Group sky value is undefined\n")
+			    call pargi (Memi[DP_APID(apsel)+i-1])
+		    }
+		}
+		 return (lstar)
+	    }
+
+	    # Re-compute the number of terms in the fitting equation.
+	    if ((DP_RECENTER(dao) == YES) && (niter >= 2))
+	        nterm = 3 * nstar
+	    else
+	        nterm = nstar
+
+	    # Zero the fitting arrays.
+	    chigrp = asumr (Memr[DP_ASUMWT(allstar)+istar-1], nstar) / nstar
+	    call aclrr (Memr[DP_APCHI(apsel)+istar-1], nstar)
+	    call aclrr (Memr[DP_ASUMWT(allstar)+istar-1], nstar)
+	    call aclrr (Memr[DP_ANUMER(allstar)+istar-1], nstar)
+	    call aclrr (Memr[DP_ADENOM(allstar)+istar-1], nstar)
+	    call aclri (Memi[DP_ANPIX(allstar)+istar-1], nstar)
+	    call aclrr (Memr[DP_AV(allstar)], nterm)
+	    call aclrr (Memr[DP_AC(allstar)], cdimen * cdimen)
+
+	    # Compute the subraster limits.
+
+	    ixmin = min (IM_LEN (im,1), max (1, int (xmin - totradius) + 1))
+	    ixmax = min (IM_LEN (im,1), max (1,int (xmax + totradius)))
+	    iymin = min (IM_LEN (im,2), max (1, int (ymin - totradius) + 1))
+	    iymax = min (IM_LEN (im,2), max (1, int (ymax + totradius)))
+
+	    # Get pointer to the required weight, scratch image and
+	    # subtracted image pixels. Need to modify this so writing
+	    # is only enabled if iter >= MIN_ITER.
+
+	    subt = dp_gst (dao, im, iymin, iymax, READ_ONLY, NO)
+	    weights = dp_gwt (dao, im, iymin, iymax, READ_WRITE, NO)
+	    data = dp_gdc (dao, im, iymin, iymax, READ_WRITE, NO)
+
+	    # Compute the fitting limits in the subraster.
+
+	    fixmin = min (ixmax, max (ixmin, int (xmin - DP_FITRAD(dao)) + 1))
+	    fixmax = min (ixmax, max (ixmin, int (xmax + DP_FITRAD(dao))))
+	    fiymin = min (iymax, max (iymin, int (ymin - DP_FITRAD(dao)) + 1))
+	    fiymax = min (iymax, max (iymin, int (ymax + DP_FITRAD(dao))))
+	    nypix = fiymax - fiymin + 1
+	    nxpix = fixmax - fixmin + 1
+
+	    # Now we deal with the pixels one by one.
+	    sumres = 0.0
+	    grpwt = 0.0
+
+	    # Accumulate the data into the matrix.
+	    call dp_alaccum (dao, im, Memr[data], DP_DNX(allstar),
+	        DP_DNY(allstar), DP_DXOFF(allstar), DP_DYOFF(allstar),
+		Memr[subt], DP_SNX(allstar), DP_SNY(allstar), DP_SXOFF(allstar),
+		DP_SYOFF(allstar), Memr[weights], DP_WNX(allstar),
+		DP_WNY(allstar), DP_WXOFF(allstar), DP_WYOFF(allstar), nxpix,
+		nypix, fixmin, fiymin, mean_sky, istar, lstar, niter, clip,
+		pererr, peakerr, sumres, grpwt, chigrp, cdimen, nterm)
+
+	    # Reflect the normal matrix across the diagonal. 
+
+	    call dp_mreflect (Memr[DP_AC(allstar)], cdimen, nterm)
+
+	    # Compute the estimate of the standard deviation of the residuals
+	    # for the group as a whole, and for each star. This estimate starts
+	    # out as SQRT (PI / 2) * [sum [weight * ABS (residual / sigma)] /
+	    # sum [weight] # and then gets corrected for bias by SQRT (# of
+	    # pixels / [# of pixel - # degrees of freedom]).
+	    #
+	    # But now we drive the value toward unity, depending upon exactly
+	    # how many pixels were involved. If CHI is based on exactly a total
+	    # weight of 3, then it is extremely poorly determined, and we just 
+	    # want to keep CHI = 1. The larger the GRPWT is, the better 
+	    # determined CHI is, and the less we want to force it toward unity.
+	    # So, just take the weighted average of CHI and unity, with weights
+	    # GRPWT = 3, and 1, respectively.
+
+	    if (grpwt > nterm) {
+	        chigrp = CHI_NORM * sumres / sqrt (grpwt * (grpwt - nterm))
+	        chigrp = ((grpwt - 3.0) * chigrp + 3.0) / grpwt
+	    } else {
+	        chigrp = 1.0
+	    }
+
+	    # CHIGRP has been pulled towards its expected value of unity to
+	    # keep the statistics of a small number of pixels from completely
+	    # dominating the error analysis. Similarly, the photometric errors
+	    # for the individual stars will be pulled toward unity now. Later
+	    # on, if the number of stars in the group is greater than one,
+	    # CHI(I) will nudged toward the group average. In order to work
+	    # optimally this requires that the gain and read noise and any
+	    # other parameters which represent the noise sources have been
+	    # properly specified.
+	    #
+	    # At the same time, make sure that every star in the group contains
+	    # at least MIN_FIT_PIXEL valid pixels if re-centroiding is being 
+	    # performed, 1 valid pixel if not. If any star in the group fails
+	    # to meet this criterion, mark that star for deletion and skip
+	    # ahead to the next group.
+
+	    if (dp_alredo (Memi[DP_APID(apsel)], Memr[DP_APMAG(apsel)],
+	        Memr[DP_APCHI(apsel)], Memr[DP_ASUMWT(allstar)],
+	        Memi[DP_ANPIX(allstar)], Memi[DP_ASKIP(allstar)],
+		Memi[DP_AIER(allstar)], istar, lstar, niter, chigrp,
+		DP_RECENTER(dao), DP_VERBOSE(dao)))
+		return (lstar)
+
+	    # Invert the matrix.
+	    if (version == 1)
+		call invers (Memr[DP_AC(allstar)], cdimen, nterm, flag)
+	    else
+		call invers2 (Memr[DP_AC(allstar)], cdimen, nterm, flag)
+
+	    if (dp_checkc (Memr[DP_AC(allstar)], cdimen, nterm,
+	        Memi[DP_APID(apsel)], Memr[DP_APMAG(apsel)],
+		Memi[DP_ASKIP(allstar)], Memi[DP_AIER(allstar)], istar, niter,
+		DP_RECENTER(dao), DP_VERBOSE(dao)))
+	        return (lstar)
+
+	    # Compute the coefficients.
+	    call mvmul (Memr[DP_AC(allstar)], cdimen, nterm,
+	        Memr[DP_AV(allstar)], Memr[DP_AX(allstar)])
+
+	    # In the beginning, the brightness of each star will be permitted
+	    # to change by no more than 2  magnitudes per iteration, and the
+	    # x,y coordinates of each centroid will be permitted to change by
+	    # no more than 0.4 pixel per iteration. Any time that the
+	    # parameter correction changes sign from one iteration to the next
+	    # the maximum permissable change will be reduced by a factor of
+	    # two. These clamps are released any time a star in the group
+	    # disappears.
+
+	    if (dp_alclamp (dao, im, Memi[DP_APID(apsel)],
+	        Memr[DP_APXCEN(apsel)], Memr[DP_APYCEN(apsel)],
+		Memr[DP_APMAG(apsel)], Memr[DP_APERR(apsel)],
+		Memr[DP_ASUMWT(allstar)], Memi[DP_ASKIP(allstar)],
+		Memr[DP_AXOLD(allstar)], Memr[DP_AXCLAMP(allstar)],
+		Memr[DP_AYOLD(allstar)], Memr[DP_AYCLAMP(allstar)], istar,
+		lstar, Memr[DP_AC(allstar)], Memr[DP_AX(allstar)], cdimen,
+		niter, clampmax, pererr, peakerr)) {
+
+	        # Flush the new data to the output buffers. Actually
+		# only data which fits this criterion need be written.
+		# However this makes sequential i/o management tricky.
+		# Leave this alone for the moment. Note only the cache-
+		# state is affected.
+	    }
+	
+	    # If there is more than one star, check to see whether any two
+	    # stars have merged.
+
+	    if (nstar > 1) {
+
+	        if (dp_almerge (Memr[DP_APXCEN(apsel)], Memr[DP_APYCEN(apsel)],
+		    Memr[DP_APMAG(apsel)], Memr[DP_APERR(apsel)],
+		    Memi[DP_ASKIP(allstar)], istar, lstar, sepcrit, sepmin,
+		    wcrit, j, k)) {
+
+		    call dp_alcentroid (Memr[DP_APXCEN(apsel)],
+		        Memr[DP_APYCEN(apsel)], Memr[DP_APMAG(apsel)], j, k)
+
+		    # Remove the k-th star from this group.
+		    if (DP_VERBOSE(dao) == YES) {
+		        call printf (
+			"REJECTING: Star ID: %d has merged with star ID: %d\n")
+			    call pargi (Memi[DP_APID(apsel)+j-1])
+			    call pargi (Memi[DP_APID(apsel)+k-1])
+		    }
+		    Memr[DP_APMAG(apsel)+j-1] = INDEFR
+		    Memi[DP_ASKIP(allstar)+j-1] = YES
+		    Memi[DP_AIER(allstar)+j-1] = ALLERR_MERGE
+
+		    # Loosen the clamps of every star in the group
+		    call aclrr (Memr[DP_AXOLD(allstar)+istar-1], nstar)
+		    call aclrr (Memr[DP_AYOLD(allstar)+istar-1], nstar)
+		    call amaxkr (Memr[DP_AXCLAMP(allstar)+istar-1], 0.5 *
+		        clampmax, Memr[DP_AXCLAMP(allstar)+istar-1], nstar)
+		    call amaxkr (Memr[DP_AYCLAMP(allstar)+istar-1], 0.5 *
+			clampmax, Memr[DP_AYCLAMP(allstar)+istar-1], nstar)
+
+	        }
+	    }
+
+	    # If the number of iterations completed is less than or equal
+	    # to 3, perform another iteration no questions asked.
+
+	    if (niter <= (MIN_ITER - 1))
+		return (lstar)
+
+	    # If not star has been removed from the group in the previous
+	    # iteration, check to see if any of the stars is too faint (more
+	    # than 12.5 magnitudes fainter thanthe PSF star). If several stars
+	    # are too faint, delete the faintest one, and set the 
+	    # brightnesses of the other faint ones exactly to 12.5
+	    # magnitudes below the PSF star. That way on the next iteration
+	    # we will see whether these  stars want to grow or disappear.
+
+	    if (dp_alfaint (Memr[DP_APMAG(apsel)], Memi[DP_ASKIP(allstar)],
+	        istar, lstar, faint, ifaint))
+		return (lstar)
+
+	    # If at least one star is more than 12.5 magnitudes fainter
+	    # than the PSF, then IFAINT is the index of the faintest of them
+	    # and faint is the relative brightness of the faintest of them.
+
+	    if (ifaint > 0) {
+
+	        if (DP_VERBOSE(dao) == YES) {
+		    call printf ("REJECTING: Star ID: %d is too faint\n")
+		        call pargi (Memi[DP_APID(apsel)+ifaint-1])
+	        }
+	        Memr[DP_APMAG(apsel)+ifaint-1] = INDEFR
+	        Memi[DP_ASKIP(allstar)+ifaint-1] = YES
+	        Memi[DP_AIER(allstar)+ifaint-1] = ALLERR_FAINT
+
+	        call aclrr (Memr[DP_AXOLD(allstar)+istar-1], nstar)
+	        call aclrr (Memr[DP_AYOLD(allstar)+istar-1], nstar)
+	        call amaxkr (Memr[DP_AXCLAMP(allstar)+istar-1], 0.5 *
+	            clampmax, Memr[DP_AXCLAMP(allstar)+istar-1], nstar)
+	        call amaxkr (Memr[DP_AYCLAMP(allstar)+istar-1], 0.5 *
+		    clampmax, Memr[DP_AYCLAMP(allstar)+istar-1], nstar)
+		
+	    # If no star is more than 12.5 magnitudes fainter than the PSF
+	    # then after the 50th iteration delete the least certain star i
+	    # if it is less than a one sigma detection; after the 100th
+	    # iteration delete the least certain star if it is less than a 
+	    # two sigma detection; after the 150th delete the least certain
+	    # star if it is less than a 3 sigma detection.
+
+	    } else if (niter >= 5) {
+
+	        call dp_almfaint (Memr[DP_APMAG(apsel)], Memr[DP_APERR(apsel)],
+		    istar, lstar, faint, ifaint)
+
+	        if (faint < wcrit) {
+
+	            if (DP_VERBOSE(dao) == YES) {
+		        call printf ("REJECTING: Star ID: %d is too faint\n")
+		            call pargi (Memi[DP_APID(apsel)+ifaint-1])
+	            }
+	            Memr[DP_APMAG(apsel)+ifaint-1] = INDEFR
+	            Memi[DP_ASKIP(allstar)+ifaint-1] = YES
+	            Memi[DP_AIER(allstar)+ifaint-1] = ALLERR_FAINT
+
+	            # Loosen the clamps of every star in the group.
+	            call aclrr (Memr[DP_AXOLD(allstar)+istar-1], nstar)
+	            call aclrr (Memr[DP_AYOLD(allstar)+istar-1], nstar)
+	            call amaxkr (Memr[DP_AXCLAMP(allstar)+istar-1], 0.5 *
+	                clampmax, Memr[DP_AXCLAMP(allstar)+istar-1], nstar)
+	            call amaxkr (Memr[DP_AYCLAMP(allstar)+istar-1], 0.5 *
+		        clampmax, Memr[DP_AYCLAMP(allstar)+istar-1], nstar)
+	        }
+	    }
+
+	    break
+	}
+
+	return (lstar)
+end
+
+
+# DP_LASTSTAR -- Find the last star in the current group.
+
+int procedure dp_laststar (last, istar, ntot)
+
+int	last[ARB]		# the grouping information array
+int	istar			# index of the first star in the group
+int	ntot			# total number of stars
+
+int	lstar
+
+begin
+	do lstar = istar, ntot, 1 {
+	    if (last[lstar] == YES)
+		break
+	}
+
+	return (lstar)
+end
+
+
+# DP_DELFAINTEST -- Delete the faintest star in the group.
+
+int procedure dp_delfaintest (mag, istar, lstar)
+
+real	mag[ARB]		# the array of magnitudes
+int	istar			# the first star in the group
+int	lstar			# the last star in the group
+
+int	i, starno
+real	faint
+
+begin
+	starno = 0
+	faint = MAX_REAL
+	do i = istar, lstar {
+	    if (mag[i] >= faint)
+		next
+	    faint = mag[i]
+	    starno = i
+	}
+
+	return (starno)
+end
+
+
+# DP_ALMSKY -- Determine the mean sky value for the current group of stars.
+
+real procedure dp_almsky (sky, nstar)
+
+real	sky[ARB]		# array of sky values
+int	nstar			# number of stars in the group
+
+int	i, nsky
+real 	sky_sum
+
+begin
+	sky_sum = 0.0
+	nsky = 0
+
+	do i = 1, nstar {
+	    if (IS_INDEFR(sky[i]))
+		next
+	    sky_sum = sky_sum + sky[i]
+	    nsky = nsky + 1
+	}
+
+	if (nsky <= 0)
+	    return (INDEFR)
+	else
+	    return (sky_sum / nsky)
+end
+
+
+# DP_ALACCUM -- Accumulate the data into the matrix.
+
+procedure dp_alaccum (dao, im, data, dnx, dny, dxoff, dyoff, subt, snx, sny,
+	sxoff, syoff, weights, wnx, wny, wxoff, wyoff, nxpix, nypix, ixmin,
+	iymin, mean_sky, istar, lstar, niter, clip, pererr, peakerr, sumres,
+	grpwt, chigrp, cdimen, nterm)
+
+pointer	dao			# pointer to the daophot structure
+pointer	im			# the input image descriptor
+real	data[dnx,dny]		# the subtracted data array
+int	dnx, dny		# dimenions of the data array
+int	dxoff, dyoff		# lower left corner of the data array
+real	subt[snx,sny]		# the scratch array
+int	snx, sny		# dimensions of the scratch array
+int	sxoff, syoff		# lower left corner of the scratch array
+real	weights[wnx,wny]	# the weight array
+int	wnx, wny		# dimensions of the weight array
+int	wxoff, wyoff		# lower left corner of the weight array
+int	nxpix, nypix		# the dimensions of the area of interest
+int	ixmin,iymin		# lower left corner of area of interest
+real	mean_sky		# the group sky value
+int	istar			# the index of the first star in current group
+int	lstar			# the index of the last star in current group
+int	niter			# the current interation
+bool	clip			# clip the errors ?
+real	pererr			# flat fielding error factor
+real	peakerr			# profile error factor
+real	sumres			# sum of the residuals
+real	grpwt			# the group weight
+real	chigrp			# the group chi value
+int	cdimen			# maximum number of maxtrix dimensions
+int	nterm			# number of terms in the matrix to fit
+
+real	fitradsq, maxgdata, fix, fiy, d, delta, sigmasq, relerr, wt, dwt
+pointer	psffit, apsel, allstar
+int	dix, diy, sxdiff, sydiff, wxdiff, wydiff
+real	sky_value, dp_alskyval()
+bool	dp_alomit()
+
+begin
+	# Set up some pointers.
+	apsel = DP_APSEL(dao)
+	psffit = DP_PSFFIT(dao)
+	allstar = DP_ALLSTAR(dao)
+
+	# These constants need to be stored more permanently in the
+	# allstar structure at some point. They should all be defined
+	# once and for all at task startup. Leave for next phase
+	# of code cleanup.
+
+	fitradsq = DP_FITRAD(dao) ** 2
+	if (IS_INDEFR(DP_MAXGDATA(dao)))
+	    maxgdata = MAX_REAL
+	else
+	    maxgdata = DP_MAXGDATA(dao)
+
+	# Compute the array offsets.
+	sxdiff = dxoff - sxoff
+	sydiff = dyoff - syoff 
+	wxdiff = dxoff - wxoff
+	wydiff = dyoff - wyoff
+
+	do diy = iymin - dyoff + 1, iymin - dyoff + nypix, 1 {
+	    fiy = real (diy + dyoff - 1)
+	    do dix = ixmin - dxoff + 1, ixmin - dxoff + nxpix, 1 {
+
+		# Skip data with negative weights.
+		if (weights[dix+wxdiff,diy+wydiff] < 0.0) 
+		    next
+		fix = real (dix + dxoff - 1)
+
+		# If this current pixel is within one fitting radius of
+		# at least one star in the current group, include it in the
+		# calculation. Otherwise skip it. While figuring this out,
+		# compute the squared distance of this pixel from the
+		# centroid of each star in the group.
+
+		if (dp_alomit (Memr[DP_APXCEN(apsel)], Memr[DP_APYCEN(apsel)],
+		    Memr[DP_ARPIXSQ(allstar)], istar, lstar, fix, fiy,
+		    fitradsq))
+		    next
+
+		call dp_alsetskip (Memr[DP_ARPIXSQ(allstar)],
+		    Memi[DP_ASKIP(allstar)], istar, lstar, fitradsq)
+
+		if (DP_GROUPSKY(dao) == NO) {
+		    sky_value = dp_alskyval (Memr[DP_APMSKY(apsel)],
+		        Memi[DP_ASKIP(allstar)], istar, lstar)
+		    if (IS_INDEFR(sky_value))
+			sky_value = mean_sky
+		} else
+		    sky_value = mean_sky
+
+		# The expected random error in the pixel is the quadratic
+		# sum of the Poisson statistics, plus the readout noise, 
+		# plus an estimated error of 0.75% of the total brightness
+		# for the difficulty of flat-fielding and bias-correcting 
+		# the chip, plus an estimated error of some fraction of the
+		# of the fourth derivative at the peak of the profile, to
+		# account for the difficulty of accurately interpolating
+		# within the PSF. The fourth derivative of the PSF is
+		# proportional to H/hwhm ** 4 (hwhm is the width of the
+		# stellar core); using the geometric mean of hwhmx and
+		# hwhmy, this becomes H/(hwhmx * hwhmy) ** 2. The ratio of
+		# the fitting error to this quantity is estimated from a
+		# good-seeing CTIO frame to be approimxately 0.027.
+
+		#d = subt[dix+sxdiff,diy+sydiff] - mean_sky
+		d = subt[dix+sxdiff,diy+sydiff] - sky_value
+		delta = max (0.0, data[dix,diy] - d)
+		if ((delta > maxgdata) && (niter >= MIN_ITER))
+		    next
+
+		# Dpos = raw data - residual
+		#      = model-predicted brightness in this pixel consisting
+		#        of sky plus all stellar profiles, which is
+		#        presumably non-negative
+		#
+		# The four noise sources in the model are
+		#        readout noise
+		#        poisson noise
+		#        flat-field errors
+		#        profile errors 
+		#
+		# Numerically the squares of these quantities are
+		#        ronoise = sigma[i,j]
+		#        poisson noise = delta / gain
+		#        flat-field error = constant * delta ** 2
+		#        profile error = constant * sum of profile ** 2
+
+		#sigmasq = weights[dix+wxdiff,diy+wydiff] + delta / 
+		    #DP_PHOTADU(dao) + (pererr * delta) ** 2 +
+		    #(peakerr * (delta - mean_sky)) ** 2
+		sigmasq = weights[dix+wxdiff,diy+wydiff] + delta / 
+		    DP_PHOTADU(dao) + (pererr * delta) ** 2 +
+		    (peakerr * (delta - sky_value)) ** 2
+		if (sigmasq <= 0.0)
+		    next
+		relerr = abs (d) / sqrt (sigmasq)
+
+		if (clip && (relerr > MAX_RELERR))
+		    next
+
+		# Now include this pixel in the fitting equation for the
+		# group.
+
+		wt = 0.0
+		call dp_alxaccum (dao, im,  Memr[DP_APXCEN(apsel)],
+		    Memr[DP_APYCEN(apsel)], Memr[DP_APMAG(apsel)],
+		    Memr[DP_ARPIXSQ(allstar)], Memi[DP_ASKIP(allstar)],
+		    Memr[DP_AX(allstar)], Memr[DP_ANUMER(allstar)],
+		    Memr[DP_ADENOM(allstar)], istar, lstar, fix, fiy, d,
+		    wt, sigmasq, nterm, fitradsq)
+
+		# At this point the vector X contains  the first
+		# derivative of the condition equation for pixel (I,J)
+		# with respect to each of the fitting parameters for
+		# all of the stars. Now these derivatives will be added 
+		# into the normal matrix and the vector of residuals.
+		# Add this residual into the weighted sum of the absolute
+		# relative residuals
+
+		dwt = wt * relerr
+		sumres = sumres + dwt
+		grpwt = grpwt + wt
+
+		# SUMRES is the weighted sum for all the pixels in the group.
+		# Now also add the weigghted value of the residuals into the
+		# accumulating sum for each of the stars.
+
+		call dp_alchi (Memi[DP_ASKIP(allstar)], Memr[DP_APCHI(apsel)],
+		    Memr[DP_ASUMWT(allstar)], Memi[DP_ANPIX(allstar)], istar,
+		    lstar, wt, dwt)
+
+		# Up until now, WT represents only the radial weighting 
+		# profile. Now figure in the anticipated standard error
+		# of the pixel.
+
+		wt = wt / sigmasq
+
+		# After the third iteration, reduce the weight of a bad pixel
+		# Note that for the first iteration, only the stellar magnitude
+		# is being solved for, which is a problem in LINEAR least 
+		# squares, and so should be solved exactly the first time.
+		# After that, the star is given two iterations to adjust it's
+		# centroid before the clipping is turned on.  After that a
+		# pixel having a residual of DP_CLIPRANGE times sigma gets
+		# reduced to half weight;  a pixel having a residual of n
+		# sigma gets a weight of 1 / [1 + (n/DP_CLIPRANGE) **
+		# DP_CLIPEXP].
+
+		if (clip)
+		    wt = wt / (1.0 + (relerr / (chigrp * DP_CLIPRANGE (dao))) **
+			 DP_CLIPEXP (dao))
+
+		# Now work this pixel into the normal matrix
+		dwt = d * wt
+		call dp_alcaccum (Memr[DP_AX(allstar)], Memr[DP_AC(allstar)],
+		    Memr[DP_AV(allstar)], Memi[DP_ASKIP(allstar)], cdimen,
+		    nterm, istar, lstar, wt, dwt)
+	    }
+	}
+end
+
+
+# DP_ALOMIT - Omit pixels which are too far from the center of any star
+# in the group.
+
+bool procedure dp_alomit (xcen, ycen, rpixsq, istar, lstar, fix, fiy,
+	fitradsq)
+
+real	xcen[ARB]		# array of x centers
+real	ycen[ARB]		# array of y centers
+real	rpixsq[ARB]		# array of pixel differences
+int	istar, lstar		# first and last star
+real	fix, fiy		# current star position
+real	fitradsq		# fit radius squared
+
+bool	omit
+int	i
+
+begin
+	omit = true
+	do i = istar, lstar {
+	    rpixsq[i] = (fix - xcen[i]) ** 2 + (fiy - ycen[i]) ** 2
+	    if (rpixsq[i] < fitradsq)
+		omit = false
+	}
+
+	return (omit)
+end
+
+
+# DP_ALSETSKIP -- Initialize the skip array.
+
+procedure dp_alsetskip (rpixsq, skip, istar, lstar, fitradsq)
+
+real	rpixsq[ARB]		# pixel distance squared
+int	skip[ARB]		# skip array
+int	istar			# the index of the first star
+int	lstar			# the index of the last star
+real	fitradsq		# the fitting radius squared
+
+int	i
+
+begin
+	do i = istar, lstar {
+	    if (rpixsq[i] < fitradsq)
+		skip[i] = NO
+	    else
+		skip[i] = YES
+	}
+end
+
+
+# DP_ALSKYVAL -- Compute the sky value to be subtracted from a given pixel
+# by averaging the sky values of all the the stars for which the pixel in
+# question is within one fitting radius.
+
+real procedure dp_alskyval (sky, skip, istar, lstar)
+
+real	sky[ARB]		# the array of sky values
+int	skip[ARB]		# the array of skip values
+int	istar			# the index of the first star
+int	lstar			# the index of the last star
+
+int	i, npts
+real	sum
+
+begin
+	sum = 0.0
+	npts = 0
+	do i = istar, lstar {
+	    if (skip[i] == YES)
+		next
+	    sum = sum + sky[i]
+	    npts = npts + 1
+	}
+
+	if (npts <= 0)
+	    return (INDEFR)
+	else
+	    return (sum / npts)
+end
+
+
+# DP_ALXACCUM - Accumulate the x vector.
+
+procedure dp_alxaccum (dao, im, xcen, ycen, mag, rpixsq, skip, x, numer1,
+	denom1, istar, lstar, fix, fiy, resid, wt, sigmasq, nterm, fitradsq)
+
+pointer	dao			# pointer to the daopot structure
+pointer	im			# pointer to the input image
+real	xcen[ARB]		# array of x centers
+real	ycen[ARB]		# array of y centers
+real	mag[ARB]		# array of relative brightnesses
+real	rpixsq[ARB]		# array of pixel distances squared
+int	skip[ARB]		# the include / exclude array
+real	x[ARB]			# the x vector to be accumulated
+real	numer1[ARB]		# first numerator array
+real	denom1[ARB]		# first denominator array
+int	istar			# index of the first star in group
+int	lstar			# index of the last star in group
+real	fix			# the x position of the current pixel
+real	fiy			# the y position of the current pixel
+real	resid			# the input data residual
+real	wt			# the output weight value
+real	sigmasq			# the sigma squared value
+int	nterm			# number of terms in the matrix
+real	fitradsq		# fitting radius squared
+
+real	psfsigsqx, psfsigsqy, dx, dy, deltax, deltay, value, radsq, rhosq
+real	dvdx, dvdy, dfdsig
+pointer	psffit
+int	nstar, i, i3, k
+real	dp_usepsf()
+
+begin
+	psffit = DP_PSFFIT(dao)
+
+	nstar = lstar - istar + 1
+	psfsigsqx = Memr[DP_PSFPARS(psffit)]
+	psfsigsqy = Memr[DP_PSFPARS(psffit)+1]
+
+	do i = istar, lstar {
+
+	    if (skip[i] == YES)
+		next
+	    radsq = rpixsq[i] / fitradsq
+	    if (radsq >= MAX_RSQ) 
+		next
+	    wt = max (wt, 5.0 / (5.0 + radsq / (1.0 - radsq)))
+		    
+	    # The condition equation for pixel[I,J] is the following.
+	    # 
+	    # data[I,J] - sum [scale * psf[I,J]] - mean_sky = residual
+	    #
+	    # Then we will adjust the scale factors, x centers and ycenters
+	    # such that sum [weight * residual ** 2] is minimized.
+	    # WT will be a function (1) of the distance of this pixel from
+	    # the center of the nearest star, (2) of the model predicted
+	    # brightness of the pixel (taking into consideration the
+	    # readout noise, the photons/ADU, and the interpolation error
+	    # of the PSF), and (3) of the size of the residual itself. (1) is
+	    # necessary to prevent the non-linear least-squares from
+	    # fit from oscillating. (2) is simply sensible weighting and (3)
+	    # is a crude attempt at making the solution more robust
+	    # against bad pixels.
+
+	    dx = fix - xcen[i]
+	    dy = fiy - ycen[i]
+	    call dp_wpsf (dao, im, xcen[i], ycen[i], deltax, deltay, 1)
+	    deltax = (deltax - 1.0) / DP_PSFX(psffit) - 1.0
+	    deltay = (deltay - 1.0) / DP_PSFY(psffit) - 1.0
+	    value = 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), deltax, deltay,
+	        dvdx, dvdy)
+
+	    if (nterm > nstar) {
+		i3 = (i - istar + 1) * 3
+		k = i3 - 2
+		x[k] = -value
+		k = i3 - 1
+		x[k] = -mag[i] * dvdx
+		x[i3] = -mag[i] * dvdy
+	    } else {
+		k = i - istar + 1
+		x[k] = -value
+	    }
+
+	    rhosq =  (dx / psfsigsqx) ** 2 + (dy / psfsigsqy) ** 2
+	    if (rhosq > MAX_RHOSQ)
+		next
+	    rhosq = 0.6931472 * rhosq
+	    dfdsig = exp (-rhosq) * (rhosq - 1.)
+	    numer1[i] = numer1[i] + dfdsig * resid / sigmasq
+	    denom1[i] = denom1[i] + dfdsig ** 2 / sigmasq
+	}
+end
+
+
+# DP_ALCHI -- Compute the weights and chis.
+
+procedure dp_alchi (skip, chi, sumwt, npix, istar, lstar, wt, dwt)
+
+int	skip[ARB]		# include / exclude array for the stars
+real	chi[ARB]		# array of chi values
+real	sumwt[ARB]		# array of weight sums
+int	npix[ARB]		# array of pixel counts
+int	istar			# index of first star
+int	lstar			# index of last star
+real	wt			# input weight
+real	dwt			# weighted residual
+
+int	i
+
+begin
+	do i = istar, lstar {
+	    if (skip[i] == YES)
+		next
+	    chi[i] = chi[i] + dwt
+	    sumwt[i] = sumwt[i] + wt
+	    npix[i] = npix[i] + 1
+	}
+end
+
+
+# DP_ALCACCUM -- Accumulate the main matrix.
+
+procedure dp_alcaccum (x, c, v, skip, cdimen, nterm, istar, lstar, wt, dwt)
+
+
+real	x[ARB]			# x array
+real	c[cdimen,ARB]		# maxtrix array
+real	v[ARB]			# v vector
+int	skip[ARB]		# include / exclude array
+int	cdimen			# maximum size of c matrix
+int	nterm			# number of terms to be fit
+int	istar			# index of the first star
+int	lstar			# index of the last star 
+real	wt			# input weight
+real	dwt			# input residual weight
+
+int	nstar, i, j, k, l, i3, i3m2, j3
+real	xkwt
+
+begin
+	nstar = lstar - istar + 1
+	do i = istar, lstar {
+	    if (skip[i] == YES)
+		next
+	    if (nterm > nstar) {
+		i3 = (i - istar + 1) * 3
+		i3m2 = i3 - 2
+		do k = i3m2, i3 
+		    v[k] = v[k] + x[k] * dwt
+		do j = istar, i {
+		    if (skip[j] == YES)
+			next
+		    j3 = (j - istar + 1) * 3
+		    do k = i3m2, i3 {
+			xkwt = x[k] * wt
+			do l = j3 - 2, min (k, j3)
+			    c[k,l] = c[k,l] + x[l] * xkwt
+		    }
+		}
+	    } else {
+		k = i - istar + 1
+		v[k] = v[k] + x[k] * dwt
+		xkwt = x[k] * wt
+		do j = istar, i {
+		    if (skip[j] == YES)
+			next
+		    l = j - istar + 1
+		    c[k,l] = c[k,l] + x[l] * xkwt
+	        }
+	    }
+	}
+end
+
+
+# DP_ALREDO -- Check to see that there are enough good pixels to fit the
+# star and compute the individual chi values.
+
+bool procedure dp_alredo (id, mag, chi, sumwt, npix, skip, aier, istar, lstar,
+	niter, chigrp, center, verbose)
+
+int	id[ARB]			# the array of star ids
+real	mag[ARB]		# the array of relative brightnesses
+real	chi[ARB]		# the array of chi values
+real	sumwt[ARB]		# the array of weight values
+int	npix[ARB]		# the array of pixel counts
+int	skip[ARB]		# the include / exclude array
+int	aier[ARB]		# the array of error codes
+int	istar			# index of the first star
+int	lstar			# index of the last star
+int	niter			# the current iteration
+real	chigrp			# chi value of the group
+int	center			# recenter the values ?
+int	verbose			# verbose mode ?
+
+bool	redo
+int	i
+
+begin
+	redo = false
+	do i = istar, lstar {
+	    if (center == YES && (niter >= 2)) {
+		if (npix[i] < MIN_NPIX) {
+		    redo = true
+		    skip[i] = YES
+		    if (verbose == YES) {
+			call printf (
+			    "REJECTING: Star ID: %d has too few good pixels\n")
+			call pargi (id[i])
+		    }
+		    aier[i] = ALLERR_NOPIX
+		    mag[i] = INDEFR
+		} else {
+		    skip[i] = NO
+		    if (sumwt[i] > MIN_SUMWT) {
+			chi[i] = CHI_NORM * chi[i] / sqrt (sumwt[i] *
+			    (sumwt[i] - MIN_SUMWT))
+			sumwt[i] = ((sumwt[i] - MIN_SUMWT) * chi[i] +
+			    MIN_SUMWT) / sumwt[i]
+		    } else
+			chi[i] = chigrp
+		}
+	    } else {
+		if (npix[i] < MIN_NPIX) {
+		    redo = true
+		    skip[i] = YES
+		    if (verbose == YES) {
+			call printf (
+			    "REJECTING: Star ID: %d has too few good pixels\n")
+			call pargi (id[i])
+		    }
+		    aier[i] = ALLERR_NOPIX
+		    mag[i] = INDEFR
+		} else {
+		    skip[i] = NO
+		    if (sumwt[i] > 1.0) {
+			chi[i] = CHI_NORM * chi[i] / sqrt (sumwt[i] *
+			    (sumwt[i] - 1.0))
+			sumwt[i] = ((sumwt[i] - MIN_SUMWT) * chi[i] +
+			    MIN_SUMWT) / sumwt[i]
+		    } else
+			chi[i] = chigrp
+		}
+	    }
+	}
+
+	return (redo)
+end
+
+
+# DP_CHECKC - Check the c matrix for valid diagonal elements.
+
+bool procedure  dp_checkc (c, cdimen, nterm, id, mag, skip, aier, istar, niter,
+	recenter, verbose)
+
+real	c[cdimen,ARB]		# the c matrix
+int	cdimen			# maximum size of the c matrix
+int	nterm			# number of terms to be fit
+int	id[ARB]			# the array of ids
+real	mag[ARB]		# the array of relative brightnesses
+int	skip[ARB]		# the include / exclude array
+int	aier[ARB]		# the array of error codes
+int	istar			# index of the first star
+int	niter			# the iteration number
+int	recenter		# recenter ?
+int	verbose			# verbose mode ?
+
+bool	redo
+int	j, starno
+
+begin
+	redo = false
+	do j = 1, nterm {
+	    if (c[j,j] > 0.0)
+		next
+	    if ((recenter == YES) && (niter >= 2))
+		starno = (j + 2) / 3
+	    else
+		starno = j
+	    starno = starno + istar - 1
+	    skip[starno] = YES
+	    if (verbose == YES) {
+		call printf (
+		    "REJECTING: Star ID: %d generated a fitting error\n")
+		    call pargi (id[starno])
+	    }
+	    aier[starno] = ALLERR_SINGULAR
+	    mag[starno] = INDEFR
+	    redo = true
+	    break
+	}
+
+	return (redo)
+end
+
+
+# DP_ALCLAMP -- Get the answers.
+
+bool procedure dp_alclamp (dao, im, id, xcen, ycen, mag, magerr, sumwt,
+	skip, xold, xclamp, yold, yclamp, istar, lstar, c, x,
+	cdimen, niter, clampmax, pererr, peakerr)
+
+pointer	dao			# pointer to the daophot structure
+pointer	im			# the input image descriptor
+int	id[ARB]			# array of star ids
+real	xcen[ARB]		# array of star x centers
+real	ycen[ARB]		# array of star y centers
+real	mag[ARB]		# array of relative brightnesses
+real	magerr[ARB]		# array of relative brightness errors
+real	sumwt[ARB]		# array of pixel distance squared values
+int 	skip[ARB]		# include / exclude array
+real	xold[ARB]		# xold array
+real	xclamp[ARB]		# xclamp array
+real	yold[ARB]		# yold array
+real	yclamp[ARB]		# yclamp array
+int	istar			# the index of the first star
+int	lstar			# the index of the last star
+real	c[cdimen,ARB]		# c matrix
+real	x[ARB]			# x vector
+int	cdimen			# the maximum size of c matrix 
+int	niter			# the current iteration
+real	clampmax		# the maximum clamp value
+real	pererr			# flat field error
+real	peakerr			# the profile error
+
+bool	redo, bufwrite
+int	i, l, k, j, lx, mx, ly, my, nx, ny
+pointer	psffit, allstar
+real	dwt, psfrad, psfradsq, df
+
+begin
+	# Get some daophot pointers.
+	allstar = DP_ALLSTAR(dao)
+	psffit = DP_PSFFIT(dao)
+
+	# Get some constants.
+	if (DP_PSFSIZE(psffit) == 0)
+	    psfrad = DP_PSFRAD(dao)
+	else
+	    psfrad = (real (DP_PSFSIZE(psffit) - 1) / 2.0 - 1.0) / 2.0
+	psfradsq = psfrad * psfrad
+
+	# Initialize.
+	bufwrite = false
+
+	do i = istar, lstar {
+
+	    if ((DP_RECENTER(dao) == YES)  && (niter >= 2)) {
+
+		l = 3 * (i - istar + 1)
+		k = l - 1
+		j = l - 2
+
+		# Note that the sign of the correction is such that it must be
+		# SUBTRACTED from the current value of the parameter to get the
+		# improved parameter value. This being the case, if the 
+		# correction to the brightness is negative a change of one
+		# magnitude is a change of factor of 2.5; if the brightness
+		# correction is positive a change of one magnitude is a change
+		# of 60%.
+
+		dwt = xold[i] * x[k]
+		if (dwt < 0.0)
+		    xclamp[i] = max (MIN_XYCLAMP, MIN_XYCLAMP_FRACTION *
+		        xclamp[i])
+		else
+		    xclamp[i] = min (clampmax, MAX_XYCLAMP_FRACTION * xclamp[i])
+		xcen[i] = xcen[i] - x[k] / (1.0 + abs (x[k] / xclamp[i]))
+		xold[i] = x[k]
+		
+		dwt = yold[i] * x[l]
+		if (dwt < 0.0)
+		    yclamp[i] = max (MIN_XYCLAMP, MIN_XYCLAMP_FRACTION *
+		        yclamp[i])
+		else
+		    yclamp[i] = min (clampmax, MAX_XYCLAMP_FRACTION * yclamp[i])
+		ycen[i] = ycen[i] - x[l] / (1.0 + abs (x[l] / yclamp[i]))
+		yold[i] = x[l]
+
+		mag[i] = mag[i] - x[j] / (1.0 + max (x[j] / (MAX_DELTA_FAINTER *
+		    mag[i]), -x[j] / (MAX_DELTA_BRIGHTER * mag[i])))
+		magerr[i] = sumwt[i] * sqrt (c[j,j])
+
+		if (niter >= MIN_ITER) {
+		    redo = false
+		    if (abs(x[j]) > max (MAX_NEW_ERRMAG * magerr[i],
+		        MAX_NEW_RELBRIGHT2 * mag[i]))
+			redo = true
+		    else {
+			df = (MAX_NEW_ERRMAG * sumwt[i]) ** 2
+			if ((x[k] ** 2) > max (df * c[k,k], MAX_PIXERR2))
+			    redo = true
+			else if ((x[l] ** 2) > max (df * c[l,l], MAX_PIXERR2))
+			    redo = true
+		    }
+		} else
+		    redo = true
+
+	    } else {
+		j = i - istar + 1
+		mag[i] = mag[i] - x[j] / (1.0 + 1.2 * abs (x[j] / mag[i]))
+		magerr[i] = sumwt[i] * sqrt (c[j,j])
+		if (niter >= 2) {
+		    redo = false
+		    if (abs(x[j]) > max (MAX_NEW_ERRMAG * magerr[i],
+			MAX_NEW_RELBRIGHT2 * mag[i]))
+			redo = true
+		} else
+		    redo = true
+	    }
+
+	    if (mag[i] < 2.0 * magerr[i])
+		redo = true
+	    if (niter >= DP_MAXITER(dao))
+		redo = false
+	    if (redo && (niter < DP_MAXITER(dao)))
+		next
+
+	    # If this star converged, write out the results for it,
+	    # flag it for deletion from the star list and subtract it
+	    # from the reference copy of the image.
+
+	    call dp_glim (xcen[i], ycen[i], psfrad, DP_DLX(allstar),
+	        DP_DMX(allstar), DP_DLY(allstar), DP_DMY(allstar),
+	        lx, mx, ly, my)
+	    nx = mx - lx + 1
+	    ny = my - ly + 1
+
+	    call dp_swstar (dao, im, Memr[DP_DBUF(allstar)], DP_DNX(allstar),
+		DP_DNY(allstar), DP_DXOFF(allstar), DP_DYOFF(allstar),
+		Memr[DP_WBUF(allstar)], DP_WNX(allstar), DP_WNY(allstar),
+		DP_WXOFF(allstar), DP_WYOFF(allstar), xcen[i], ycen[i],
+		mag[i], lx, ly, nx, ny, psfradsq, DP_PHOTADU(dao), pererr,
+		peakerr)
+
+	    skip[i] = YES
+
+	    if (DP_VERBOSE(dao) == YES) {
+		call dp_wout (dao, im, xcen[i], ycen[i], xcen[i], ycen[i], 1)
+		call printf (
+		"FITTING:   ID: %5d  XCEN: %8.2f  YCEN: %8.2f  MAG: %8.2f\n")
+		    call pargi (id[i])
+		    call pargr (xcen[i])
+		    call pargr (ycen[i])
+		if (mag[i] <= 0.0)
+		    call pargr (INDEFR)
+		else
+		    call pargr (DP_PSFMAG(psffit) - 2.5 * log10 (mag[i]))
+
+	    }
+	    bufwrite = true
+	}
+
+	return (bufwrite)
+end
+
+
+# DP_SWSTAR -- Subtract the fitted star for the data and adjust the
+# weights.
+
+procedure dp_swstar (dao, im, data, dnx, dny, dxoff, dyoff, weights, wnx, wny,
+	wxoff, wyoff, xcen, ycen, mag, lx, ly, nx, ny, psfradsq, gain,
+	pererr, peakerr)
+
+pointer	dao				# pointer to the daophot structure
+pointer	im				# pointer to the input image
+real	data[dnx,dny]			# the data array
+int	dnx, dny			# dimensions of the data array
+int	dxoff, dyoff			# lower left corner of data array
+real	weights[wnx,wny]		# the weights array
+int	wnx, wny			# dimensions of the weights array
+int	wxoff, wyoff			# lower left corner of weights array
+real	xcen, ycen			# the position of the star
+real	mag				# relative brightness of the star
+int	lx, ly				# lower left corner region of interest
+int	nx, ny				# size of region of interest
+real	psfradsq			# the psf radius squared
+real	gain				# the gain in photons per adu
+real	pererr				# the flat field error factor
+real	peakerr				# the peak error factor
+
+real	deltax, deltay, dx, dy, dysq, diff, dvdx, dvdy
+pointer	psffit
+int	di, dj, wxdiff, wydiff
+real	dp_usepsf()
+
+begin
+	psffit = DP_PSFFIT(dao)
+
+	wxdiff = dxoff - wxoff
+	wydiff = dyoff - wyoff
+	call dp_wpsf (dao, im, xcen, ycen, deltax, deltay, 1)
+	deltax = (deltax - 1.0) / DP_PSFX(psffit) - 1.0
+	deltay = (deltay - 1.0) / DP_PSFY(psffit) - 1.0
+
+	do dj = ly - dyoff + 1, ly - dyoff + ny {
+	    dy = (dj + dyoff - 1) - ycen
+	    dysq = dy * dy
+	    do di = lx - dxoff + 1, lx - dxoff + nx {
+		if (weights[di+wxdiff,dj+wydiff] <= -MAX_REAL)
+		    next
+	        dx = (di + dxoff - 1) - xcen
+		if ((dx * dx + dysq) >= psfradsq)
+		    next
+		diff = mag * 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), deltax, deltay,
+		    dvdx, dvdy)
+		data[di,dj] = data[di,dj] - diff
+		if (diff > 0.0) {
+		    diff = diff / gain + pererr *
+		        2.0 * max (0.0, data[di,dj]) * diff +
+			(peakerr * diff) ** 2
+		    if (weights[di+wxdiff,dj+wydiff] >= 0.0)
+			weights[di+wxdiff,dj+wydiff] = weights[di+wxdiff,
+			    dj+wydiff] + diff
+		    else
+			weights[di+wxdiff,dj+wydiff] = - (abs (weights[di+
+			    wxdiff,dj+wydiff]) + diff)
+		}
+	    }
+	}
+end
+
+
+# DP_ALMERGE -- Determine whether or not two stars should merge.
+
+bool procedure dp_almerge (xcen, ycen, mag, magerr, skip, istar, lstar, 
+	sepcrit, sepmin, wcrit, k, l)
+
+real	xcen[ARB]		# array of x centers
+real	ycen[ARB]		# array of y centers
+real	mag[ARB]		# array of relative brightnesses
+real	magerr[ARB]		# array of relative brightness errors
+int	skip[ARB]		# the include / exclude array
+int	istar			# index to the first star
+int	lstar			# index to the last star
+real	sepcrit			# the critical separation
+real	sepmin			# the minimum separation
+real	wcrit			# the critical error
+int	k, l			# indices of stars to be merged
+
+bool	merge
+int	i, j
+real	rsq, sep
+
+begin
+	# Initialize.
+	rsq = sepcrit
+	merge = false
+	k = 0
+	l = 0
+
+	# Find the closest two stars within the critical separation.
+	do i = istar + 1, lstar {
+
+	    if (skip[i] == YES)
+		next
+
+	    do j = istar, i - 1 {
+		if (skip[j] == YES)
+		    next
+		sep = (xcen[i] - xcen[j]) ** 2 + (ycen[i] - ycen[j]) ** 2
+		if (sep >= rsq)
+		    next
+		    
+		# Two stars are overlapping. Identify the fainter of the two.
+		rsq = sep
+		if (mag[i] < mag[j]) {
+		    k = i
+		    l = j
+		} else {
+		    k = j
+		    l = i
+		}
+
+		merge = true
+	    }
+	}
+
+	# No stars overlap.
+	if (! merge)
+	    return (merge)
+
+	# The stars are not close enough.
+	if ((rsq > sepmin) && (mag[k] > wcrit * magerr[k])) 
+	    merge = false
+
+	return (merge)
+end
+
+
+# DP_ALCENTROID -- Merge two stars by adding their relative brightnesses
+# and replacing their individual positions with the relative brightness
+# weighted mean position.
+
+procedure dp_alcentroid (xcen, ycen, mag, k, l)
+
+real	xcen[ARB]		# array of x centers
+real	ycen[ARB]		# array of y centers
+real	mag[ARB]		# array of magnitudes
+int	k, l			# input indices, k is fainter
+
+begin
+	xcen[l] = xcen[l] * mag[l] + xcen[k] * mag[k]
+	ycen[l] = ycen[l] * mag[l] + ycen[k] * mag[k]
+	mag[l] = mag[l] + mag[k]
+	xcen[l] = xcen[l] / mag[l]
+	ycen[l] = ycen[l] / mag[l]
+end
+
+
+# DP_ALFAINT -- Find all the stars with realtive brightnesss less than
+# a given number and set their magnitudes to a given minimum relative
+# brightness.
+
+bool procedure dp_alfaint (mag, skip, istar, lstar, faint, ifaint)
+
+real	mag[ARB]		# array of magnitudes
+int	skip[ARB]		# skip array
+int	istar, lstar		# first and last stars
+real	faint			# faintest star
+int	ifaint			# index of faintest star
+
+int	i
+
+begin
+	# Initialize
+	faint = MIN_REL_BRIGHT
+	ifaint = 0
+
+	do i = istar, lstar {
+	    if (skip[i] == YES)
+		return (true)
+	    if (mag[i] < faint) {
+		faint = mag[i] 
+		ifaint = i
+	    }
+	    if (mag[i] < MIN_REL_BRIGHT)
+	        mag[i] = MIN_REL_BRIGHT
+	}
+
+	return (false)
+end
+
+
+# DP_ALMFAINT -- Find the star with the greatest error in its relative
+# brightness.
+
+procedure dp_almfaint (mag, magerr, istar, lstar, faint, ifaint)
+
+real	mag[ARB]		# array of relative brightnesses
+real	magerr[ARB]		# array of relative brightness errors
+int	istar			# index of first star
+int	lstar			# index of last star
+real	faint			# computed faintness limit
+int	ifaint			# index of faintest star
+
+int	i
+real	wt
+
+begin
+	faint = MAX_REAL
+	ifaint = 0
+	do i = istar, lstar {
+	    wt = mag[i] / magerr[i]
+	    if (wt < faint) {
+		faint = wt
+		ifaint = i
+	    }
+	}
+end
diff --git a/noao/digiphot/daophot/allstar/dpalwrite.x b/noao/digiphot/daophot/allstar/dpalwrite.x
new file mode 100644
index 00000000..c4688924
--- /dev/null
+++ b/noao/digiphot/daophot/allstar/dpalwrite.x
@@ -0,0 +1,556 @@
+include	<tbset.h>
+include <time.h>
+include "../lib/daophotdef.h"
+include "../lib/apseldef.h"
+include "../lib/allstardef.h"
+
+# DP_TNEWAL -- Create an new ALLSTAR output ST table.
+
+procedure dp_tnewal (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 table definition.
+	call smark (sp)
+	call salloc (colnames, ALL_NOUTCOLUMN * (SZ_COLNAME + 1), TY_CHAR)
+	call salloc (colunits, ALL_NOUTCOLUMN * (SZ_COLUNITS + 1), TY_CHAR)
+	call salloc (colformat, ALL_NOUTCOLUMN * (SZ_COLFMT + 1), TY_CHAR)
+	call salloc (col_dtype, ALL_NOUTCOLUMN, TY_INT)
+	call salloc (col_len, ALL_NOUTCOLUMN, 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 ("%6d", 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 ("ADC", 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 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, ALL_NOUTCOLUMN 
+	    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], ALL_NOUTCOLUMN)
+	call tbtcre (tp)
+
+	# Write out the header parameters.
+	call dp_talpars (dao, tp)
+
+	call sfree (sp)
+	
+end
+
+
+define AL_NAME1STR "#N%4tID%10tXCENTER%20tYCENTER%30tMAG%42tMERR%56tMSKY%71t\
+NITER%80t\\\n"
+define AL_UNIT1STR "#U%4t##%10tpixels%20tpixels%30tmagnitudes%42tmagnitudes%56t\
+counts%71t##%80t\\\n"
+define AL_FORMAT1STR "#F%4t%%-9d%10t%%-10.3f%20t%%-10.3f%30t%%-12.3f%42t\
+%%-14.3f%56t%%-15.7g%71t%%-6d%80t \n"
+
+define AL_NAME2STR "#N%12tSHARPNESS%24tCHI%36tPIER%42tPERROR%80t\\\n"
+define AL_UNIT2STR "#U%12t##%24t##%36t##%42tperrors%80t\\\n"
+define AL_FORMAT2STR "#F%12t%%-23.3f%24t%%-12.3f%36t%%-6d%42t%%-13s%80t \n"
+
+
+# DP_XNEWAL -- Initialize a new output ALLSTAR text file.
+
+procedure dp_xnewal (dao, tp)
+
+pointer	dao		# pointer to the daophot structure
+int	tp		# the output file descriptor
+
+begin
+	# Write out the header parameters.
+	call dp_xalpars (dao, tp)
+
+	# Write out the banner.
+	call fprintf (tp, "#\n")
+	call fprintf (tp, AL_NAME1STR)
+	call fprintf (tp, AL_UNIT1STR)
+	call fprintf (tp, AL_FORMAT1STR)
+	call fprintf (tp, "#\n")
+	call fprintf (tp, AL_NAME2STR)
+	call fprintf (tp, AL_UNIT2STR)
+	call fprintf (tp, AL_FORMAT2STR)
+	call fprintf (tp, "#\n")
+end
+
+
+# DP_TALPARS -- Add various parameters to the header of the ALLSTAR table.
+
+procedure dp_talpars (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 workin 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", "allstar")
+
+	# 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, "ALLSTARFILE", 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])
+	}
+	call dp_imroot (DP_OUTIMAGE(dao), Memc[outstr], SZ_LINE)
+	call tbhadt (tp, "SUBIMAGE", 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, "GRPSKY", itob (DP_GROUPSKY(dao)))
+	call tbhadb (tp, "FITSKY", itob (DP_FITSKY(dao)))
+	call tbhadr (tp, "PSFMAG", DP_PSFMAG(psffit))
+	call tbhadr (tp, "PSFRAD", DP_SPSFRAD(dao))
+	call tbhadr (tp, "FITRAD", DP_SFITRAD(dao))
+	call tbhadr (tp, "ANNULUS", DP_SANNULUS(dao))
+	call tbhadr (tp, "DANNULUS", DP_SDANNULUS(dao))
+	call tbhadi (tp, "MAXITER", DP_MAXITER(dao))
+	call tbhadi (tp, "MAXGROUP", DP_MAXGROUP(dao))
+	#call tbhadi (tp, "MAXNSTAR", DP_MAXNSTAR(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 tbhadr (tp, "MERGERAD", DP_SMERGERAD(dao))
+
+	call sfree(sp)
+end
+
+
+# DP_XALPARS -- Add various parameters to the header of the PEAK table
+
+procedure dp_xalpars (dao, tp)
+
+pointer	dao			# pointer to the daophot structure
+int	tp			# the output file descriptor
+
+pointer	sp, psffit, outstr, date, time, dummy
+bool	itob()
+int	envfind()
+
+begin
+	# Define some daophot pointers.
+	psffit = DP_PSFFIT(dao)
+
+	# Allocate workin 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)
+	call salloc (dummy, 2, TY_CHAR)
+	Memc[dummy] = EOS
+
+	# 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", Memc[dummy])
+	if (envfind ("userid", Memc[outstr], SZ_LINE) <= 0)
+	    call dp_sparam (tp, "USER", Memc[outstr], "name", Memc[dummy])
+	call gethost (Memc[outstr], SZ_LINE)
+	call dp_sparam (tp, "HOST", Memc[outstr], "computer", Memc[dummy])
+	call dp_date (Memc[date], Memc[time], SZ_DATE)
+	call dp_sparam (tp, "DATE", Memc[date], "yyyy-mm-dd", Memc[dummy])
+	call dp_sparam (tp, "TIME", Memc[time], "hh:mm:ss", Memc[dummy])
+	call dp_sparam (tp, "PACKAGE", "daophot", "name", Memc[dummy])
+	call dp_sparam (tp, "TASK", "allstar", "name", Memc[dummy])
+
+	# Write out the files names.
+	call dp_imroot (DP_INIMAGE(dao), Memc[outstr], SZ_LINE)
+	call dp_sparam (tp, "IMAGE", Memc[outstr], "imagename", Memc[dummy])
+	call dp_froot (DP_INPHOTFILE(dao), Memc[outstr], SZ_LINE)
+	call dp_sparam (tp, "PHOTFILE", Memc[outstr], "filename", Memc[dummy])
+	call dp_imroot (DP_PSFIMAGE(dao), Memc[outstr], SZ_LINE)
+	call dp_sparam (tp, "PSFIMAGE", Memc[outstr], "imagename", Memc[dummy])
+	call dp_froot (DP_OUTPHOTFILE(dao), Memc[outstr], SZ_LINE)
+	call dp_sparam (tp, "ALLSTARFILE", Memc[outstr], "filename",
+	    Memc[dummy])
+	if (DP_OUTREJFILE(dao) == EOS)
+	    call dp_sparam (tp, "REJFILE", "\"\"", "filename", Memc[dummy])
+	else {
+	    call dp_froot (DP_OUTREJFILE(dao), Memc[outstr], SZ_LINE)
+	    call dp_sparam (tp, "REJFILE", Memc[outstr], "filename",
+	        Memc[dummy])
+	}
+	call dp_imroot (DP_OUTIMAGE(dao), Memc[outstr], SZ_LINE)
+	call dp_sparam (tp, "SUBIMAGE", Memc[outstr], "imagename",
+	    Memc[dummy])
+
+	# Define the data characteristics.
+	call dp_rparam (tp, "SCALE", DP_SCALE(dao), "units/pix", Memc[dummy])
+	call dp_rparam (tp, "DATAMIN", DP_MINGDATA(dao), "counts", Memc[dummy])
+	call dp_rparam (tp, "DATAMAX", DP_MAXGDATA(dao), "counts", Memc[dummy])
+	call dp_rparam (tp, "GAIN", DP_PHOTADU(dao), "e-/adu", Memc[dummy])
+	call dp_rparam (tp, "READNOISE", DP_READNOISE(dao), "electrons",
+	    Memc[dummy])
+
+	# Define the observing parameters.
+	call dp_sparam (tp, "OTIME", DP_OTIME(dao), "timeunit", Memc[dummy])
+	call dp_rparam (tp, "XAIRMASS", DP_XAIRMASS(dao), "number", Memc[dummy])
+	call dp_sparam (tp, "IFILTER", DP_IFILTER(dao), "filter", Memc[dummy])
+
+	# Define the fitting parameters.
+	call dp_bparam (tp, "RECENTER", itob (DP_RECENTER(dao)), "switch",
+	    Memc[dummy])
+	call dp_bparam (tp, "GRPSKY", itob (DP_GROUPSKY(dao)), "switch",
+	    Memc[dummy])
+	call dp_bparam (tp, "FITSKY", itob (DP_FITSKY(dao)), "switch",
+	    Memc[dummy])
+	call dp_rparam (tp, "PSFMAG", DP_PSFMAG(psffit), "magnitude",
+	    Memc[dummy])
+	call dp_rparam (tp, "PSFRAD", DP_SPSFRAD(dao), "scaleunit", Memc[dummy])
+	call dp_rparam (tp, "FITRAD", DP_SFITRAD(dao), "scaleunit", Memc[dummy])
+	call dp_rparam (tp, "ANNULUS", DP_SANNULUS(dao), "scaleunit",
+	    Memc[dummy])
+	call dp_rparam (tp, "DANNULUS", DP_SDANNULUS(dao), "scaleunit",
+	    Memc[dummy])
+	call dp_iparam (tp, "MAXITER", DP_MAXITER(dao), "number", Memc[dummy])
+	call dp_iparam (tp, "MAXGROUP", DP_MAXGROUP(dao), "number", Memc[dummy])
+	#call dp_iparam (tp, "MAXNSTAR", DP_MAXNSTAR(dao), "number",
+	    #Memc[dummy])
+	call dp_rparam (tp, "FLATERROR", DP_FLATERR(dao), "percentage",
+	    Memc[dummy])
+	call dp_rparam (tp, "PROFERROR", DP_PROFERR(dao), "percentage",
+	    Memc[dummy])
+	call dp_iparam (tp, "CLIPEXP", DP_CLIPEXP(dao), "number", Memc[dummy])
+	call dp_rparam (tp, "CLIPRANGE", DP_CLIPRANGE(dao), "sigma",
+	    Memc[dummy])
+	call dp_rparam (tp, "MERGERAD", DP_SMERGERAD(dao), "scaleunit",
+	    Memc[dummy])
+
+	call sfree (sp)
+end
+
+
+# DP_TALWRITE -- Write out a photometry record to an ALLSTAR ST table.
+
+procedure dp_talwrite (tpout, tprout, colpoint, id, x, y, mag, magerr, sky,
+	chi, numer, denom, skip, aier, niter, istar, lstar, star, rstar,
+	psfmag, csharp)
+
+pointer	tpout		# the output photometry file descriptor
+pointer	tprout		# the output rejections file descriptor
+int	colpoint[ARB]	# array of column pointers
+int	id[ARB]		# array of ids
+real	x[ARB]		# array of xpositions
+real	y[ARB]		# array of y positions
+real	mag[ARB]	# array of magnitude values
+real	magerr[ARB]	# array of mangitudes
+real	sky[ARB]	# array of sky values
+real	chi[ARB]	# array of chi values
+real	numer[ARB]	# array of first numerator values
+real	denom[ARB]	# array of first denominator values
+int	skip[ARB]	# array of skipped stars
+int	aier[ARB]	# array of error codes
+int	niter		# number of iterations
+int	istar, lstar	# first and last stars
+int	star		# photometry file row number
+int	rstar		# rejections file row number
+real	psfmag		# magnitude of the psf
+real	csharp		# the sharpness constant
+
+int	i, iter, pier, plen
+pointer	sp, perror
+real	err, sharp
+int	dp_gallerr()
+
+begin
+	call smark (sp)
+	call salloc (perror, SZ_FNAME, TY_CHAR)
+
+	do i = istar, lstar {
+
+	    # Skip the star ?
+	    if (skip[i] == NO)
+		next
+	    if (IS_INDEFR(x[i]) || IS_INDEFR(y[i]))
+		next
+
+	    # Get the results.
+	    if (IS_INDEFR(mag[i]) || mag[i] <= 0.0 || denom[i] == 0.) {
+		mag[i] = INDEFR
+		iter = 0
+		err = INDEFR
+		sharp = INDEFR
+		chi[i] = INDEFR
+	    } else {
+		if (magerr[i] <= 0.0)
+		    err = 0.0
+		else
+	            err = 1.085736 * magerr[i] / mag[i]
+		sharp = csharp * numer[i] / (mag[i] * denom[i])
+		sharp = max (MIN_SHARP, min (sharp, MAX_SHARP))
+	        mag[i] = psfmag - 2.5 * log10 (mag[i])
+		iter = niter
+	    }
+	    pier = aier[i]
+	    plen = dp_gallerr (pier, Memc[perror], SZ_FNAME)
+
+	    # Write the output row.
+	    if ((tprout != NULL) && (pier != ALLERR_OK)) {
+	        rstar = rstar + 1
+	        call tbrpti (tprout, colpoint[1], id[i], 1, rstar)
+	        call tbrptr (tprout, colpoint[2], x[i], 1, rstar)
+	        call tbrptr (tprout, colpoint[3], y[i], 1, rstar)
+	        call tbrptr (tprout, colpoint[4], mag[i], 1, rstar)
+	        call tbrptr (tprout, colpoint[5], err, 1, rstar)
+	        call tbrptr (tprout, colpoint[6], sky[i], 1, rstar)
+	        call tbrpti (tprout, colpoint[7], iter, 1, rstar)
+	        call tbrptr (tprout, colpoint[8], sharp, 1, rstar)
+	        call tbrptr (tprout, colpoint[9], chi[i], 1, rstar)
+	        call tbrpti (tprout, colpoint[10], pier, 1, rstar)
+		call tbrptt (tprout, colpoint[11], Memc[perror], plen,
+		    1, rstar)
+	    } else {
+	        star = star + 1
+	        call tbrpti (tpout, colpoint[1], id[i], 1, star)
+	        call tbrptr (tpout, colpoint[2], x[i], 1, star)
+	        call tbrptr (tpout, colpoint[3], y[i], 1, star)
+	        call tbrptr (tpout, colpoint[4], mag[i], 1, star)
+	        call tbrptr (tpout, colpoint[5], err, 1, star)
+	        call tbrptr (tpout, colpoint[6], sky[i], 1, star)
+	        call tbrpti (tpout, colpoint[7], iter, 1, star)
+	        call tbrptr (tpout, colpoint[8], sharp, 1, star)
+	        call tbrptr (tpout, colpoint[9], chi[i], 1, star)
+	        call tbrpti (tpout, colpoint[10], pier, 1, star)
+		call tbrptt (tpout, colpoint[11], Memc[perror], plen,
+		    1, star)
+	    }
+	}
+
+	call sfree (sp)
+end
+
+
+define AL_DATA1STR "%-9d%10t%-10.3f%20t%-10.3f%30t%-12.3f%42t%-14.3f%56t%-15.7g%71t%-6d%80t\\\n"
+define AL_DATA2STR "%12t%-12.3f%24t%-12.3f%36t%-6d%42t%-13.13s%80t \n"
+
+# DP_XALWRITE -- Write out photometry record to an ALLSTAR text file.
+
+procedure dp_xalwrite (tpout, tprout, id, x, y, mag, magerr, sky, chi, numer,
+	denom, skip, aier, niter, istar, lstar, psfmag, csharp)
+
+int	tpout		# the output photometry file descriptor
+int	tprout		# the output rejections file descriptor
+int	id[ARB]		# array of ids
+real	x[ARB]		# array of xpositions
+real	y[ARB]		# array of y positions
+real	mag[ARB]	# array of magnitude values
+real	magerr[ARB]	# array of magnitudes
+real	sky[ARB]	# array of sky values
+real	chi[ARB]	# array of chi values
+real	numer[ARB]	# array of first numerator values
+real	denom[ARB]	# array of first denominator values
+int	skip[ARB]	# array of skipped stars
+int	aier[ARB]	# array of error codes
+int	niter		# number of iterations
+int	istar, lstar	# first and last stars
+real	psfmag		# magnitude of the psf
+real	csharp		# sharpness constant
+
+int	i, iter, pier, plen
+pointer	sp, perror
+real	err, sharp
+int	dp_gallerr()
+
+begin
+	call smark (sp)
+	call salloc (perror, SZ_FNAME, TY_CHAR)
+
+	do i = istar, lstar {
+
+	    # Skip the star ?
+	    if (skip[i] == NO)
+		next
+	    if (IS_INDEFR(x[i]) || IS_INDEFR(y[i]))
+		next
+
+	    if (IS_INDEFR(mag[i]) || mag[i] <= 0.0 || denom[i] == 0.) {
+		mag[i] = INDEFR
+		iter = 0
+		err = INDEFR
+		sharp = INDEFR
+		chi[i] = INDEFR
+	    } else {
+		if (magerr[i] <= 0.0)
+		    err = 0.0
+		else
+	            err = 1.085736 * magerr[i] / mag[i]
+		sharp = csharp * numer[i] / (mag[i] * denom[i])
+		sharp = max (MIN_SHARP, min (sharp, MAX_SHARP))
+	        mag[i] = psfmag - 2.5 * log10 (mag[i])
+	        iter = niter
+	    }
+	    pier = aier[i]
+	    plen = dp_gallerr (pier, Memc[perror], SZ_FNAME)
+
+	    if ((tprout != NULL) && (pier != ALLERR_OK)) {
+	        call fprintf (tprout, AL_DATA1STR)
+	            call pargi (id[i])
+	            call pargr (x[i])
+	            call pargr (y[i])
+	            call pargr (mag[i])
+	            call pargr (err)
+	            call pargr (sky[i])
+	            call pargi (iter)
+	        call fprintf (tprout, AL_DATA2STR)
+	            call pargr (sharp)
+	            call pargr (chi[i])
+		    call pargi (pier)
+		    call pargstr (Memc[perror])
+	    } else {
+	        call fprintf (tpout, AL_DATA1STR)
+	            call pargi (id[i])
+	            call pargr (x[i])
+	            call pargr (y[i])
+	            call pargr (mag[i])
+	            call pargr (err)
+	            call pargr (sky[i])
+	            call pargi (iter)
+	        call fprintf (tpout, AL_DATA2STR)
+	            call pargr (sharp)
+	            call pargr (chi[i])
+		    call pargi (pier)
+		    call pargstr (Memc[perror])
+	    }
+	}
+
+	call sfree (sp)
+end
+
+
+# DP_GALLERR -- Set the ALLSTAR task error code string.
+ 
+int procedure dp_gallerr (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 ALLERR_OK:
+            plen = gstrcpy ("No_error", perror, maxch)
+	case ALLERR_BIGGROUP:
+            plen = gstrcpy ("Big_group", perror, maxch)
+        case ALLERR_INDEFSKY:
+            plen = gstrcpy ("Bad_sky", perror, maxch)
+        case ALLERR_NOPIX:
+            plen = gstrcpy ("Npix_too_few", perror, maxch)
+        case ALLERR_SINGULAR:
+            plen = gstrcpy ("Singular", perror, maxch)
+        case ALLERR_FAINT:
+            plen = gstrcpy ("Too_faint", perror, maxch)
+        case ALLERR_MERGE:
+            plen = gstrcpy ("Merged", perror, maxch)
+        case ALLERR_OFFIMAGE:
+            plen = gstrcpy ("Off_image", perror, maxch)
+        default:
+            plen = gstrcpy ("No_error", perror, maxch)
+        }
+ 
+        return (plen)
+end
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
diff --git a/noao/digiphot/daophot/allstar/dpcache.x b/noao/digiphot/daophot/allstar/dpcache.x
new file mode 100644
index 00000000..0782b98f
--- /dev/null
+++ b/noao/digiphot/daophot/allstar/dpcache.x
@@ -0,0 +1,244 @@
+include	<imhdr.h>
+include <mach.h>
+include	"../lib/allstardef.h"
+include "../lib/daophotdef.h"
+
+define	FUDGE	1.2		# fudge factor for memory allocation
+
+# DP_CACHE -- Determine whether it is possible to store all the data
+# in memory or not.
+
+procedure dp_cache (dao, im, subim, cache)
+
+pointer	dao			# pointer to the daophot structure
+pointer	im			# pointer to the input image
+pointer	subim			# pointer to the output subtracted image
+int	cache			# cache the data ?
+
+int	npix, req_mem1, req_mem2, req_mem3, old_mem, max_mem
+pointer	sp, temp, allstar, data, weights, subt
+int	sizeof(), begmem()
+pointer	immap()
+errchk	begmem(), calloc()
+
+begin
+	# Allocate some working memory.
+	call smark (sp)
+	call salloc (temp, SZ_FNAME, TY_CHAR)
+
+	# Define the allstar pointer.
+	allstar = DP_ALLSTAR(dao)
+
+	# Store the old working set size.
+	DP_SZOLDCACHE(allstar) = begmem (0, old_mem, max_mem)
+
+	# Figure out the memory requirements.
+	npix = IM_LEN(im,1) * IM_LEN(im,2)
+	req_mem1 = FUDGE *  npix * sizeof (TY_REAL)
+	req_mem2 = req_mem1 + req_mem1
+	req_mem3 = req_mem2 + req_mem1
+
+	# Initialize.
+	DP_CACHE (allstar,A_DCOPY) = NO
+	DP_CACHE (allstar,A_SUBT) = NO
+	DP_CACHE (allstar,A_WEIGHT) = NO
+
+	DP_DBUF(allstar) = NULL
+	DP_SBUF(allstar) = NULL
+	DP_WBUF(allstar) = NULL
+
+	# Use IRAF images as temporary storage space.
+	if (cache == NO) {
+
+	    # The output subtracted image.
+	    DP_DATA (allstar) = subim
+
+	    # The scratch image.
+	    call mktemp ("subt", Memc[temp], SZ_FNAME)
+	    DP_SUBT (allstar) = immap (Memc[temp], NEW_COPY, im)
+	    IM_NDIM(DP_SUBT(allstar)) = 2
+	    IM_PIXTYPE(DP_SUBT(allstar)) = TY_REAL
+
+	    # The weights image.
+	    call mktemp ("wt", Memc[temp], SZ_FNAME)
+	    DP_WEIGHTS (allstar) = immap (Memc[temp], NEW_COPY, im)
+	    IM_NDIM(DP_WEIGHTS(allstar)) = 2
+	    IM_PIXTYPE(DP_WEIGHTS(allstar)) = TY_REAL
+
+	    # Set the type of caching to be used.
+	    DP_ISCACHE(allstar) = NO
+	    DP_ALLCACHE(allstar) = NO
+
+	} else {
+
+	    # Is the memory available.
+	    if (old_mem >= req_mem3) {
+		DP_CACHE(allstar,A_DCOPY) = YES
+		DP_CACHE(allstar,A_WEIGHT) = YES
+	        DP_CACHE (allstar,A_SUBT) = YES
+	    } else {
+	        if (begmem (req_mem1, old_mem, max_mem) >= req_mem1)
+		    DP_CACHE(allstar,A_SUBT) = YES
+	        if (begmem (req_mem2, old_mem, max_mem) >= req_mem2)
+		    DP_CACHE(allstar,A_WEIGHT) = YES
+	        if (begmem (req_mem3, old_mem, max_mem) >= req_mem3)
+		    DP_CACHE(allstar,A_DCOPY) = YES
+	    }
+
+	    # Allocate space for the scratch image.
+	    subt = NULL
+	    if (DP_CACHE(allstar, A_SUBT) == YES) {
+		iferr {
+	            call calloc (subt, npix, TY_REAL)
+		} then {
+		    if (subt != NULL)
+		        call mfree (subt, TY_REAL)
+		    call mktemp ("subt", Memc[temp], SZ_FNAME)
+		    subt = immap (Memc[temp], NEW_COPY, im)
+		    IM_NDIM(subt) = 2
+		    IM_PIXTYPE(subt) = TY_REAL
+		    DP_SUBT(allstar) = subt
+		    DP_CACHE(allstar,A_SUBT) = NO
+		} else
+		    DP_SUBT(allstar) = subt
+	    } else {
+		call mktemp ("subt", Memc[temp], SZ_FNAME)
+		subt = immap (Memc[temp], NEW_COPY, im)
+		IM_NDIM(subt) = 2
+		IM_PIXTYPE(subt) = TY_REAL
+		DP_SUBT(allstar) = subt
+	    }
+
+	    # Allocate space for the weights image.
+	    weights = NULL
+	    if (DP_CACHE(allstar, A_WEIGHT) == YES) {
+		iferr {
+	            call calloc (weights, npix, TY_REAL)
+		} then {
+		    if (weights != NULL)
+		        call mfree (weights, TY_REAL)
+		    call mktemp ("wt", Memc[temp], SZ_FNAME)
+		    weights = immap (Memc[temp], NEW_COPY, im)
+		    IM_NDIM(weights) = 2
+		    IM_PIXTYPE(weights) = TY_REAL
+		    DP_WEIGHTS(allstar) = weights
+		    DP_CACHE(allstar,A_WEIGHT) = NO
+		} else
+		    DP_WEIGHTS(allstar) = weights
+	    } else {
+		call mktemp ("wt", Memc[temp], SZ_FNAME)
+		weights = immap (Memc[temp], NEW_COPY, im)
+		IM_NDIM(weights) = 2
+		IM_PIXTYPE(weights) = TY_REAL
+		DP_WEIGHTS(allstar) = weights
+	    }
+
+	    # Allocate space for the output subtracted image.
+	    data = NULL
+	    if (DP_CACHE(allstar, A_DCOPY) == YES) {
+		iferr {
+	            call calloc (data, npix, TY_REAL)
+		} then {
+		    if (data != NULL)
+		        call mfree (data, TY_REAL)
+		    DP_DATA(allstar) = subim
+		    DP_CACHE(allstar,A_DCOPY) = NO
+		} else {
+		    DP_DATA(allstar) = data
+		}
+	    } else
+		DP_DATA(allstar) = subim
+
+	    # Set the type of caching to be used.
+	    if (DP_CACHE(allstar,A_DCOPY) == NO && DP_CACHE(allstar,A_SUBT) ==
+		NO && DP_CACHE(allstar,A_WEIGHT) == NO) {
+		DP_ISCACHE(allstar) = NO
+		DP_ALLCACHE(allstar) = NO
+	    } else if (DP_CACHE(allstar,A_DCOPY) == YES && DP_CACHE(allstar,
+	        A_SUBT) == YES && DP_CACHE(allstar,A_WEIGHT) == YES) {
+		DP_ISCACHE(allstar) = YES
+		DP_ALLCACHE(allstar) = YES
+	    } else {
+		DP_ISCACHE(allstar) = YES
+		DP_ALLCACHE(allstar) = NO
+	    }
+
+	}
+
+	call sfree (sp)
+end
+
+
+# DP_UNCACHE -- Release all the stored memory.
+
+procedure dp_uncache (dao, subim, savesub)
+
+pointer	dao		# pointer to the daophot strucuture
+pointer	subim		# pointer to the output subtracted image
+int	savesub		# save the subtracted image ?
+
+int	j, ncol, nline
+pointer	sp, imname, v, data, buf, allstar
+pointer	impnlr()
+
+begin
+	# Allocate working space.
+	call smark (sp)
+	call salloc (imname, SZ_FNAME, TY_CHAR)
+	call salloc (v, IM_MAXDIM, TY_LONG)
+
+	# Define the allstar pointer.
+	allstar = DP_ALLSTAR(dao)
+
+	ncol = IM_LEN(subim,1)
+	nline = IM_LEN(subim,2)
+
+	# Write out the subtracted image.
+	if (DP_CACHE(allstar,A_DCOPY) == YES && savesub == YES) {
+	    data = DP_DBUF(allstar)
+	    call amovkl (long(1), Meml[v], IM_MAXDIM)
+	    do j = 1, nline {
+		if (impnlr (subim, buf, Meml[v]) != EOF)
+		    call amovr (Memr[data], Memr[buf], ncol)
+	        data = data + ncol
+	    }
+	}
+	DP_DATA(allstar) = NULL
+
+	# Release any memory used by the subtracted image.
+	if (DP_DBUF(allstar) != NULL)
+	    call mfree (DP_DBUF(allstar), TY_REAL)
+	DP_DBUF(allstar) = NULL
+
+	# Delete the scratch image if any.
+	if (DP_CACHE(allstar,A_SUBT) == NO) {
+	    call strcpy (IM_HDRFILE(DP_SUBT(allstar)), Memc[imname], SZ_FNAME)
+	    call imunmap (DP_SUBT(allstar)) 
+	    call imdelete (Memc[imname])
+	}
+	DP_SUBT(allstar) = NULL
+
+	# Release any memory used by the scratch image.
+	if (DP_SBUF(allstar) != NULL)
+	    call mfree (DP_SBUF(allstar), TY_REAL)
+	DP_SBUF(allstar) = NULL
+
+	# Delete the weights image if any.
+	if (DP_CACHE(allstar,A_WEIGHT) == NO) {
+	    call strcpy (IM_HDRFILE(DP_WEIGHTS(allstar)), Memc[imname],
+	        SZ_FNAME)
+	    call imunmap (DP_WEIGHTS(allstar)) 
+	    call imdelete (Memc[imname])
+	}
+	DP_WEIGHTS(allstar) = NULL
+
+	# Release any memory used by the weights buffer.
+	if (DP_WBUF(allstar) != NULL)
+	    call mfree (DP_WBUF(allstar), TY_REAL)
+	DP_WBUF(allstar) = NULL
+
+	# Reset the working set size to the previous value.
+	call fixmem (DP_SZOLDCACHE(allstar))
+
+	call sfree (sp)
+end
diff --git a/noao/digiphot/daophot/allstar/dpglim.x b/noao/digiphot/daophot/allstar/dpglim.x
new file mode 100644
index 00000000..d5f1fedb
--- /dev/null
+++ b/noao/digiphot/daophot/allstar/dpglim.x
@@ -0,0 +1,17 @@
+# DP_GLIM -- Get the lower and upper limits of a section around a specified
+# center.
+
+procedure dp_glim (xc, yc, radius, ixmin, ixmax, iymin, iymax, lx, mx, ly, my)
+
+real	xc, yc			# the x and y center points
+real	radius			# the radial distance
+int	ixmin, ixmax		# absolute x boundaries
+int	iymin, iymax		# absolute y boundaries
+int	lx, mx, ly, my		# the returned limits
+
+begin
+        lx = max (ixmin - 1, min (ixmax, int (xc - radius))) + 1
+        mx = max (ixmin, min (ixmax, int (xc + radius)))
+        ly = max (iymin - 1, min (iymax, int (yc - radius))) + 1
+        my = max (iymin, min (iymax, int (yc + radius)))
+end
diff --git a/noao/digiphot/daophot/allstar/dprectify.x b/noao/digiphot/daophot/allstar/dprectify.x
new file mode 100644
index 00000000..5fd7c0db
--- /dev/null
+++ b/noao/digiphot/daophot/allstar/dprectify.x
@@ -0,0 +1,74 @@
+# DP_RECTIFY -- Shuffle a real vector based upon an input index vector using
+# dynamic memory.
+
+procedure dp_rectify (x, index, nitem)
+
+real	x[ARB]
+int	index[ARB]
+int	nitem
+
+pointer	sp, hold
+
+begin
+	call smark (sp)
+	call salloc (hold, nitem, TY_REAL)
+	call dp_dorect (x, Memr[hold], index, nitem)
+	call sfree (sp)
+end
+
+
+# DP_IRECTIFY -- Shuffle an integer vector based upon an input index vector
+# using dynamic memory.
+
+procedure dp_irectify (x, index, nitem)
+
+int	x[ARB]
+int	index[ARB]
+int	nitem
+
+pointer	sp, hold
+
+begin
+	call smark (sp)
+	call salloc (hold, nitem, TY_INT)
+	call dp_idorect (x, Memi[hold], index, nitem)
+	call sfree (sp)
+end
+
+
+# DP_DORECT -- Shuffle a vector based upon an input index vector using
+# a preallocated dummy array.
+
+procedure dp_dorect (x, hold, index, nitem)
+
+real	x[ARB]
+real	hold[ARB]
+int	index[ARB]
+int	nitem
+
+int	i
+
+begin
+	call amovr (x, hold, nitem)
+	do i = 1, nitem
+	    x[i] = hold[index[i]]
+end
+
+
+# DP_IDORECT -- Shuffle an integer  vector based upon an input index vector
+# using a preallocated dummy array.
+
+procedure dp_idorect (x, hold, index, nitem)
+
+int	x[ARB]
+int	hold[ARB]
+int	index[ARB]
+int	nitem
+
+int	i
+
+begin
+	call amovi (x, hold, nitem)
+	do i = 1, nitem
+	    x[i] = hold[index[i]]
+end
diff --git a/noao/digiphot/daophot/allstar/dpregroup.x b/noao/digiphot/daophot/allstar/dpregroup.x
new file mode 100644
index 00000000..a7561fb8
--- /dev/null
+++ b/noao/digiphot/daophot/allstar/dpregroup.x
@@ -0,0 +1,219 @@
+# DP_ALSORT -- Sort the stars on the y coordinate.
+
+procedure dp_alsort (id, x, y, mag, sky, sumwt, dxold, dyold, xclamp, yclamp,
+	nstar)
+
+int	id[ARB]			# array if star ids
+real	x[ARB]			# array of star x values
+real	y[ARB]			# array of y values
+real	mag[ARB]		# array of star magnitudes
+real	sky[ARB]		# array of star sky values
+real	sumwt[ARB]		# array of weighted sums
+real	dxold[ARB]		# the previous star x array
+real	dyold[ARB]		# the previous star y array
+real	xclamp[ARB]		# the x array clamps
+real	yclamp[ARB]		# the y array clamps
+int	nstar			# the number of stars
+
+int	ier
+pointer	sp, index
+
+begin
+	# Allocate some working memory.
+	call smark (sp)
+	call salloc (index, nstar, TY_INT)
+
+	# Sort the star list on y.
+	call quick (y, nstar, Memi[index], ier)
+
+	# Recitfy the remaining arrays.
+	call dp_irectify (id, Memi[index], nstar)
+	call dp_rectify (x, Memi[index], nstar)
+	call dp_rectify (mag, Memi[index], nstar)
+	call dp_rectify (sky, Memi[index], nstar)
+	call dp_rectify (sumwt, Memi[index], nstar)
+	call dp_rectify (dxold, Memi[index], nstar)
+	call dp_rectify (dyold, Memi[index], nstar)
+	call dp_rectify (xclamp, Memi[index], nstar)
+	call dp_rectify (yclamp, Memi[index], nstar)
+
+	# Release memory.
+	call sfree (sp)
+end
+
+
+# DP_REGROUP -- Group stars into physical associations based on proximity.
+
+procedure dp_regroup (id, x, y, mag, sky, chi, dxold, dyold, xclamp, yclamp,
+	nstar, radius, last)
+
+int	id[ARB]			# array if star ids
+real	x[ARB]			# array of star x values
+real	y[ARB]			# array of y values
+real	mag[ARB]		# array of star magnitudes
+real	sky[ARB]		# array of star sky values
+real	chi[ARB]		# array of chis
+real	dxold[ARB]		# the previous star x array
+real	dyold[ARB]		# the previous star y array
+real	xclamp[ARB]		# the x array clamps
+real	yclamp[ARB]		# the y array clamps
+int	nstar			# the number of stars
+real	radius			# the fitting radius
+int	last[ARB]		# the last array (NO or YES for last star)
+
+int	itop, itest, j, k, i, ier
+pointer sp, index
+real	critrad, critradsq, xtest, ytest, dx, dy
+
+begin
+	# If there is only one stars its value of last is YES.
+	if (nstar <= 1) {
+	    last[1] = YES
+	    return
+	}
+
+	# Allocate some working memory.
+	call smark (sp)
+	call salloc (index, nstar, TY_INT)
+
+	# Compute the critical radius.
+	critrad = 2.0 * radius
+	critradsq = critrad * critrad
+
+	# Sort the star list on y and rectify the accompanying x array.
+	call quick (y, nstar, Memi[index], ier)
+	call dp_rectify (x, Memi[index], nstar)
+
+	# At this point the stars are in a stack NSTAR stars long. The
+	# variable ITEST will point to the position in the stack
+	# occupied by the star which is currently the center of a
+	# circle of the CRITRAD pixels, within which we are 
+	# looking for other stars. ITEST starts out with a value of one.
+	# ITOP points to the top position in the stack of the 
+	# stars which have not yet been assigned to groups. ITOP starts
+	# out with the value of 2. Each time through, the program goes 
+	# down through the stack from ITOP and looks for stars within 
+	# the critrad pixels from the star at stack position ITEST.
+	# When such a star is found, it changes places in the 
+	# stack with the star at ITOP and ITOP is incremented by one.
+	# When the search reaches a star J such that Y(J) - Y(ITEST)
+	# > CRITRAD we know that no further stars will be found 
+	# within the CRITRAD pixels, the pointer ITEST is incremented
+	# by one, and the search proceeds again from the new
+	# value of ITOP. If the pointer ITEST catches up with the
+	# pointer ITOP, then the current group is complete
+	# with the star at the position ITEST-1, and
+	# ITOP is incremented by one. If ITOP catches up with NSTAR
+	# the grouping process is complete.
+
+	itop = 2
+	do itest = 1, nstar {
+
+	    # Initialize the last array.
+	    last[itest] = NO
+
+	    # ITEST has reached ITOP; so no other unassigned stars are
+	    # within CRITRADIUS pixels of any member of the current
+	    # group. The group is therefore complete. Signify this by
+	    # setting LAST[ITEST-1] = YES, and incrementing the value of
+	    # ITOP by one.  If ITOP is greater than NSTAR at this point
+	    # list is finished.
+
+	    if (itest == itop) {
+		j = itest - 1
+		if (j > 0)
+		    last[j] = YES
+		itop = itop + 1
+		if (itop > nstar) {
+		    last[itest] = YES
+		 	break
+		 }
+	    }
+
+	    # Now go through the list of unassigned stars, occupying
+	    # positions ITOP through NSTAR in the stack, to look for
+	    # stars within CRITRADIUS pixels of the star at 
+	    # position ITEST inthe stack. If one is found, move it
+ 	    # up to stack position ITOP and increment ITOP by one.
+
+	    xtest = x[itest]
+	    ytest = y[itest]
+	    j = itop
+
+	    do i = j, nstar {
+
+		# Check the distance between the stars.
+		dy = y[i] - ytest
+		if (dy > critrad)
+		    break
+		dx = x[i] - xtest
+		if (abs (dx) > critrad)
+		    next
+		if ((dx * dx + dy * dy) > critradsq)
+		    next
+
+		# This star is within CRITRAD pixels of the
+		# star at stack position ITEST. Therefore it is
+		# added to the current group by moving it up to position
+		# ITOP in the stack, where the pointer ITEST may reach it.
+
+		call dp_rgswap (itop, i, x, y, Memi[index])
+
+		# Now increment ITOP by 1 to point at the top most 
+		# unassigned star in the stack.
+		itop = itop + 1
+		if (itop > nstar) {
+		    do k = itest, nstar - 1
+			last[k] = NO
+		    last[nstar] = YES
+		    break
+		}
+	    }
+
+	    if (itop > nstar)
+		break
+	}
+
+	# Reorder the remaining arrays to match the x and y arrays.
+	call dp_irectify (id, Memi[index], nstar)
+   	call dp_rectify (mag, Memi[index], nstar)
+	call dp_rectify (sky, Memi[index], nstar)
+	call dp_rectify (chi, Memi[index], nstar)
+	call dp_rectify (dxold, Memi[index], nstar)
+	call dp_rectify (dyold, Memi[index], nstar)
+	call dp_rectify (xclamp, Memi[index], nstar)
+	call dp_rectify (yclamp, Memi[index], nstar)
+
+	call sfree (sp)
+end
+
+
+# DP_RGSWAP -- Swap the i-th and j-th stars in the stack without otherwise
+# altering the order of the stars..
+
+procedure dp_rgswap (i, j, x, y, index)
+
+int	i, j		# indices of the two stars to be swapped
+real	x[ARB]		# the array of x values
+real	y[ARB]		# the array of y values
+int	index[ARB]	# the index array
+
+int	ihold, k,l
+real	xhold, yhold
+
+begin
+	xhold = x[j]
+	yhold = y[j]
+	ihold = index[j]
+
+	do k = j, i + 1, -1 {
+	    l = k - 1
+ 	    x[k] = x[l]
+	    y[k] = y[l]
+	    index[k] = index[l]
+	}
+
+	x[i] = xhold
+	y[i] = yhold
+	index[i] = ihold
+end
diff --git a/noao/digiphot/daophot/allstar/mkpkg b/noao/digiphot/daophot/allstar/mkpkg
new file mode 100644
index 00000000..d607592b
--- /dev/null
+++ b/noao/digiphot/daophot/allstar/mkpkg
@@ -0,0 +1,32 @@
+# ALLSTAR task
+
+$checkout libpkg.a ".."
+$update libpkg.a
+$checkin libpkg.a ".."
+$exit
+
+libpkg.a:
+	dpabuf.x	<imhdr.h>		../lib/daophotdef.h       \
+			../lib/allstardef.h
+	dpaconfirm.x	../lib/daophotdef.h
+	dpcache.x  	 <imhdr.h>		<mach.h>                  \
+			 ../lib/daophotdef.h	../lib/allstardef.h
+	dpregroup.x
+	dprectify.x
+	dpalwrite.x	<tbset.h>               <time.h>                  \
+			../lib/daophotdef.h     ../lib/apseldef.h         \
+			../lib/allstardef.h
+	dpastar.x	<mach.h>                <imhdr.h>                 \
+			../lib/daophotdef.h	../lib/allstardef.h       \
+			../lib/apseldef.h
+	dpalphot.x	<mach.h>		<imhdr.h>		  \
+			../lib/daophotdef.h     ../lib/apseldef.h         \
+			../lib/allstardef.h
+	dpalinit.x 	<imhdr.h>		<mach.h>                  \
+			<math.h>                ../lib/daophotdef.h       \
+			../lib/allstardef.h
+	dpglim.x
+	dpalmemstar.x	../lib/daophotdef.h	../lib/allstardef.h
+	t_allstar.x  	<fset.h>		../lib/daophotdef.h       \
+			<imhdr.h>
+	;
diff --git a/noao/digiphot/daophot/allstar/t_allstar.x b/noao/digiphot/daophot/allstar/t_allstar.x
new file mode 100644
index 00000000..b3991179
--- /dev/null
+++ b/noao/digiphot/daophot/allstar/t_allstar.x
@@ -0,0 +1,355 @@
+include <imhdr.h>
+include	<fset.h>
+include "../lib/daophotdef.h"
+
+# T_ALLSTAR  -- Group a list of stars into physical associations, fit the PSF
+# to the stars in each group simultaneously, and produce an output subtracted
+# image.
+
+procedure t_allstar ()
+
+pointer	image			# the input image
+pointer	psfimage		# the input psf image
+pointer	photfile		# the input photometry file
+pointer	allstarfile		# the output photometry file
+pointer	subimage		# the output subtracted image
+pointer	rejfile			# the output rejections file
+int	cache			# cache the data in memory
+int	verbose			# verbose mode ?
+int	verify			# verify critical task parameters ?
+int	update			# update the task parameters ?
+int	version			# version number
+
+pointer	sp, outfname, im, subim, dao, str
+int	imlist, limlist, alist, lalist, pimlist, lpimlist, olist, lolist
+int	simlist, lsimlist, rlist, lrlist, photfd, psffd, allfd, root, savesub
+int	rejfd, wcs
+bool	ap_text
+
+pointer	immap(), tbtopn()
+int	clgwrd(), clgeti()
+int	open(), fnldir(), strncmp(), strlen(), btoi(), access()
+int	fstati(), imtopen, imtlen(), imtgetim(), fntopnb(), fntlenb()
+int	fntgfnb()
+bool	itob(), clgetb()
+
+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 (psfimage, SZ_FNAME, TY_CHAR)
+	call salloc (photfile, SZ_FNAME, TY_CHAR)
+	call salloc (allstarfile, SZ_FNAME, TY_CHAR)
+	call salloc (rejfile, SZ_FNAME, TY_CHAR)
+	call salloc (subimage, SZ_FNAME, TY_CHAR)
+	call salloc (outfname, SZ_FNAME, TY_CHAR)
+	call salloc (str, SZ_FNAME, TY_CHAR)
+
+	# Get the task input and output file names.
+	call clgstr ("image", Memc[image], SZ_FNAME)
+ 	call clgstr ("photfile", Memc[photfile], SZ_FNAME)
+	call clgstr ("psfimage", Memc[psfimage], SZ_FNAME)
+	call clgstr ("allstarfile", Memc[allstarfile], SZ_FNAME)
+	call clgstr ("subimage", Memc[subimage], SZ_FNAME)
+	call clgstr ("rejfile", Memc[rejfile], SZ_FNAME)
+
+	# Get the task mode parameters.
+	verbose = btoi (clgetb ("verbose"))
+	verify = btoi (clgetb ("verify"))
+	update = btoi (clgetb ("update"))
+	cache = btoi (clgetb ("cache"))
+	version = clgeti ("version")
+	version = 2
+
+	# 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[allstarfile], NO)
+	lolist = fntlenb (olist)
+	simlist = imtopen (Memc[subimage])
+	lsimlist = imtlen (simlist)
+	rlist = fntopnb (Memc[rejfile], NO)
+	lrlist = fntlenb (rlist)
+
+	# Test the lengths of the photometry file, psf image and subtracted
+	# image lists are the same as the length of the input image.
+
+	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 imtclose (simlist)
+	    call sfree (sp)
+	    call error (0,
+	        "Incompatible image and 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 imtclose (simlist)
+	    call sfree (sp)
+	    call error (0,
+	        "Incompatible image and psf file list lengths")
+	}
+
+	if ((limlist != lsimlist) && (strncmp (Memc[subimage], DEF_DEFNAME,
+	    DEF_LENDEFNAME) != 0)) {
+	    call imtclose (imlist)
+	    call fntclsb (alist)
+	    call imtclose (pimlist)
+	    call fntclsb (olist)
+	    call fntclsb (rlist)
+	    call imtclose (simlist)
+	    call sfree (sp)
+	    call error (0,
+	        "Incompatible image and subtracted image list lengths")
+	}
+
+	if ((limlist != lolist) && (strncmp (Memc[allstarfile], DEF_DEFNAME,
+	    DEF_LENDEFNAME) != 0)) {
+	    call imtclose (imlist)
+	    call fntclsb (alist)
+	    call imtclose (pimlist)
+	    call fntclsb (olist)
+	    call fntclsb (rlist)
+	    call imtclose (simlist)
+	    call sfree (sp)
+	    call error (0,
+	        "Incompatible image and subtracted image list lengths")
+	}
+
+	if ((lrlist > 0) && (limlist != lrlist) && (strncmp (Memc[rejfile],
+	    DEF_DEFNAME, DEF_LENDEFNAME) != 0)) {
+	    call imtclose (imlist)
+	    call fntclsb (alist)
+	    call imtclose (pimlist)
+	    call fntclsb (olist)
+	    call fntclsb (rlist)
+	    call imtclose (simlist)
+	    call sfree (sp)
+	    call error (0,
+	        "Incompatible image and subtracted image list lengths")
+	}
+
+	# Open the input image, initialize the daophot structure, get the
+	# pset parameters
+
+	call dp_gppars (dao)	
+
+	# Set some parameters.
+	call dp_seti (dao, VERBOSE, verbose)
+
+	# Verify and update the parameters as appropriate.
+	if (verify == YES) {
+	    call dp_aconfirm (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 daophot fitting structure.
+	call dp_apsetup (dao)
+
+	# Initialize the allstar structure.
+	call dp_allstarsetup (dao)
+
+	# Loop over the images. 
+	while (imtgetim (imlist, Memc[image], SZ_FNAME) != EOF) {
+
+	    # Open the image and store some header parameters.
+	    im = immap (Memc[image], READ_ONLY, 0)		
+	    call dp_imkeys (dao, im)
+	    call dp_sets (dao, INIMAGE, Memc[image])
+
+	    # 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)
+	        photfd = open (Memc[outfname], READ_ONLY, TEXT_FILE)
+	    else
+	        photfd = tbtopn (Memc[outfname], READ_ONLY, 0)
+	    call dp_sets (dao, INPHOTFILE, Memc[outfname])
+	    call dp_wgetapert (dao, im, photfd, DP_MAXNSTAR(dao), ap_text)
+
+	    # Open the PSF image and read in the PSF.
+	    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)
+	    psffd = immap (Memc[outfname], READ_ONLY, 0)
+	    call dp_sets (dao, PSFIMAGE, Memc[outfname])
+	    call dp_readpsf (dao, psffd)
+	
+	    # Open the output ALLSTAR file. If the output is DEF_DEFNAME,
+	    # dir$default or a directory specification then the extension
+	    # "als" is added to the image name and a suitable version
+	    # number if appended to the output name.
+
+	    if (fntgfnb (olist, Memc[allstarfile], SZ_FNAME) == EOF)
+		call strcpy (DEF_DEFNAME, Memc[allstarfile], SZ_FNAME)
+	    root = fnldir (Memc[allstarfile], Memc[outfname], SZ_FNAME)
+	    if (strncmp (DEF_DEFNAME, Memc[allstarfile+root],
+	        DEF_LENDEFNAME) == 0 || root == strlen (Memc[allstarfile]))
+	        call dp_outname (Memc[image], Memc[outfname], "als",
+		    Memc[outfname], SZ_FNAME)
+	    else 
+	        call strcpy (Memc[allstarfile], Memc[outfname], SZ_FNAME)
+	    if (DP_TEXT(dao)  == YES)
+	        allfd = open (Memc[outfname], NEW_FILE, TEXT_FILE)
+	    else
+	        allfd = tbtopn (Memc[outfname], NEW_FILE, 0)
+	    call dp_sets (dao, OUTPHOTFILE, Memc[outfname])
+
+	    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], "arj",
+		        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])
+
+	    # Open the subtracted image.
+	    savesub = YES
+	    if (imtgetim (simlist, Memc[subimage], SZ_FNAME) == EOF)
+		call strcpy (DEF_DEFNAME, Memc[subimage], SZ_FNAME)
+	    root = fnldir (Memc[subimage], Memc[outfname], SZ_FNAME)
+	    if (strncmp (DEF_DEFNAME, Memc[subimage+root],
+	        DEF_LENDEFNAME) == 0 || root == strlen (Memc[subimage])) {
+	        call dp_oimname (Memc[image], Memc[outfname], "sub",
+		    Memc[outfname], SZ_FNAME)
+	        subim = immap (Memc[outfname], NEW_COPY, im)
+		IM_NDIM(subim) = 2
+	        IM_PIXTYPE(subim) = TY_REAL
+	    } else {
+	        call strcpy (Memc[subimage], Memc[outfname], SZ_FNAME)
+	        subim = immap (Memc[outfname], NEW_COPY, im)
+		IM_NDIM(subim) = 2
+	        IM_PIXTYPE(subim) = TY_REAL
+	    }
+	    call dp_sets (dao, OUTIMAGE, Memc[outfname])
+
+	    # Fit the stars.
+	    call dp_astar (dao, im, subim, allfd, rejfd, cache, savesub,
+	        version)
+
+	    # Close the input image.
+	    call imunmap (im)
+
+	    # Close the input photometry file.
+	    if (ap_text)
+		call close (photfd)
+	    else
+	        call tbtclo (photfd)
+
+	    # Close PSF image.
+	    call imunmap (psffd)
+
+	    # Close the output photometry file.
+	    if (DP_TEXT(dao) == YES)
+		call close (allfd)
+	    else
+	        call tbtclo (allfd)
+
+	    # Close the output rejections files.
+	    if (rejfd != NULL) {
+	        if (DP_TEXT(dao) == YES)
+		    call close (rejfd)
+	        else
+	            call tbtclo (rejfd)
+	    }
+
+	    # Close the output subtracted image.
+	    call strcpy (IM_HDRFILE(subim), Memc[subimage], SZ_FNAME)
+	    call imunmap (subim)
+	    if (savesub == NO)
+	        call imdelete (Memc[subimage])
+	}
+
+	# Close the file/image lists.
+	call imtclose (imlist)
+	call fntclsb (alist)
+	call imtclose (pimlist)
+	call fntclsb (olist)
+	call fntclsb (rlist)
+	call imtclose (simlist)
+
+	# Close the allstar structure.
+	call dp_alclose (dao)
+
+	# Close the photometry structure.
+	call dp_apclose (dao)
+
+	# Close the PSF structure.
+	call dp_fitclose (dao)
+
+	# Close up the daophot structure.
+	call dp_free (dao)
+
+	call sfree(sp)
+end