Initial commit

1 files changed, 1383 insertions, 0 deletions

diff --git a/noao/digiphot/daophot/nstar/dpnstarfit.x b/noao/digiphot/daophot/nstar/dpnstarfit.x
new file mode 100644
index 00000000..8b864c86
--- /dev/null
+++ b/noao/digiphot/daophot/nstar/dpnstarfit.x
@@ -0,0 +1,1383 @@
+include <mach.h>
+include	<imhdr.h>
+include "../lib/daophotdef.h"
+include "../lib/apseldef.h"
+include "../lib/nstardef.h"
+
+
+# DP_NSTARFIT -- Fit the stellar group.
+
+int procedure dp_nstarfit (dao, im, nin_group, mean_sky, cdimen, iter,
+	converge)
+
+pointer	dao			# pointer to the daophot structure
+pointer im			# pointer to the input image
+int	nin_group		# original group size
+real	mean_sky		# the mean sky for the group
+int	cdimen			# dimensions of the coefficient matrix
+int	iter			# the current iteration
+bool	converge		# did the fit converge ?
+
+bool 	clip, refit
+int	i, j, k, xpix, ypix, group_size, nterm, ncols, mindex, ifaint, tifaint
+int	ntmin, ixmin, ixmax, iymin, iymax, ifxmin, ifxmax, ifymin, ifymax, flag
+pointer	apsel, psffit, nstar, subim, ypixel, pixel
+real	fitradsq, cutoff, psfradsq, sepcrit, sepmin, perr, peakerr, sky_value
+real	read_noise, mingdata, maxgdata, chigrp, wcrit, xmin, xmax, ymin, ymax
+real	datum, sumres, grpwt, xtemp, ytemp, weight, ds, pred_pixval, sigmasq
+real	relerr, dswt, faint, tfaint, noise
+
+bool	dp_nstmerge(), dp_ntomit(), dp_ntmin(), dp_ncheckc()
+pointer imgs2r()
+real	dp_ntskyval(), dp_ntsubtract()
+
+begin
+	# Define the daophot pointers.
+	psffit = DP_PSFFIT (dao)
+	apsel = DP_APSEL(dao)
+	nstar = DP_NSTAR (dao)
+
+	# Set up some daophot constants. At some point these will be computed
+	# when the NSTAR task is started up instead of at the beginning of
+	# each group fit. For the moment it is convenient and not too
+	# costly to compute them here. Also initialize the fit.
+
+	if (iter == 1) {
+
+	    # Compute the fitting and psf radii.
+	    fitradsq = DP_FITRAD (dao) ** 2
+	    psfradsq = DP_PSFRAD(dao) ** 2
+	    cutoff = CUT_FACTOR * fitradsq
+
+	    # Compute the merging parameters.
+	    if (IS_INDEFR(DP_MERGERAD(dao))) {
+	        sepcrit = 2.0 * (Memr[DP_PSFPARS(psffit)] ** 2 +
+	            Memr[DP_PSFPARS(psffit)+1] ** 2)
+	        sepmin = min (1.0, FRACTION_MINSEP * sepcrit)
+	    } else {
+		sepcrit = DP_MERGERAD(dao) ** 2
+	        sepmin = min (1.0, FRACTION_MINSEP * sepcrit)
+	    }
+
+	    # Compute the noise parameters.
+	    read_noise = (DP_READNOISE(dao) / DP_PHOTADU(dao)) ** 2
+	    perr = 0.01 * DP_FLATERR(dao)
+	    peakerr = 0.01 * DP_PROFERR(dao) / (Memr[DP_PSFPARS(psffit)] *
+	        Memr[DP_PSFPARS(psffit)+1])
+
+	    # Compute the minimum and maximum good pixel values.
+	    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)
+
+	    # Initialize the fit.
+	    chigrp = 1.0
+	    clip = false
+	    if (DP_RECENTER(dao) == YES) {
+	        nterm = 3 * nin_group
+		ntmin = 3
+	    } else {
+	        nterm = nin_group
+		ntmin = 1
+	    }
+	    if (DP_FITSKY(dao) == YES) {
+		nterm = nterm + 1
+		ntmin = ntmin + 1
+	    }
+	    call aclrr (Memr[DP_NXOLD(nstar)], nterm) 
+	    call amovkr (1.0, Memr[DP_NXCLAMP(nstar)], nterm)
+	}
+
+	# Start a new iteration.
+	group_size = nin_group
+	if (DP_RECENTER(dao) == YES)
+	    nterm = 3 * group_size
+	else
+	    nterm = group_size
+	if (DP_FITSKY(dao) == YES)
+	    nterm = nterm + 1
+
+	# Begin fitting the current group of stars.
+	repeat {
+
+	    # Initialize the convergence criteria.
+	    converge = false
+
+	    # Check that there is at least 1 star in the group.
+	    if (group_size < 1)
+	        return (group_size)
+
+	    # Set up the critical error for star rejection.
+	    if (iter >= NITER_MAX)
+	        wcrit = WCRIT_MAX
+	    else if (iter >= NITER_MED)
+	        wcrit = WCRIT_MED
+	    else if (iter >= NITER_MIN)
+	        wcrit = WCRIT_MIN
+	    else
+	        wcrit = MAX_REAL
+
+	    # Set the sky fitting derivative.
+	    if (DP_FITSKY(dao) == YES)
+	        Memr[DP_NX(nstar)+nterm-1] = -1.0
+
+	    # Initialize arrays.
+	    call aclrr (Memr[DP_APCHI(apsel)], group_size)
+	    call aclrr (Memr[DP_NSUMWT(nstar)], group_size)
+	    call aclrr (Memr[DP_NNUMER(nstar)], group_size)
+	    call aclrr (Memr[DP_NDENOM(nstar)], group_size)
+	    call amovki (NSTERR_OK, Memi[DP_NIER(nstar)], group_size)
+
+	    # Compute the minimum and maximum x and y values.
+	    call alimr (Memr[DP_APXCEN(apsel)], group_size, xmin, xmax) 
+	    call alimr (Memr[DP_APYCEN(apsel)], group_size, ymin, ymax) 
+
+	    # Check to see whether any two stars are within the critical
+	    # difference from each other.
+
+	    if ((group_size > 1) && dp_nstmerge (Memr[DP_APXCEN(apsel)],
+		Memr[DP_APYCEN(apsel)], Memr[DP_APMAG(apsel)],
+		Memr[DP_APERR(apsel)], group_size, sepcrit, sepmin, wcrit,
+		i, j, k)) {
+
+		# Compute the new centroid and brightness.
+		call dp_nstcen (Memr[DP_APXCEN(apsel)], Memr[DP_APYCEN(apsel)],
+		    Memr[DP_APMAG(apsel)], i, j, k)
+
+		# Print out verbose comments.
+		if (DP_VERBOSE (dao) == YES)  {
+		    call printf ("\tStar %-5d has merged with star %-5d\n")
+		        call pargi (Memi[DP_APID(apsel)+k-1])
+		        call pargi (Memi[DP_APID(apsel)+i-1])
+		}
+
+		# Recompute the mean sky.
+		if ((DP_FITSKY(dao) == NO) && (DP_GROUPSKY(dao) == YES))
+		    mean_sky = (mean_sky * group_size -
+		        Memr[DP_APMSKY(apsel)+k-1]) / (group_size - 1)
+
+		# Now remove the k-th star from the group.
+		call dp_remove (k, group_size, Memi[DP_APID(apsel)],
+		    Memr[DP_APXCEN(apsel)], Memr[DP_APYCEN(apsel)],
+		    Memr[DP_APMAG(apsel)], Memr[DP_APMSKY(apsel)],
+		    Memi[DP_NIER(nstar)], NSTERR_MERGE)
+
+		# After deleting a star, resize the matrix, release all of
+		# the clamps and back up the iteration counter.
+
+		if (DP_RECENTER(dao) == YES)
+		    nterm = 3 * group_size
+		else
+		    nterm = group_size
+		if (DP_FITSKY(dao) == YES)
+		    nterm = nterm + 1
+		clip = false
+		call aclrr (Memr[DP_NXOLD(nstar)], nterm)
+		call amovkr (1.0, Memr[DP_NXCLAMP(nstar)], nterm)
+		iter = max (1, iter - 1)
+		next
+	    }
+
+	    # Now we can proceed with the iteration. If this is the first
+	    # iteration read in the subraster. Determine the size of the
+	    # subraster we need to extract from the image for this group.
+
+	    if (iter == 1) {
+	        ixmin = max (1, int (xmin - DP_PSFRAD(dao) -
+		    DP_FITRAD(dao)) + 1)
+	        iymin = max (1, int (ymin - DP_PSFRAD(dao) -
+		    DP_FITRAD(dao)) + 1)
+	        ixmax = min (IM_LEN(im,1), int (xmax + DP_PSFRAD(dao) +
+		    DP_FITRAD(dao)))
+	        iymax = min (IM_LEN(im,2), int (ymax + DP_PSFRAD(dao) +
+		    DP_FITRAD(dao)))
+	        subim = imgs2r (im, ixmin, ixmax, iymin, iymax)
+	        ncols = ixmax - ixmin + 1
+	    }
+
+	    # Compute the area on the subraster that is off interest to
+	    # the current iteration.
+
+	    ifxmin = max (ixmin, int (xmin - DP_FITRAD(dao)) + 1)
+	    ifymin = max (iymin, int (ymin - DP_FITRAD(dao)) + 1)
+	    ifxmax = min (ixmax, int (xmax + DP_FITRAD(dao)))
+	    ifymax = min (iymax, int (ymax + DP_FITRAD(dao)))
+
+	    # Zero the normal matrix and the vector of residuals.
+
+	    call aclrr (Memr[DP_NV(nstar)], nterm)
+	    call aclrr (Memr[DP_NC(nstar)], cdimen * cdimen)
+	    call aclri (Memi[DP_NNPIX(nstar)], group_size)
+
+	    sumres = 0.0
+	    grpwt = 0.0
+
+	    # Loop over the pixels in the subraster.
+
+	    ypixel = subim + (ifymin - iymin) * ncols + ifxmin - ixmin - 1
+	    do ypix = ifymin, ifymax {
+
+		ytemp = ypix
+	        pixel = ypixel
+
+	        do xpix = ifxmin, ifxmax {
+
+		    xtemp = xpix
+		    pixel = pixel + 1
+		    datum = Memr[pixel]
+
+		    # Reject pixel if not in good data range.
+		    if ((datum < mingdata) || (datum > maxgdata))
+		        next
+
+		    # If this pixel is within one fitting radius of at
+		    # least one star then include it in the solution.
+		    # While figuring this out, compute the squared distance
+		    # of this pixel from the centroid of each star in the
+		    # group, and sum its contribution to the number
+		    # of pixels within one fitting radius of each star.
+
+		    if (dp_ntomit (Memr[DP_APXCEN(apsel)],
+		        Memr[DP_APYCEN(apsel)], Memr[DP_NRPIXSQ(nstar)],
+			Memi[DP_NSKIP(nstar)], group_size, xtemp, ytemp,
+			cutoff))
+		        next
+
+		    if ((DP_FITSKY(dao) == NO) && (DP_GROUPSKY(dao) == NO)) {
+			sky_value = dp_ntskyval (Memr[DP_APMSKY(apsel)],
+			    Memi[DP_NSKIP(nstar)], group_size)
+			if (IS_INDEFR(sky_value))
+			    sky_value = mean_sky
+		    } else
+			sky_value = mean_sky
+
+		    # Subtract the mean sky from the pixel.
+		    #ds = datum - mean_sky
+		    ds = datum - sky_value
+
+		    # Now loop over the stars and subtract from this pixel
+		    # the light contribution from each star within one psf
+		    # radius.
+		    #
+		    # The condition equation for pixel[xpix,ypix] is
+		    #
+		    # residual = data[xpix,ypix] - sky - sum (scale *
+		    #            psf[xpix-xc,ypix-yc])
+		    #
+		    # The scale, xc, and yc are iterated until
+		    #
+		    #          sum (weight * residual ** 2)
+		    #
+		    # is minimized. Weight will be a function of 1) the
+		    # distance of the pixel from the center of the nearest
+		    # star, 2) the model-predicted brightness of the pixel
+		    # (taking into consideration the readout noise, the
+		    # photons/ADU, and the interpolating error of the PSF,
+		    # and, 3) the size of the residual itself. One is
+		    # necessary to prevent the non-linear least squares
+		    # solution from oscillating: oftimes it will come
+		    # to pass that if you include a pixel in the solution
+		    # then the predicted shift of the centroid will cause
+		    # that pixel to be excluded in the next iteration, and the 
+		    # new predicted shift of the centroid will cause that
+		    # pixel to be included again. This could go on ad
+		    # infinitum. The cure is to have the weight of the pixel
+		    # go asymptotically to zero as its distance from the
+		    # stellar centroid approaches the fitting radius. In
+		    # the case like that just described, the solution can
+		    # then find a real minimum of the sum of the weighted
+		    # squared residuals with that pixel at some low-weight
+		    # position just inside the fitting radius. Two is
+		    # just sensible weighting. Three is a crude attempt
+		    # at making the solution more robust against bad pixels.
+
+		    weight = dp_ntsubtract (dao, im, Memr[DP_APXCEN(apsel)],
+		        Memr[DP_APYCEN(apsel)], Memr[DP_NRPIXSQ(nstar)],
+		        Memr[DP_APMAG(apsel)], Memi[DP_NSKIP(nstar)],
+		        Memr[DP_NX(nstar)], group_size, xtemp, ytemp, ds,
+			psfradsq, fitradsq, DP_RECENTER(dao))
+
+		    # At this point the vector X contains the first
+		    # derivative of the condition equation, for the pixel
+		    # under consideration, with respect to each of the
+		    # fitting parameters for all of these stars. We now need
+		    # to add these values into the normal matrix and the vector
+		    # of residuals.
+
+		    # 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
+		    # of the total brightness for the difficulty of flat-
+		    # fielding and bias subtraction, plus an estimated error
+		    # of the same fraction of the fourth derivative at the
+		    # peak of the profile, to account for the difficulty 
+		    # of accurately interpolating within the point-spread
+		    # function. The fourth derivative of the PSF is
+		    # proportional to H / sigma ** 4 (sigma is the Gaussian
+		    # width parameter for the stellar core); using the geometric
+		    # mean of sigma(x) and sigma(y), this becomes H / [sigma(x)
+		    # * sigma(y)] ** 2. The ratio of the fitting error to this
+		    # quantity is estimated to be approximately 0.027 from a
+		    # good-seeing CTIO frame. (This is probably a function of
+		    # seeing, sampling etc.)
+
+		    # Get the residual from the PSF fit and the pixel
+		    # intensity as predicted by the fit. Pred_pixval = raw data
+		    # minus residual = model predicted value of the intensity at
+		    # this point.
+
+		    pred_pixval = max (0.0, datum - ds)
+		    if ((pred_pixval > maxgdata) && (iter >= 4))
+			next
+
+		    #sigmasq = pred_pixval / DP_PHOTADU (dao) + read_noise + 
+		        #(perr * pred_pixval) ** 2 + (peakerr *
+		        #(pred_pixval - mean_sky)) ** 2
+		    sigmasq = pred_pixval / DP_PHOTADU (dao) + read_noise + 
+		        (perr * pred_pixval) ** 2 + (peakerr *
+		        (pred_pixval - sky_value)) ** 2
+		    if (sigmasq <= 0.0)
+			next
+		    relerr = abs (ds) / sqrt (sigmasq)
+
+		    # Add this residual into the weighted sum of the
+		    # absolute relative residuals.
+
+		    sumres = sumres + relerr * weight
+		    grpwt = grpwt + weight
+
+		    # Add into the accumulating sums of the weighted 
+		    # absolute relative residuals and of the image sharpness
+		    # parameter for each of the stars. Include in the
+		    # sharpness index only those pixels within NCORE_SIGMA
+		    # sigma of the centroid of the object. This saves time
+		    # and floating underflows by excluding pixels
+		    # which contribute less than about one part in a
+		    # million to the index.
+
+		    call dp_acsharp (Memr[DP_APXCEN(apsel)],
+		        Memr[DP_APYCEN(apsel)], Memi[DP_NSKIP(nstar)],
+			Memi[DP_NNPIX(nstar)], Memr[DP_NNUMER(nstar)],
+			Memr[DP_NDENOM(nstar)], Memr[DP_NSUMWT(nstar)],
+			Memr[DP_APCHI(apsel)], group_size, xtemp, ytemp,
+			Memr[DP_PSFPARS(psffit)], Memr[DP_PSFPARS(psffit)+1],
+			ds, sigmasq, relerr, weight)
+
+		    # If clipping is in effect, reduce the weight of a bad
+		    # pixel. A pixel having a residual of 2.5 sigma gets
+		    # reduced to half weight and one with a rersidual of 5
+		    # sigma gets weight of 1 / 257.
+
+		    weight = weight / sigmasq
+		    if (clip) 
+		        weight = weight / (1.0 + (relerr / (DP_CLIPRANGE(dao) *
+			    chigrp)) ** DP_CLIPEXP(dao))
+		    dswt = ds * weight
+
+		    # Work this pixel into the normal matrix.
+		    call dp_mataccum (Memr[DP_NX(nstar)], Memi[DP_NSKIP(nstar)],
+		        group_size, Memr[DP_NC(nstar)], Memr[DP_NV(nstar)],
+			cdimen, nterm, weight, dswt, DP_RECENTER(dao),
+			DP_FITSKY(dao))
+
+	        }
+
+	        ypixel = ypixel + ncols
+	    }
+
+	    # Make sure that every star in the group has at least MIN_FIT_PIXEL
+	    # pixels within one fitting radius.
+
+	    refit = dp_ntmin (Memi[DP_APID(apsel)], Memr[DP_APXCEN(apsel)],
+	        Memr[DP_APYCEN(apsel)], Memr[DP_APMAG(apsel)],
+	        Memr[DP_APMSKY(apsel)], Memi[DP_NNPIX(nstar)],
+		Memi[DP_NIER(nstar)], group_size, nterm, DP_RECENTER(dao),
+		DP_FITSKY(dao), DP_GROUPSKY(dao), mean_sky, DP_VERBOSE(dao))
+	    if (group_size < 1)
+	        return (group_size)
+	    if (refit)
+	        next
+
+	    # Reflect the normal matrix across the diagonal.
+	    call dp_mreflect (Memr[DP_NC(nstar)], cdimen, nterm)
+
+	    # Compute the robust estimate of the standard deviation of the
+	    # residuals for the group as a whole, and for each star. This
+	    # estimate is sqrt (PI/2) * weighted mean absolute relative
+	    # residual. Do you like that "absolute relative" stuff? (PBS)
+	    # NO! (LED)
+	    #
+	    #      CHI = CHI_NORM * SUM (weight * resid) / (# of pixels)
+	    # 
+	    # This gets corrected for bias by being multiplied by
+	    #
+	    #	SQRT (# of pixels) / (# of pixels - 3)
+
+	    # Then the value is driven towards unity, depending on
+	    # exactly how many pixels were involved: if CHIOLD is based
+	    # on a total weight of 3, then it is extremely poorly determined
+	    # and we just want to keep CHIOLD = 1. The larger GRPWT is, the
+	    # better determined CHIOLD is, and the less we want to force it
+	    # toward unity. So, just take the weighted average of CHIOLD and 
+	    # unity, with weights GRPWT - 3 and 1, respectively.
+
+	    if (grpwt > ntmin) {
+	        chigrp = CHI_NORM * sumres * sqrt (1.0 / (grpwt * (grpwt -
+	            ntmin)))
+	        chigrp = ((grpwt - ntmin) * chigrp + ntmin) / grpwt
+	    }
+
+	    # CHIOLD has been pulled toward 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 
+	    # greated than one, CHI will be nudged toward the group average.
+	    # In order to work optimally, of course, this requires that	
+	    # the # of photons / ADU, the READNOISE and the other noise
+	    # contributors are properly specified.
+
+	    call dp_ntchi (Memr[DP_NSUMWT(nstar)], Memr[DP_APCHI(apsel)],
+	        group_size, ntmin, chigrp, grpwt)
+
+	    # Invert the matrix.
+	    call invers (Memr[DP_NC(nstar)], cdimen, nterm, flag)
+
+	    # Check for a singular matrix.
+	    refit = dp_ncheckc (Memr[DP_NC(nstar)], cdimen, nterm,
+		Memi[DP_APID(apsel)], Memr[DP_APXCEN(apsel)],
+		Memr[DP_APYCEN(apsel)], Memr[DP_APMAG(apsel)],
+		Memr[DP_APMSKY(apsel)], Memi[DP_NIER(nstar)], group_size,
+		DP_RECENTER(dao), DP_FITSKY(dao), DP_GROUPSKY(dao), mean_sky,
+		DP_VERBOSE(dao)) 
+	    if (group_size < 1)
+	        return (group_size)
+	    if (refit)
+	        next
+
+	    # Solve for position and scale factor increments.
+	    call mvmul (Memr[DP_NC(nstar)], cdimen, nterm, Memr[DP_NV(nstar)],
+	        Memr[DP_NX(nstar)])
+
+	    if (iter <= 1)
+	        refit = true
+	    else
+	        refit = false
+
+	    # Fit the sky.
+	    if (DP_FITSKY(dao) == YES) {
+	        noise = sqrt (abs (mean_sky / DP_PHOTADU(dao)) + read_noise)
+	        mean_sky = mean_sky - max (-3.0 * noise,
+		    min (Memr[DP_NX(nstar)+nterm-1], 3.0 * noise))
+		if (abs (Memr[DP_NX(nstar)+nterm-1]) > (1.0e-4 * mean_sky))
+		    refit = true
+	    }
+
+	    # 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 pixels per iteration. Any time that the 
+	    # parameter correction changes sign from one iteration to the
+	    # next, the maximum permissible change will be reduced by a factor
+	    # of two. These clamps are released any time a star disappears.
+
+	    call dp_ntclamp (Memr[DP_APXCEN(apsel)], Memr[DP_APYCEN(apsel)],
+	        Memr[DP_APMAG(apsel)], Memr[DP_APERR(apsel)],
+	        Memr[DP_NSUMWT(nstar)], group_size, Memr[DP_NC(nstar)], cdimen,
+		Memr[DP_NX(nstar)], Memr[DP_NXOLD(nstar)],
+		Memr[DP_NXCLAMP(nstar)], DP_RECENTER(dao), clip, refit)
+
+	    # Check whether the estimated centroid of the any star has
+	    # moved so far out of the limits of the picture that it has fewer
+	    # than 4 or 5 pixels within one fitting radius.
+
+	    call dp_ntcentroid (Memi[DP_APID(apsel)], Memr[DP_APXCEN(apsel)],
+	        Memr[DP_APYCEN(apsel)], Memr[DP_APMAG(apsel)],
+	        Memr[DP_APMSKY(apsel)], Memi[DP_NIER(nstar)], group_size,
+		ixmin, ixmax, iymin, iymax, fitradsq, DP_FITSKY(dao),
+		DP_GROUPSKY(dao), mean_sky, refit, DP_VERBOSE(dao))
+
+	    if (group_size < 1)
+	        return (group_size)
+
+	    # Update matrix dimensions.
+	    if (DP_RECENTER(dao) == YES)
+	        nterm = 3 * group_size
+	    else
+		nterm = group_size
+	    if (DP_FITSKY(dao) == YES)
+		nterm = nterm + 1
+
+	    # Now check whether any of the stars is too faint (more than 12.5
+	    # magnitudes fainter than the PSF star). If several stars are too
+	    # faint, delete the faintest one, and set the brightness of the 
+	    # other faint ones to 12.5 magnitudes below the PSF star. That way
+	    # on the next iteration we will see whether these stars want to 
+	    # grow or to disappear.
+
+	    faint = 0.0
+	    ifaint = 0
+	    call dp_ntfmag (Memr[DP_APMAG(apsel)], group_size, tfaint, tifaint)
+
+	    # If at least one star is more than 12.5 magnitudes fainter
+	    # than the PSF then ifaint is the relative index of the faintest
+	    # of them, and faint is the relative brightness of the
+	    # faintest of them.
+
+	    # No very faint star was detected.
+	    if (tifaint <= 0) {
+
+	        # If the solution has not converged and if the number of
+	        # iterations is less than MIN_ITER, perform another iteration
+	        # with no questions asked.
+
+	        if ((refit) && (iter < MIN_ITER))
+		    return (group_size)
+	    
+	        # If the solution doesn't think it has converged, after the
+	        # fourth iteration delete the least certain star if it is less
+	        # less than a one-sigma detection; after the eighth iteration
+	        # delete the least certain star if it is less than a 1.5 sigma
+	        # detection; after the twelfth iteration OR if the solution
+	        # thinks it has converged, delete the least certain star if it 
+	        # is less than a two-sigma detection.
+
+	        call dp_fsnoise (Memr[DP_APMAG(apsel)], Memr[DP_APERR(apsel)],
+		    group_size, faint, ifaint)
+
+	        if ((refit) && (iter < DP_MAXITER (dao)) && (faint < wcrit))
+		    return (group_size)
+	    }
+
+	    # Reject the appropriate star.
+	    if ((tifaint > 0) || (faint >= MIN_FAINT)) {
+
+	        # Either (a) the solution thinks it has not converged
+	        # and the faintest star is more uncertain than sqrt(wcrit)
+	        # or (b) the solution thinks it has converged and the
+	        # faintest star is more uncertain than two-sigma.
+
+	        if (DP_VERBOSE (dao) == YES) {
+		    mindex = max (tifaint, ifaint)
+		    call printf (
+		    "\tStar %-5d has been deleted because it is too faint\n")
+		        call pargi (Memi[DP_APID(apsel)+mindex-1])
+	        }
+
+		if ((DP_FITSKY(dao) == NO) && (DP_GROUPSKY(dao) == YES) &&
+		    (group_size > 1))
+		    mean_sky = (mean_sky * group_size -
+		        Memr[DP_APMSKY(apsel)+max(tifaint,ifaint)-1]) /
+			    (group_size - 1)
+
+	        call dp_remove (max (tifaint, ifaint), group_size,
+	            Memi[DP_APID(apsel)], Memr[DP_APXCEN(apsel)],
+		    Memr[DP_APYCEN(apsel)], Memr[DP_APMAG(apsel)],
+		    Memr[DP_APMSKY(apsel)], Memi[DP_NIER(nstar)],
+		    NSTERR_FAINT)
+
+	        if (group_size < 1)
+		    return (group_size)
+
+		if (DP_RECENTER(dao) == YES)
+	            nterm = 3 * group_size
+		else
+		    nterm = group_size
+		if (DP_FITSKY(dao) == YES)
+		    nterm = nterm + 1
+	        call aclrr (Memr[DP_NXOLD(nstar)], nterm)
+	        call amovkr (1.0, Memr[DP_NXCLAMP(nstar)], nterm)
+	        clip = false
+	        iter = max (1, iter - 1)
+	        next
+	    }
+
+	    # Solution has either converged or gone to the maximum number
+	    # of iterations.
+
+	    if ((iter < DP_MAXITER(dao)) && (! clip)) {
+
+	        # The first convergence milestone has been reached. Turn on the
+	        # clipper, loosen the clamps and keep on going.
+
+		if (DP_CLIPEXP(dao) > 0)
+	            clip = true
+	        converge = false
+	        call aclrr (Memr[DP_NXOLD(nstar)], nterm)
+	        call amaxkr (Memr[DP_NXCLAMP(nstar)], 0.25,
+	            Memr[DP_NXCLAMP(nstar)], nterm)
+	        return (group_size)
+	    }
+
+	    converge = true		    
+
+	} until (converge)
+
+	return (group_size)
+end
+
+
+# DP_NSTMERGE -- Decide whether two stars in a group should merge.
+
+bool procedure dp_nstmerge (xcen, ycen, mag, magerr, group_size, sepcrit,
+	sepmin, wcrit, i, j, k)
+
+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
+int	group_size		# group size
+real	sepcrit			# the critical separation squared
+real	sepmin			# the minimum separation
+real	wcrit			# critical error for rejection
+int	i, j, k			# output indices
+
+real	separation
+
+begin
+	do i = 1, group_size {
+  	    do j = 1,  i - 1 {
+
+		# Compute the separation.
+	        separation = (xcen[j] - xcen[i]) ** 2 +
+		    (ycen[j] - ycen[i]) ** 2
+		if (separation > sepcrit)
+		    next
+
+		# Find the fainter of the two stars.
+		k = j
+		if (mag[i] < mag[j]) 
+		    k = i
+		if ((separation <  sepmin) || ((magerr[k] / mag[k])  > wcrit))
+		        return (true)
+	    }
+	}
+
+	return (false)
+end
+
+
+# DP_NSTCEN --  Recompute the centroid and brightness of the i-th star.
+
+procedure dp_nstcen (xcen, ycen, mag, i, j, k)
+
+real	xcen[ARB]		# the x centers
+real	ycen[ARB]		# the y centers
+real	mag[ARB]		# the magnitudes
+int	i, j, k			# array indices
+
+begin
+	# Now eliminate the fainter of the two stars. The k-th
+	# star is now the fainter of the two, the i-th the
+	# brighter.
+
+        if (mag[i] < mag[j])
+            i = j
+
+	# Now replace the centroid of the i-th star with the
+	# weighted mean of the most recent estimates of the 
+	# centroids of the i-th and the k-th stars, and the
+	# brightness of i-th with the sum of the brightnesses.
+
+	xcen[i] = xcen[i] * mag[i] + xcen[k] * mag[k]
+	ycen[i] = ycen[i] * mag[i] + ycen[k] * mag[k]
+	mag[i] = mag[i] + mag[k]
+	xcen[i] = xcen[i] / mag[i]
+	ycen[i] = ycen[i] / mag[i]
+end
+
+
+# DP_NTOMIT -- Check whether a pixel is within one fitting radius of another
+# star. 
+
+bool procedure dp_ntomit (xcen, ycen, rpixsq, skip, group_size, fx, fy, cutoff)
+
+real	xcen[ARB]		# array of x centers
+real	ycen[ARB]		# array of y centers
+real	rpixsq[ARB]		# array of distances squared
+int	skip[ARB]		# array of skip values
+int	group_size		# the group size
+real	fx, fy			# pixel position in image
+real	cutoff			# radius cuttoff for pixel inclusion
+
+bool	omit
+int	i
+
+begin
+	omit = true
+	do i = 1, group_size {
+	    skip[i] = YES
+	    rpixsq[i] = (fx - xcen[i]) ** 2 + (fy - ycen[i]) ** 2
+	    if (rpixsq[i] > cutoff)
+		next
+	    skip[i] = NO
+	    omit = false
+	}	 
+
+	return (omit)
+end
+
+
+# DP_NTSKYVAL -- Compute the average sky value to use for a particular
+# pixel by averaging the sky values of all stars for which the
+# pixel is within one fitting radius.
+
+real procedure dp_ntskyval (sky, skip, nstar)
+
+real	sky[ARB]		# array of sky values
+int	skip[ARB]		# array of skip values
+int	nstar			# the number of stars
+
+int	i, npts
+real	sum
+
+begin
+	sum = 0.0
+	npts = 0
+	do i = 1, nstar {
+	    if (skip[i] == YES)
+		next
+	    sum = sum + sky[i]
+	    npts = npts + 1
+	}
+	if (npts <= 0)
+	    return (INDEFR)
+	else
+	    return (sum / npts)
+end
+
+
+# DP_NTSUBTRACT -- Procedure to subtract the contribution of a particular
+# pixel from a particular star.
+
+real procedure dp_ntsubtract (dao, im, xcen, ycen, rpixsq, mag, skip, x,
+	group_size, fx, fy, ds, psfradsq, fitradsq, recenter)
+
+pointer	dao			# pointer to the daophot structure
+pointer	im			# the input image descriptor
+real	xcen[ARB]		# array of x centers
+real	ycen[ARB]		# array of y centers
+real	rpixsq[ARB]		# array of distances squared
+real	mag[ARB]		# array of magnitudes
+int	skip[ARB]		# array of skip values
+real	x[ARB]			# x accumulation vector array
+int	group_size		# size of the group
+real	fx, fy			# center of pixel in image
+real	ds			# pixel value
+real	psfradsq		# psf radius squared
+real	fitradsq		# fit radius squared
+int	recenter		# recenter the coordinates
+
+int	i, i3, k
+pointer	psffit
+real	weight, dx, dy, deltax, deltay, val, dvdx, dvdy, rsq
+real	dp_usepsf()
+
+begin
+	psffit = DP_PSFFIT(dao)
+
+	weight = 0.0
+	do i = 1, group_size {
+	    if (rpixsq[i] >= psfradsq)
+		next
+	    dx = fx - xcen[i]
+	    dy = fy - 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
+	    val = 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)
+	    ds = ds - mag[i] * val
+	    if (skip[i] == YES)
+		next
+	    rsq = rpixsq[i] / fitradsq
+	    weight = max (weight, 5.0 / (5.0 + rsq / (1.0 - rsq)))
+	    if (recenter == YES) {
+	        i3 = 3 * i
+	        k = i3 - 2
+	        x[k] = -val
+	        k = i3 - 1
+	        x[k] = -mag[i] * dvdx
+	        x[i3] = -mag[i] * dvdy
+	    } else
+		x[i] = -val
+	}
+
+	return (weight)
+end
+
+
+# DP_ACSHARP -- Procedure to accumulate sums of the weighted absolute
+# relative residuals and the image sharpness parameter for each of the
+# stars.
+
+procedure dp_acsharp (xcen, ycen, skip, npix, numer, denom, sumwt, chi,
+	group_size, fx, fy, fwhmx, fwhmy, ds, sigmasq, relerr, weight)
+
+real	xcen[ARB]		# array of object x centers
+real	ycen[ARB]		# array of object y centers
+int	skip[ARB]		# array of skip values
+int	npix[ARB]		# array of numbers of pixels
+real	numer[ARB]		# numerator array
+real	denom[ARB]		# denominator array
+real	sumwt[ARB]		# array of summed weights
+real	chi[ARB]		# array of chis
+int	group_size		# group size paramter.
+real	fx, fy			# position of data in image
+real	fwhmx, fwhmy		# gaussian core widths in x and y
+real	ds			# the data residual
+real	sigmasq			# the sigma squared
+real	relerr			# the relative error
+real	weight			# the weight
+
+int	i
+real	rhosq, dfdsig
+
+begin
+	do i = 1, group_size {
+
+	    if (skip[i] == YES)
+		next
+
+	    # Accumulate the number of pixels, chi and sum of the weights.
+	    npix[i] = npix[i] + 1
+	    chi[i] = chi[i] + relerr * weight
+	    sumwt[i] = sumwt[i] + weight
+
+	    # Include in the sharpness index only those pixels
+	    # within NCORE_SIGMASQ of the centroid of the
+	    # object. (This saves time and floating underflows
+	    # by excluding pixels which contribute very little
+	    # to the index.
+
+	    rhosq = ((xcen[i] - fx) / fwhmx) ** 2 + ((ycen[i] - fy) /
+	        fwhmy) ** 2
+	    if (rhosq > NCORE_SIGMASQ)
+		next
+	    rhosq = 0.6931472 * rhosq
+	    dfdsig = exp (-rhosq) * (rhosq - 1.0)
+	    numer[i] = numer[i] + dfdsig * ds / sigmasq
+	    denom[i] = denom[i] + (dfdsig ** 2) / sigmasq
+	}			
+end
+
+
+# DP_MATACCUM -- Procedure to accumulate the data into the matrices.
+
+procedure dp_mataccum (x, skip, group_size, c, v, cdimen, nterm, weight, dswt,
+	recenter, fitsky)
+
+real	x[ARB]			# x array
+int	skip[ARB]		# skip vector
+int	group_size		# size of the group
+real	c[cdimen,ARB]		# coefficient matrix
+real	v[ARB]			# vector array
+int	cdimen			# dimensions of the coefficient matrix
+int	nterm			# the number of terms
+real	weight			# weight
+real	dswt			# data weight
+int	recenter		# recenter the coordinates
+int	fitsky			# fit the sky value
+
+int	i, i3, i3m2, k, j, l
+
+begin
+	if (fitsky == YES) {
+	    c[nterm,nterm] = c[nterm,nterm] + weight
+	    v[nterm] = v[nterm] - dswt
+	}
+
+	do i = 1, group_size {
+	    if (skip[i] == YES)
+		next
+	    if (recenter == YES) {
+	        i3 = i * 3
+	        i3m2 = i3 - 2
+	        do k = i3m2, i3 {
+		    if (fitsky == YES)
+			c[nterm,k] = c[nterm,k] - x[k] * weight
+		    v[k] = v[k] + x[k] * dswt
+		}
+	        do j = 1, i {
+		    if (skip[j] == YES)
+		        next
+		    do k = i3m2, i3 {
+		        do l = 3 * j - 2, min (k, 3 * j)
+			    c[k,l] = c[k,l] + x[k] * x[l] * weight
+		    }
+	        }
+	    } else {
+		v[i] = v[i] + x[i] * dswt
+		if (fitsky == YES)
+			c[nterm,i] = c[nterm,i] - x[i] * weight
+		do j = 1, i {
+		    if (skip[j] == YES)
+		        next
+		    c[i,j] = c[i,j] + x[i] * x[j] * weight
+		}
+	    }
+	}
+end
+
+
+# DP_NTMIN -- Make sure that every star in the group has at least 
+# MIN_NPIX pixels within one fitting radius.
+
+bool procedure dp_ntmin (ids, xcen, ycen, mag, sky, npix, nier, group_size,
+	nterm, recenter, fitsky, groupsky, mean_sky, verbose)
+
+int	ids[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	sky[ARB]		# array of sky values
+int	npix[ARB]		# array of pixel numbers
+int	nier[ARB]		# array of error codes
+int	group_size		# size of the group
+int	nterm			# number of terms
+int	recenter		# recenter the objects
+int	fitsky			# fit the sky value
+int	groupsky		# use group sky value
+real	mean_sky		# the current mean sky value
+int	verbose			# verbose flag
+
+bool	redo
+int	i
+
+begin
+	redo = false
+	do i = 1, group_size {
+	    if (npix[i] >= MIN_NPIX)
+	        next
+	    redo = true
+	    if (verbose == YES) {
+		call printf (
+		    "\tStar %-5d has been deleted: too few good pixels\n")
+		    call pargi (ids[i])
+	    }
+	    if ((fitsky == NO) && (groupsky == YES) && (group_size > 1))
+		mean_sky = (mean_sky * group_size - sky[i]) / (group_size - 1)
+	    call dp_remove (i, group_size, ids, xcen, ycen, mag, sky,
+		nier, NSTERR_NOPIX)
+	    if (group_size <= 0)
+		return (redo)
+	    if (recenter == YES)
+	        nterm = 3 * group_size
+	    else
+		nterm = group_size
+	    if (fitsky == YES)
+		nterm = nterm + 1
+	}
+
+	return (redo)
+end
+
+
+# DP_MREFLECT -- Reflect the normal matrix around the diagonal.
+
+procedure dp_mreflect (c, cdimen, nterm)
+
+real	c[cdimen,ARB]		# coefficient matrix
+int	cdimen			# dimension of the c matrix
+int	nterm			# number of terms
+
+int	l, k
+
+begin
+	# Reflect the normal matrix across the diagonal.
+	do l = 2, nterm {
+	    do k = 1, l - 1
+		c[k,l] = c[l,k]
+	}
+end
+
+
+# DP_NTCHI -- Compute the chi value for each star.
+
+procedure dp_ntchi (sumwt, chi, group_size, ntmin, chigrp, grpwt)
+
+real	sumwt[ARB]		# sum of the weights
+real	chi[ARB]		# the chis:wq
+int	group_size		# size of the group
+int	ntmin			# minimum number of points
+real	chigrp			# the group chi
+real	grpwt			# the group weight
+
+int	i
+
+begin
+	do i = 1, group_size {
+	    if (sumwt[i] > ntmin) {
+		chi[i] = CHI_NORM * chi[i] / sqrt ((sumwt[i] - ntmin) *
+		    sumwt[i])
+		sumwt[i] = ((sumwt[i] - ntmin) * chi[i] + MIN_SUMWT) /
+		    sumwt[i]
+	    } else { 
+		chi[i] = chigrp
+		sumwt[i] = grpwt
+	    }
+	}
+end
+
+
+# DP_NCHECKC -- Check the inverted matrix for singularity.
+
+bool procedure dp_ncheckc (c, cdimen, nterm, ids, xcen, ycen, mag, sky,
+	nier, group_size, recenter, fitsky, groupsky, mean_sky, verbose)
+
+real	c[cdimen,ARB]		# coefficient matrix
+int	cdimen			# dimension of the c matrix
+int	nterm			# number of terms
+int	ids[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	sky[ARB]		# array of sky values
+int	nier[ARB]		# array of error codes
+int	group_size		# size of the group
+int	recenter		# recenter the objects
+int	fitsky			# fit the sky value
+int	groupsky		# use group sky value
+real	mean_sky		# the current mean sky value
+int	verbose			# verbose flag
+
+bool	redo
+int	j, starno, i
+
+begin
+	redo = false
+	do j = 1, nterm {
+	    if (c[j,j] > 0.0)
+		next
+	    redo = true
+	    if ((j == nterm) && (fitsky == YES))
+		starno = 0
+	    else if (recenter == YES)
+		starno = (j + 2) / 3
+	    else
+		starno = j
+	    if (starno == 0) {
+		if (verbose == YES) {
+		    do i = 1, group_size {
+		        call printf (
+			    "\tStar %-5d has been deleted: singular matrix\n")
+			    call pargi (ids[i])
+		    }
+		}
+		call amovki (NSTERR_SINGULAR, nier, group_size)
+		group_size = 0
+	    } else {
+		if (verbose == YES) {
+		    call printf (
+			"\tStar %-5d has been deleted: singular matrix\n")
+			call pargi (ids[starno])
+		}
+		if ((fitsky == NO) && (groupsky == YES) && (group_size > 1))
+		    mean_sky = (mean_sky * group_size - sky[starno]) /
+			(group_size - 1)
+		call dp_remove (starno, group_size, ids, xcen, ycen, mag,
+		    sky, nier, NSTERR_SINGULAR)
+	    }
+	    if (group_size <= 0)
+		return (redo)
+	    if (recenter == YES)
+		nterm = 3 * group_size
+	    else
+		nterm = group_size
+	    if (fitsky == YES)
+		nterm = nterm + 1
+	}
+
+	return (redo)
+end
+
+
+# DP_NTCLAMP -- Restrict the amount the solution can vary on each iteration.
+
+procedure dp_ntclamp (xcen, ycen, mag, magerr, sumwt, group_size, c, cdimen,
+	x, xold, clamp, recenter, clip, redo)
+
+real	xcen[ARB]		# x centers array
+real	ycen[ARB]		# y centers array
+real	mag[ARB]		# magnitude array
+real	magerr[ARB]		# magnitude errors array
+real	sumwt[ARB]		# array of weight sums
+int	group_size		# size of the group
+real	c[cdimen, ARB]		# coefficient matrix
+int	cdimen			# dimensions of c
+real	x[ARB]			# x vector
+real	xold[ARB]		# old x vector
+real	clamp[ARB]		# clamp on solution matrix
+int	recenter		# recenter the objects
+bool	clip			# clip the matrix
+bool	redo			# redo the solution
+
+int	i, l, j, k
+real	df
+
+begin
+	do i = 1, group_size {
+
+
+	    # If any correction has changed sign since the last
+	    # iteration, reduce the maximum permissible change by
+	    # a factor of two.
+
+	    # 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 (the least-squares thinks that the star should
+	    # be brighter) a change of 1 magnitude is a change of a
+	    # factor of 2.5; if the brightness correction is positive
+	    # (the star should be fainter) a change of 1 magnitude
+	    # is a change of 60%.
+
+	    if (recenter == YES) {
+
+	        l = 3 * i
+	        k = l - 1
+	        j = l - 2
+
+	        if ((xold[j] * x[j]) < 0.0) 
+		    clamp[j] = 0.5 * clamp[j]
+	        mag[i] = mag[i] - x[j] / (1.0 + max (x[j] /
+		    (MAX_DELTA_FAINTER * mag[i]), -x[j] / (MAX_DELTA_BRIGHTER *
+		    mag[i])) / clamp[j])
+	        xold[j] = x[j]
+
+	        if ((xold[k] * x[k]) < 0.0)
+		    clamp[k] = 0.5 * clamp[k]
+	        if ((xold[l] * x[l]) < 0.0)
+		    clamp[l] = 0.5 * clamp[l]
+	        xcen[i] = xcen[i] - x[k] / (1.0 + abs(x[k]) / (clamp[k] *
+	            MAX_DELTA_PIX))
+	        ycen[i] = ycen[i] - x[l] / (1.0 + abs(x[l]) / (clamp[l] *
+	            MAX_DELTA_PIX))
+	        xold[k] = x[k]
+	        xold[l] = x[l]
+	        magerr[i] =sumwt[i] *  sqrt (c[j,j])
+
+	    } else {
+
+	        if ((xold[i] * x[i]) < 0.0) 
+		    clamp[i] = 0.5 * clamp[i]
+	        mag[i] = mag[i] - x[i] / (1.0 + max (x[i] /
+		    (MAX_DELTA_FAINTER * mag[i]), -x[i] / (MAX_DELTA_BRIGHTER *
+		    mag[i])) / clamp[i])
+	        xold[i] = x[i]
+	        magerr[i] =sumwt[i] *  sqrt (c[i,i])
+	    }
+
+
+	    # There are two milestones in the convergence process: the fits
+	    # proceed normally until each star's magnitude changes by less
+	    # than its standard error or MAX_NEW_ERRMAG magnitudes, whichever
+	    # is greater, and its x and y centroids change by less than 0.02 
+	    # pixel. At this point the least squares begins to apply
+	    # down-weighting of pixels with large residuals as described
+	    # above.  The fits then continue until each star's
+	    # magnitude changes by less than MAX (MAX_NEW_ERRMAG * std. error,
+	    # MAX_NEW_RELBRIGHT2 magnitude), ad its centroids change by
+	    # less than 0.002 pixel.
+
+	    if (redo)
+		next
+
+	    if (clip) {
+		if (abs (x[j]) > max (MAX_NEW_ERRMAG * magerr[i],
+		        MAX_NEW_RELBRIGHT2 * mag[i])) {
+		    redo = true
+		} else if (recenter == YES) {
+		    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 {
+		if (abs (x[j]) > max (magerr[i], MAX_NEW_RELBRIGHT1 *
+		        mag[i])) {
+		    redo = true
+		} else if (recenter == YES) {
+		    df = sumwt[i] ** 2
+		    if (x[k] ** 2 > max (df * c[k,k], MAX_PIXERR1))
+			redo = true
+		    else if (x[l] ** 2 > max (df * c[l,l], MAX_PIXERR1))
+			redo = true
+		}
+	    }
+	}
+end
+
+
+# DP_NTCENTROID -- Check the new centroids to see if they have moved too
+# far off the edge of the image.
+
+procedure dp_ntcentroid (ids, xcen, ycen, mag, sky, nier, group_size, ixmin,
+	ixmax, iymin, iymax, fitradsq, fitsky, groupsky, mean_sky, redo,
+	verbose)
+
+int	ids[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	sky[ARB]			# array of sky values
+int	nier[ARB]			# array of error codes
+int	group_size			# size of the group
+int	ixmin,ixmax			# subraster x limits
+int	iymin,iymax			# subraster y limits
+real	fitradsq			# fit radius squared
+int	fitsky				# fit the sky value
+int	groupsky			# use the group sky value
+real	mean_sky			# the mean sky value
+bool	redo				# redo fit
+int	verbose				# verbose mode
+
+int	i
+real	dx, dy
+
+begin
+	# Check whether the centroid of any star has moved so far outside
+	# the picture that it has fewer than four or five pixels within
+	# one fitting radius.
+
+	do i = 1, group_size {
+
+	    # If the centroid of the star is outside the picture in x or
+	    # y, then DX or DY is its distance from the center of the edge 
+	    # pixel; otherwise DX and DY are zero.
+
+	    dx = max (ixmin - xcen[i], xcen[i] - ixmax, 0.0)
+	    dy = max (iymin - ycen[i], ycen[i] - iymax, 0.0)
+	    if ((dx <= MAX_PIX_INCREMENT) && (dy <= MAX_PIX_INCREMENT))
+		    next
+	    if (((dx + 1.0) ** 2 + (dy + 1.0) ** 2) < fitradsq)
+		next
+
+	    # Print a warning message about the star.
+	    if (verbose == YES) {
+		call printf (
+		"\tStar %-5d has been deleted: new center too far off image\n")
+		    call pargi (i)
+	    }
+
+	    # Adjust the sky.
+	    if ((fitsky == NO) && (groupsky == YES) && (group_size > 1))
+		mean_sky = (mean_sky * group_size - sky[i]) / (group_size - 1)
+
+	    # Delete it.
+	    call dp_remove (i, group_size, ids, xcen, ycen, mag, sky, nier,
+	        NSTERR_OFFIMAGE)
+	    if (group_size < 1)
+		break
+	    redo = true
+	}
+end
+
+
+# DP_NTFMAG -- Check for faint stars.
+
+procedure dp_ntfmag (mag, group_size, faint, ifaint)
+
+real	mag[ARB]		# array of magnitudes
+int	group_size		# size of the group
+real	faint			# faintest magnitude
+int	ifaint			# index of faintest magnitude
+
+int	i
+
+begin
+	faint = 1.0
+	ifaint = 0
+	do i = 1, group_size {
+	    if (mag[i] > MIN_REL_BRIGHT)
+		next
+	    if (mag[i] <= faint) {
+		faint = mag[i]
+		ifaint = i
+	    }
+	    mag[i] = MIN_REL_BRIGHT
+	}
+end
+
+
+# DP_FSNOISE -- Compute the smallest signal to noise ratio.
+
+procedure dp_fsnoise (mag, magerr, group_size, faint, ifaint)
+
+real	mag[ARB]		# array of magnitudes
+real	magerr[ARB]		# array of magnitude errors
+int	group_size		# size of group
+real	faint			# faint value
+int	ifaint			# faint index
+
+int	i
+real	weight
+
+begin
+	faint = 0.0
+	ifaint = 0
+	do i = 1, group_size {
+	    weight = magerr[i] / mag[i]
+	    if (weight < faint)
+		next
+	    faint = weight
+	    ifaint = i
+	}
+end
+
+
+# DP_REMOVE -- Remove the i-th star from the list of stars in the current
+# group by moving it to the end of the group.
+
+procedure dp_remove (i, nstar, ids, xcen, ycen, mag, sky, nier, pier)
+
+int	i			# the star to be removed
+int	nstar			# the number of stars in the group
+int	ids[ARB]		# the array of star ids
+real	xcen[ARB]		# the array of star x positions
+real	ycen[ARB]		# the array of star y positions
+real	mag[ARB]		# the array of magnitudes
+real	sky[ARB]		# the array of sky values.
+int	nier[ARB]		# array of error codes
+int	pier			# error code for deleted star
+
+int	ihold, phold
+real	xhold, yhold, shold, mhold
+
+begin
+	nier[i] = pier
+	if (i != nstar) {
+
+	    ihold = ids[i]
+	    xhold = xcen[i]
+	    yhold = ycen[i]
+	    shold = sky[i]
+	    mhold = mag[i]
+	    phold = nier[i]
+
+	    ids[i] = ids[nstar]
+	    xcen[i] = xcen[nstar]
+	    ycen[i] = ycen[nstar]
+ 	    sky[i] = sky[nstar]
+	    mag[i] = mag[nstar]
+	    nier[i] = nier[nstar]
+
+	    ids[nstar] = ihold
+	    xcen[nstar] = xhold
+	    ycen[nstar] = yhold
+	    sky[nstar] = shold
+	    mag[nstar] = mhold
+	    nier[nstar] = phold
+	}
+	nstar = nstar - 1
+end