Index: branches/eam_branches/ipp-20110710/psphot/src/Makefile.am
===================================================================
--- branches/eam_branches/ipp-20110710/psphot/src/Makefile.am	(revision 32157)
+++ branches/eam_branches/ipp-20110710/psphot/src/Makefile.am	(revision 32205)
@@ -130,4 +130,5 @@
 	psphotVisual.c                 \
 	psphotCullPeaks.c	       \
+	psphotFootprintSaddles.c       \
 	psphotVersion.c		       \
 	psphotModelGroupInit.c	       \
@@ -183,4 +184,5 @@
 	psphotRadialPlot.c	       \
 	psphotKronMasked.c             \
+	psphotKronIterate.c            \
 	psphotDeblendSatstars.c	       \
 	psphotMosaicSubimage.c	       \
Index: branches/eam_branches/ipp-20110710/psphot/src/psphot.h
===================================================================
--- branches/eam_branches/ipp-20110710/psphot/src/psphot.h	(revision 32157)
+++ branches/eam_branches/ipp-20110710/psphot/src/psphot.h	(revision 32205)
@@ -482,3 +482,9 @@
 bool psphotSetRadiusMomentsExact (float *fitRadius, float *windowRadius, pmReadout *readout, pmSource *source, psImageMaskType markVal);
 
+bool psphotFootprintSaddles(pmReadout *readout, psArray *footprints);
+bool psphotMaskFootprint (pmReadout *readout, pmSource *source, psImageMaskType markVal);
+
+bool psphotKronIterate (pmConfig *config, const pmFPAview *view, const char *filerule);
+bool psphotKronIterateReadout(pmConfig *config, psMetadata *recipe, const pmFPAview *view, pmReadout *readout, psArray *sources, pmPSF *psf);
+
 #endif
Index: branches/eam_branches/ipp-20110710/psphot/src/psphotExtendedSourceFits.c
===================================================================
--- branches/eam_branches/ipp-20110710/psphot/src/psphotExtendedSourceFits.c	(revision 32157)
+++ branches/eam_branches/ipp-20110710/psphot/src/psphotExtendedSourceFits.c	(revision 32205)
@@ -319,18 +319,8 @@
 	psphotSetRadiusMomentsExact(&fitRadius, &windowRadius, readout, source, markVal);
 
-	// Recalculate the source moments using the larger extended-source moments radius.  At
-	// this stage, skip Gaussian windowing, and do not clip pixels by S/N.  Save the
-	// psf-based moments for output
-	psfMoments = *source->moments;
-	if (!pmSourceMoments (source, windowRadius, 0.0, 0.0, 0.0, maskVal & PS_NOT_IMAGE_MASK(markVal))) {
-	    // subtract the best fit from the object, leave local sky
-	    fprintf (stderr, "skipping (2) %f, %f\n", source->peak->xf, source->peak->yf);
-	    *source->moments = psfMoments;
-
-	    // probably need to bump up the size for subtraction
-	    pmSourceSub (source, PM_MODEL_OP_FULL, maskVal);
-	    // XXX raise an error flag of some kind
-	    continue;
-	}
+	// UPDATE : we have changed the moments calculation.  There is now an iteration within 
+	// psphotKronMasked to determine moments appropriate for a larger object.  The values
+	// Mrf, KronFlux, and KronFluxErr are calculated for the iterated radius.  The other
+	// values are left at the psf-based values.
 
         // allocate the array to store the model fits
@@ -338,8 +328,10 @@
             source->modelFits = psArrayAllocEmpty (models->list->n);
         }
-	// extFitPars are the non-parametric data measured for this model fit (moments, kron, etc)
-        if (source->extFitPars == NULL) {
-            source->extFitPars = psArrayAllocEmpty (models->list->n);
-        }
+
+	// XXX deprecate?
+	// XXX // extFitPars are the non-parametric data measured for this model fit (moments, kron, etc)
+        // XXX if (source->extFitPars == NULL) {
+        // XXX     source->extFitPars = psArrayAllocEmpty (models->list->n);
+        // XXX }
 
         // loop here over the models chosen for each source (exclude by S/N)
@@ -405,17 +397,18 @@
           }
 
-	  pmSourceExtFitPars *extFitPars = pmSourceExtFitParsAlloc();
-	  extFitPars->Mxx = source->moments->Mxx;
-	  extFitPars->Mxy = source->moments->Mxy;
-	  extFitPars->Myy = source->moments->Myy;
-	  extFitPars->Mrf = source->moments->Mrf;
-	  extFitPars->Mrh = source->moments->Mrh;
-	  extFitPars->peakMag = (source->peak->rawFlux > 0) ? -2.5*log10(source->peak->rawFlux) : NAN;
+	  // XXX deprecate?
+	  // XXX pmSourceExtFitPars *extFitPars = pmSourceExtFitParsAlloc();
+	  // XXX extFitPars->Mxx = source->moments->Mxx;
+	  // XXX extFitPars->Mxy = source->moments->Mxy;
+	  // XXX extFitPars->Myy = source->moments->Myy;
+	  // XXX extFitPars->Mrf = source->moments->Mrf;
+	  // XXX extFitPars->Mrh = source->moments->Mrh;
+	  // XXX extFitPars->peakMag = (source->peak->rawFlux > 0) ? -2.5*log10(source->peak->rawFlux) : NAN;
 	  
 	  // save kron mag, but assign apMag & psfMag on output (not yet calculated)
-	  extFitPars->krMag = -2.5*log10(source->moments->KronFlux);
-
-          psArrayAdd (source->extFitPars, 4, extFitPars);
-	  psFree (extFitPars);
+	  // XXX extFitPars->krMag = -2.5*log10(source->moments->KronFlux);
+
+          // psArrayAdd (source->extFitPars, 4, extFitPars);
+	  // psFree (extFitPars);
 
           // test for fit quality / result
Index: branches/eam_branches/ipp-20110710/psphot/src/psphotFindFootprints.c
===================================================================
--- branches/eam_branches/ipp-20110710/psphot/src/psphotFindFootprints.c	(revision 32157)
+++ branches/eam_branches/ipp-20110710/psphot/src/psphotFindFootprints.c	(revision 32205)
@@ -58,4 +58,15 @@
     psphotCullPeaks(readout, significance, recipe, detections->footprints);
     detections->peaks = pmFootprintArrayToPeaks(detections->footprints);
+
+    psphotFootprintSaddles (readout, detections->footprints);
+
+    int nSaddle = 0;
+    for (int i = 0; i < detections->peaks->n; i++) {
+	pmPeak *peak = detections->peaks->data[i];
+	
+	if (peak->saddlePoints) nSaddle += peak->saddlePoints->n;
+    }
+    fprintf (stderr, "%d saddle points for %ld peaks\n", nSaddle, detections->peaks->n);
+
     psLogMsg ("psphot", PS_LOG_INFO, "%ld peaks, %ld total footprints: %f sec\n", detections->peaks->n, detections->footprints->n, psTimerMark ("psphot.footprints"));
 
Index: branches/eam_branches/ipp-20110710/psphot/src/psphotFootprintSaddles.c
===================================================================
--- branches/eam_branches/ipp-20110710/psphot/src/psphotFootprintSaddles.c	(revision 32205)
+++ branches/eam_branches/ipp-20110710/psphot/src/psphotFootprintSaddles.c	(revision 32205)
@@ -0,0 +1,129 @@
+# include "psphotInternal.h"
+# ifndef ROUND
+# define ROUND(X) ((int) ((X) + 0.5*SIGN(X)))
+# endif
+
+psVector *getSaddlePoint (pmReadout *readout, int X1, int Y1, int X2, int Y2, bool swapcoords);
+
+// for each footprint with multiple detections, find the complete set of saddles, that is, the
+// local minima which lie on the lines between the peaks (same as a 'Coll' but less homophonic
+// with 'cull')
+bool psphotFootprintSaddles(pmReadout *readout, psArray *footprints) {  // array of pmFootprints
+    
+# if (PM_PEAKS_CULL_WITH_SMOOTHED_IMAGE)
+    psLogMsg ("psphot", PS_LOG_INFO, "Finding saddle points from footprints using the smoothed image");
+# else
+    psLogMsg ("psphot", PS_LOG_INFO, "Finding saddle points from footprints using the raw (unsmoothed) image");
+# endif
+
+    int xMax = readout->image->numCols - 1;
+    int yMax = readout->image->numRows - 1;
+
+    // FILE *f = fopen ("saddlepoints.dat", "w");
+
+    for (int i = 0; i < footprints->n; i++) {
+	pmFootprint *fp = footprints->data[i];
+	if (fp->peaks == NULL) continue;
+	if (fp->peaks->n < 2) continue;
+
+	// loop over all peaks in the footprint:
+	for (int is = 0; is < fp->peaks->n; is++) {
+
+	    pmPeak *startPeak = fp->peaks->data[is];
+	    if (!startPeak->saddlePoints) startPeak->saddlePoints = psArrayAllocEmpty (4);
+	
+	    // loop over all peaks in the footprint, skipping this one
+	    for (int ie = is + 1; ie < fp->peaks->n; ie++) {
+		if (ie == is) continue;
+
+		pmPeak *endPeak = fp->peaks->data[ie];
+		if (!endPeak->saddlePoints) endPeak->saddlePoints = psArrayAllocEmpty (4);
+
+		// step across all pixels between the startPeak and the endPeak, using a
+		// bresenham-like line-drawing operation
+		
+		int X1 = PS_MIN(xMax, PS_MAX(0, ROUND(startPeak->xf)));
+		int Y1 = PS_MIN(yMax, PS_MAX(0, ROUND(startPeak->yf)));
+		int X2 = PS_MIN(xMax, PS_MAX(0, ROUND(endPeak->xf)));
+		int Y2 = PS_MIN(yMax, PS_MAX(0, ROUND(endPeak->yf)));
+		
+		int dX = X2 - X1;
+		int dY = Y2 - Y1;
+
+		bool FlipCoords = (abs(dX) < abs(dY));
+		bool FlipDirect = FlipCoords ? (Y1 > Y2) : (X1 > X2);
+
+		psVector *saddlePoint = NULL;
+		if (!FlipDirect && !FlipCoords) saddlePoint = getSaddlePoint (readout, X1, Y1, X2, Y2, false);
+		if ( FlipDirect && !FlipCoords) saddlePoint = getSaddlePoint (readout, X2, Y2, X1, Y1, false);
+		if (!FlipDirect &&  FlipCoords) saddlePoint = getSaddlePoint (readout, Y1, X1, Y2, X2, true);
+		if ( FlipDirect &&  FlipCoords) saddlePoint = getSaddlePoint (readout, Y2, X2, Y1, X1, true);
+		psAssert (saddlePoint, "failed to get saddlePoint");
+
+		// fprintf (f, "%f %f | %d %d | %f %f\n", startPeak->xf, startPeak->yf, saddlePoint->data.S32[0], saddlePoint->data.S32[1], endPeak->xf, endPeak->yf);
+
+		psArrayAdd (startPeak->saddlePoints, 4, saddlePoint);
+		psArrayAdd (endPeak->saddlePoints, 4, saddlePoint);
+		psFree (saddlePoint);
+		// NOTE : all saddlePoints for the startPeak are also saddlePoints for the endPeak
+	    }
+	}
+    }
+    // fclose (f);
+    return true;
+}
+
+// psLogMsg ("psphot", PS_LOG_INFO, "%ld peaks, %ld total footprints: %f sec\n", detections->peaks->n, detections->footprints->n, 
+
+psVector *getSaddlePoint (pmReadout *readout, int X1, int Y1, int X2, int Y2, bool swapcoords) {
+
+    psAssert (X1 >= 0, "invalid");
+    psAssert (Y1 >= 0, "invalid");
+    psAssert (X1 < readout->image->numCols, "invalid");
+    psAssert (Y1 < readout->image->numRows, "invalid");
+
+    psAssert (X2 >= 0, "invalid");
+    psAssert (Y2 >= 0, "invalid");
+    psAssert (X2 < readout->image->numCols, "invalid");
+    psAssert (Y2 < readout->image->numRows, "invalid");
+
+    psVector *saddlePoint = psVectorAlloc(2, PS_TYPE_F32);
+
+    int dX = X2 - X1;
+    int dY = Y2 - Y1;
+		
+    float minValue = swapcoords ? readout->image->data.F32[X1][Y1] : readout->image->data.F32[Y1][X1];
+    saddlePoint->data.S32[0] = swapcoords ? Y1 : X1;
+    saddlePoint->data.S32[1] = swapcoords ? X1 : Y1;
+
+    int Y = Y1;
+    int e = 0;
+    for (int X = X1; X <= X2; X++) {
+	if (swapcoords) {
+	    float newValue = readout->image->data.F32[X][Y];
+	    if (newValue < minValue) {
+		minValue = newValue;
+		saddlePoint->data.S32[0] = Y;
+		saddlePoint->data.S32[1] = X;
+	    }
+	} else {
+	    float newValue = readout->image->data.F32[Y][X];
+	    if (newValue < minValue) {
+		minValue = newValue;
+		saddlePoint->data.S32[0] = X;
+		saddlePoint->data.S32[1] = Y;
+	    }
+	}
+	e += dY;
+	int e2 = 2 * e;
+	if (e2 > dX) {
+	    Y++;
+	    e -= dX;
+	} 
+	if (e2 < -dX) {
+	    Y--;
+	    e += dX;
+	}
+    }
+    return (saddlePoint);
+}
Index: branches/eam_branches/ipp-20110710/psphot/src/psphotKronIterate.c
===================================================================
--- branches/eam_branches/ipp-20110710/psphot/src/psphotKronIterate.c	(revision 32205)
+++ branches/eam_branches/ipp-20110710/psphot/src/psphotKronIterate.c	(revision 32205)
@@ -0,0 +1,421 @@
+# include "psphotInternal.h"
+# ifndef ROUND
+# define ROUND(X) ((int) ((X) + 0.5*SIGN(X)))
+# endif
+
+bool psphotKronWindowSetSource(pmSource *source, psImage *kronWindow, bool useKronRadius, bool insert);
+bool psphotKronWindowMag (pmSource *source, psImage *kronWindow, float radius, float minKronRadius, psImageMaskType maskVal);
+
+bool psphotKronIterate (pmConfig *config, const pmFPAview *view, const char *filerule)
+{
+    bool status = true;
+
+    // select the appropriate recipe information
+    psMetadata *recipe  = psMetadataLookupPtr (&status, config->recipes, PSPHOT_RECIPE);
+    psAssert (recipe, "missing recipe?");
+
+    int num = psphotFileruleCount(config, filerule);
+
+    // skip the chisq image (optionally?)
+    int chisqNum = psMetadataLookupS32 (&status, config->arguments, "PSPHOT.CHISQ.NUM");
+    if (!status) chisqNum = -1;
+
+    // loop over the available readouts
+    for (int i = 0; i < num; i++) {
+        if (i == chisqNum) continue; // skip chisq image
+
+        // find the currently selected readout
+        pmFPAfile *file = pmFPAfileSelectSingle(config->files, filerule, i); // File of interest
+        psAssert (file, "missing file?");
+
+        pmReadout *readout = pmFPAviewThisReadout(view, file->fpa);
+        psAssert (readout, "missing readout?");
+
+        pmDetections *detections = psMetadataLookupPtr (&status, readout->analysis, "PSPHOT.DETECTIONS");
+        psAssert (detections, "missing detections?");
+
+        psArray *sources = detections->allSources;
+        psAssert (sources, "missing sources?");
+
+        pmPSF *psf = psMetadataLookupPtr (&status, readout->analysis, "PSPHOT.PSF");
+        // psAssert (psf, "missing psf?");
+
+        if (!psphotKronIterateReadout (config, recipe, view, readout, sources, psf)) {
+            psError (PSPHOT_ERR_CONFIG, false, "failed to measure magnitudes for %s entry %d", filerule, i);
+            return false;
+        }
+    }
+    return true;
+}
+
+int psphotKapaChannel (int channel);
+bool psphotVisualShowMask (int kapaFD, psImage *inImage, const char *name, int channel);
+bool psphotVisualRangeImage (int kapaFD, psImage *inImage, const char *name, int channel, float min, float max);
+
+bool psphotKronIterateReadout(pmConfig *config, psMetadata *recipe, const pmFPAview *view, pmReadout *readout, psArray *sources, pmPSF *psf) {
+
+    bool status = false;
+
+    if (!sources->n) {
+        psLogMsg ("psphot", PS_LOG_INFO, "no sources, skipping masked kron");
+        return true;
+    }
+
+    psTimerStart ("psphot.kron");
+
+    // determine the number of allowed threads
+    int nThreads = psMetadataLookupS32(&status, config->arguments, "NTHREADS"); // Number of threads
+    if (!status) {
+        nThreads = 0;
+    }
+
+    float RADIUS = psMetadataLookupF32 (&status, readout->analysis, "PSF_MOMENTS_RADIUS");
+    if (!status) {
+        RADIUS = psMetadataLookupF32 (&status, recipe, "PSF_MOMENTS_RADIUS");
+    }
+
+    float MIN_KRON_RADIUS = psMetadataLookupF32 (&status, readout->analysis, "MOMENTS_MIN_KRON");
+    if (!status) {
+        MIN_KRON_RADIUS = 0.25*RADIUS;
+    }
+
+    float EXT_FIT_MAX_RADIUS = psMetadataLookupF32 (&status, recipe, "EXT_FIT_MAX_RADIUS");
+
+    // bit-masks to test for good/bad pixels
+    psImageMaskType maskVal = psMetadataLookupImageMask(&status, recipe, "MASK.PSPHOT");
+    assert (maskVal);
+
+    // bit-mask to mark pixels not used in analysis
+    psImageMaskType markVal = psMetadataLookupImageMask(&status, recipe, "MARK.PSPHOT");
+    assert (markVal);
+
+    // maskVal is used to test for rejected pixels, and must include markVal
+    maskVal |= markVal;
+
+    // psphotSaveImage (NULL, readout->mask, "kron.unmasked.fits");
+
+    pmSourcePhotometryMode photMode = PM_SOURCE_PHOT_PSFONLY;
+
+    // XXX tmp visualization
+    // int kapa = psphotKapaChannel (1);
+
+    // source analysis is done in S/N order (brightest first)
+    sources = psArraySort (sources, pmSourceSortByFlux);
+    if (!sources->n) {
+        psLogMsg ("psphot", PS_LOG_INFO, "no sources, skipping blend");
+        return true;
+    }
+
+    // generate the window image: multiply the flux by this to downweight neighbors
+    psImage *kronWindow = psImageAlloc (readout->image->numCols, readout->image->numRows, PS_TYPE_F32);
+    psImageInit (kronWindow, 1.0);
+
+    // start with the currently known moments (Mxx, Myy, Mxy) and generate a window image
+    for (int i = 0; i < sources->n; i++) {
+
+        pmSource *source = sources->data[i];
+        if (!source->peak) continue; // XXX how can we have a peak-less source?
+
+        // allocate space for moments
+        if (!source->moments) continue;
+
+	// set a window function for each source based on the moments
+	psphotKronWindowSetSource (source, kronWindow, false, true);
+    }
+    // psphotVisualRangeImage (kapa, kronWindow, "kronwin", 0, 0.0, 1.0);
+
+    psphotSaveImage (NULL, kronWindow, "kron.window.v0.fits");
+
+    for (int i = 0; i < sources->n; i++) {
+
+        pmSource *source = sources->data[i];
+        if (!source->peak) continue; // XXX how can we have a peak-less source?
+
+        // allocate space for moments
+        if (!source->moments) continue;
+
+	// replace object in image
+	if (source->tmpFlags & PM_SOURCE_TMPF_SUBTRACTED) {
+	    pmSourceAdd (source, PM_MODEL_OP_FULL, maskVal);
+	}
+
+	// iterate to the window radius
+	float windowRadius = PS_MIN(PS_MAX(RADIUS, 4.0*source->moments->Mrf), EXT_FIT_MAX_RADIUS);
+
+	// re-allocate image, weight, mask arrays for each peak with box big enough to fit BIG_RADIUS
+	pmSourceRedefinePixels (source, readout, source->peak->x, source->peak->y, windowRadius + 2);
+
+	// clear the window function for this source based on the moments
+	psphotKronWindowSetSource (source, kronWindow, (i > 0), false);
+	// psphotKronWindowSetSource (source, kronWindow, false, false);
+	// psphotVisualRangeImage (kapa, kronWindow, "kronwin", 1, 0.0, 1.0);
+
+	// this function populates moments->Mrf,KronFlux,KronFluxErr
+	psphotKronWindowMag (source, kronWindow, windowRadius, MIN_KRON_RADIUS, maskVal);
+	psImageMaskPixels (source->maskObj, "AND", PS_NOT_IMAGE_MASK(markVal));
+
+	// set a window function for each source based on the moments
+	psphotKronWindowSetSource (source, kronWindow, true, true);
+	// psphotKronWindowSetSource (source, kronWindow, false, true);
+
+	// test source fluxes
+        pmSourceMagnitudes (source, psf, photMode, maskVal, markVal, source->apRadius);
+	float kmag = -2.5*log10(source->moments->KronFlux);
+	if (source->psfMag - kmag > 0.25) {
+	    // fprintf (stderr, "continue\n");
+	}
+
+	// re-subtract the object, leave local sky
+	pmSourceSub (source, PM_MODEL_OP_FULL, maskVal);
+    }
+
+    psphotSaveImage (NULL, kronWindow, "kron.window.v1.fits");
+
+    for (int i = 0; i < sources->n; i++) {
+
+        pmSource *source = sources->data[i];
+        if (!source->peak) continue; // XXX how can we have a peak-less source?
+
+        // allocate space for moments
+        if (!source->moments) continue;
+
+	// replace object in image
+	if (source->tmpFlags & PM_SOURCE_TMPF_SUBTRACTED) {
+	    pmSourceAdd (source, PM_MODEL_OP_FULL, maskVal);
+	}
+
+	// iterate to the window radius
+	float windowRadius = PS_MIN(PS_MAX(RADIUS, 4.0*source->moments->Mrf), EXT_FIT_MAX_RADIUS);
+
+	// re-allocate image, weight, mask arrays for each peak with box big enough to fit BIG_RADIUS
+	pmSourceRedefinePixels (source, readout, source->peak->x, source->peak->y, windowRadius + 2);
+
+	// clear the window function for this source based on the moments
+	psphotKronWindowSetSource (source, kronWindow, (i > 0), false);
+	// psphotKronWindowSetSource (source, kronWindow, false, false);
+	// psphotVisualRangeImage (kapa, kronWindow, "kronwin", 1, 0.0, 1.0);
+
+	// this function populates moments->Mrf,KronFlux,KronFluxErr
+	psphotKronWindowMag (source, kronWindow, windowRadius, MIN_KRON_RADIUS, maskVal);
+	psImageMaskPixels (source->maskObj, "AND", PS_NOT_IMAGE_MASK(markVal));
+
+	// set a window function for each source based on the moments
+	psphotKronWindowSetSource (source, kronWindow, true, true);
+	// psphotKronWindowSetSource (source, kronWindow, false, true);
+
+	// test source fluxes
+        pmSourceMagnitudes (source, psf, photMode, maskVal, markVal, source->apRadius);
+	float kmag = -2.5*log10(source->moments->KronFlux);
+	if (source->psfMag - kmag > 0.25) {
+	    // fprintf (stderr, "continue\n");
+	}
+
+	// re-subtract the object, leave local sky
+	pmSourceSub (source, PM_MODEL_OP_FULL, maskVal);
+    }
+
+    psphotSaveImage (NULL, kronWindow, "kron.window.v2.fits");
+    psFree (kronWindow);
+
+    psLogMsg ("psphot.kron", PS_LOG_DETAIL, "measure masked kron magnitudes : %f sec for %ld objects\n", psTimerMark ("psphot.kron"), sources->n);
+
+    return true;
+}
+
+bool psphotKronWindowMag (pmSource *source, psImage *kronWindow, float radius, float minKronRadius, psImageMaskType maskVal) {
+
+    PS_ASSERT_PTR_NON_NULL(source, false);
+    PS_ASSERT_PTR_NON_NULL(source->peak, false);
+    PS_ASSERT_PTR_NON_NULL(source->pixels, false);
+    PS_ASSERT_FLOAT_LARGER_THAN(radius, 0.0, false);
+
+    psF32 R2 = PS_SQR(radius);
+
+    // a note about coordinates: coordinates of objects throughout psphot refer to the primary
+    // image coordinates.  the source->pixels image has an offset relative to its parent of
+    // col0,row0: a pixel (x,y) in the primary image has coordinates of (x-col0, y-row0) in
+    // this subimage.  we subtract off the peak coordinates, adjusted to this subimage, to have
+    // minimal round-off error in the sums.  since these values are subtracted just to minimize
+    // the dynamic range and are added back below, the exact value does not matter. these are
+    // (int) so they can be used in the image index below.
+
+    // Now calculate higher-order moments, using the above-calculated first moments to adjust coordinates
+    // Xn  = SUM (x - xc)^n * (z - sky)
+
+    psF32 RF = 0.0;
+    psF32 RS = 0.0;
+
+    // the peak position is less accurate but less subject to extreme deviations
+    float dX = source->moments->Mx - source->peak->xf;
+    float dY = source->moments->My - source->peak->yf;
+    float dR = hypot(dX, dY);
+    float Xo = (dR < 2.0) ? source->moments->Mx : source->peak->xf;
+    float Yo = (dR < 2.0) ? source->moments->My : source->peak->yf;
+
+    // center of mass in subimage.  Note: the calculation below uses pixel index, so we correct
+    // xCM, yCM from pixel coords to pixel index here.
+    psF32 xCM = Xo - 0.5 - source->pixels->col0; // coord of peak in subimage
+    psF32 yCM = Yo - 0.5 - source->pixels->row0; // coord of peak in subimage
+
+    int Xwo = source->pixels->col0;
+    int Ywo = source->pixels->row0;
+
+    for (psS32 row = 0; row < source->pixels->numRows ; row++) {
+
+	psF32 yDiff = row - yCM;
+	if (fabs(yDiff) > radius) continue;
+
+	psF32 *vPix = source->pixels->data.F32[row];
+	psF32 *vWin = &kronWindow->data.F32[row + Ywo][Xwo];
+
+	psImageMaskType *vMsk = (source->maskObj == NULL) ? NULL : source->maskObj->data.PS_TYPE_IMAGE_MASK_DATA[row];
+
+	for (psS32 col = 0; col < source->pixels->numCols ; col++, vPix++, vWin++) {
+	    if (vMsk) {
+		if (*vMsk & maskVal) {
+		    vMsk++;
+		    continue;
+		}
+		vMsk++;
+	    }
+	    if (isnan(*vPix)) continue;
+
+	    psF32 xDiff = col - xCM;
+	    if (fabs(xDiff) > radius) continue;
+
+	    // radius is just a function of (xDiff, yDiff)
+	    psF32 r2  = PS_SQR(xDiff) + PS_SQR(yDiff);
+	    if (r2 > R2) continue;
+
+	    // flux * window
+	    float weight  = *vWin;
+	    // float weight  = 1.0;
+	    psF32 pDiff = *vPix * weight;
+
+	    // Kron Flux uses the 1st radial moment (maybe Gaussian windowed?)
+	    psF32 rf = pDiff * sqrt(r2);
+	    psF32 rs = pDiff;
+
+	    RF  += rf;
+	    RS  += rs;
+	}
+    }
+
+    // Saturate the 1st radial moment
+    float Mrf = MAX(minKronRadius, RF/RS);
+    if (sqrt(source->peak->detValue) < 10.0) {
+	Mrf = MIN (radius, Mrf);
+    }
+
+    // Calculate the Kron magnitude (make this block optional?)
+    float radKron  = 2.5*Mrf;
+    float radKron2 = radKron*radKron;
+
+    int nKronPix = 0;
+    float Sum = 0.0;
+    float Var = 0.0;
+    float Win = 0.0;
+
+    for (psS32 row = 0; row < source->pixels->numRows ; row++) {
+
+	psF32 yDiff = row - yCM;
+	if (fabs(yDiff) > radKron) continue;
+
+	psF32 *vPix = source->pixels->data.F32[row];
+	psF32 *vWgt = source->variance->data.F32[row];
+	psF32 *vWin = &kronWindow->data.F32[row + Ywo][Xwo];
+
+	psImageMaskType *vMsk = (source->maskObj == NULL) ? NULL : source->maskObj->data.PS_TYPE_IMAGE_MASK_DATA[row];
+
+	for (psS32 col = 0; col < source->pixels->numCols ; col++, vPix++, vWgt++, vWin++) {
+	    if (vMsk) {
+		if (*vMsk & maskVal) {
+		    vMsk++;
+		    continue;
+		}
+		vMsk++;
+	    }
+	    if (isnan(*vPix)) continue;
+
+	    psF32 xDiff = col - xCM;
+	    if (fabs(xDiff) > radKron) continue;
+
+	    // radKron is just a function of (xDiff, yDiff)
+	    psF32 r2  = PS_SQR(xDiff) + PS_SQR(yDiff);
+	    if (r2 > radKron2) continue;
+
+	    float weight  = *vWin;
+	    // float weight  = 1.0;
+	    psF32 pDiff = *vPix * weight;
+	    psF32 wDiff = *vWgt * weight;
+
+	    Sum += pDiff;
+	    Var += wDiff;
+	    Win += weight;
+	    nKronPix ++;
+	}
+    }
+
+    source->moments->Mrf = Mrf;
+    source->moments->KronFlux    = Sum;
+    source->moments->KronFluxErr = Sum * (nKronPix / Win);
+    // source->moments->KronFluxErr = sqrt(Var / Win);
+
+    return true;
+}
+
+bool psphotKronWindowSetSource(pmSource *source, psImage *kronWindow, bool useKronRadius, bool insert) {
+
+    if (!source) return false;
+    if (!source->moments) return false;
+    psAssert(kronWindow, "need a window");
+
+    // we have a source with moments Mx, My, Mxx, Mxy, Myy.  we just need to define a Gaussian that has 
+    // these values (and peak of 1.0)
+
+    int Nx = kronWindow->numCols;
+    int Ny = kronWindow->numRows;
+
+    float Xo = source->moments->Mx;
+    float Yo = source->moments->My;
+
+    float Mxx = source->moments->Mxx;
+    float Mxy = source->moments->Mxy;
+    float Myy = source->moments->Myy;
+
+    float Mmajor = 0.5*(Mxx + Myy) + 0.5*sqrt(PS_SQR(Mxx - Myy) + 4.0*PS_SQR(Mxy));
+    float Mminor = 0.5*(Mxx + Myy) - 0.5*sqrt(PS_SQR(Mxx - Myy) + 4.0*PS_SQR(Mxy));
+
+    // float scale = useKronRadius ? 2.0 * source->moments->Mrf / Mmajor : 2.0;
+    float scale = 3.0 * source->moments->Mrf / Mmajor;
+    // float scale = 2.0;
+
+    float Sxx = scale * Mmajor * Mminor / Myy; // sigma_x^2
+    float Sxy = Mxy / (scale * Mmajor * Mminor);
+    float Syy = scale * Mmajor * Mminor / Mxx; // sigma_y^2
+
+    float Smajor = sqrt(Mmajor);
+
+    int minX = PS_MIN(PS_MAX(Xo - 5*Smajor, 0), Nx - 1);
+    int maxX = PS_MIN(PS_MAX(Xo + 5*Smajor, 0), Nx - 1);
+    int minY = PS_MIN(PS_MAX(Yo - 5*Smajor, 0), Ny - 1);
+    int maxY = PS_MIN(PS_MAX(Yo + 5*Smajor, 0), Ny - 1);
+
+    float rMxx = 0.5 / Sxx;
+    float rMyy = 0.5 / Syy;
+    
+    for (int iy = minY; iy < maxY; iy++) {
+	for (int ix = minX; ix < maxX; ix++) {
+
+	    float dX = (ix + 0.5 - Xo);
+	    float dY = (iy + 0.5 - Yo);
+
+	    float z = rMxx * PS_SQR(dX) + rMyy * PS_SQR(dY) + Sxy*dX*dY;
+
+	    float f = insert ? 1.001 - exp(-z) : 1.0 / (1.001 - exp(-z));
+
+	    kronWindow->data.F32[iy][ix] *= f;
+	}
+    }
+    return true;
+}
Index: branches/eam_branches/ipp-20110710/psphot/src/psphotKronMasked.c
===================================================================
--- branches/eam_branches/ipp-20110710/psphot/src/psphotKronMasked.c	(revision 32157)
+++ branches/eam_branches/ipp-20110710/psphot/src/psphotKronMasked.c	(revision 32205)
@@ -1,4 +1,8 @@
 # include "psphotInternal.h"
-
+# ifndef ROUND
+# define ROUND(X) ((int) ((X) + 0.5*SIGN(X)))
+# endif
+
+bool psphotKronMask (pmSource *source, psImage *kronMask, psImageMaskType markVal);
 bool psphotKronMag (pmSource *source, float radius, float minKronRadius, psImageMaskType maskVal);
 
@@ -45,4 +49,8 @@
 }
 
+int psphotKapaChannel (int channel);
+bool psphotVisualShowMask (int kapaFD, psImage *inImage, const char *name, int channel);
+bool psphotVisualRangeImage (int kapaFD, psImage *inImage, const char *name, int channel, float min, float max);
+
 bool psphotKronMaskedReadout(pmConfig *config, psMetadata *recipe, const pmFPAview *view, pmReadout *readout, psArray *sources, pmPSF *psf) {
 
@@ -72,4 +80,6 @@
     }
 
+    float EXT_FIT_MAX_RADIUS = psMetadataLookupF32 (&status, recipe, "EXT_FIT_MAX_RADIUS");
+
     // bit-masks to test for good/bad pixels
     psImageMaskType maskVal = psMetadataLookupImageMask(&status, recipe, "MASK.PSPHOT");
@@ -83,9 +93,14 @@
     maskVal |= markVal;
 
+    // psphotSaveImage (NULL, readout->mask, "kron.unmasked.fits");
+
+    pmSourcePhotometryMode photMode = PM_SOURCE_PHOT_PSFONLY;
+
+    // XXX tmp visualization
+    // int kapa = psphotKapaChannel (1);
+
     // generate the mask image: increment counter for every source overlapping the pixel
-    psImage *kronMask = psImageAlloc (readout->image->numCols, readout->image->numRows, PS_TYPE_U8);
+    psImage *kronMask = psImageAlloc (readout->image->numCols, readout->image->numRows, PS_TYPE_S32);
     psImageInit (kronMask, 0);
-    int Nx = kronMask->numCols;
-    int Ny = kronMask->numRows;
     for (int i = 0; i < sources->n; i++) {
 
@@ -96,56 +111,59 @@
         if (!source->moments) continue;
 
+	// replace object in image
+	if (source->tmpFlags & PM_SOURCE_TMPF_SUBTRACTED) {
+	    pmSourceAdd (source, PM_MODEL_OP_FULL, maskVal);
+	}
+
+	// save the PSF-based values
+	// source->moments->KronFluxPSF    = source->moments->KronFlux;
+	// source->moments->KronFluxPSFErr = source->moments->KronFluxErr;
+	// source->moments->KronRadiusPSF  = source->moments->Mrf;
+
+	// iterate to the window radius
 	float windowRadius = RADIUS;
-	if (source->moments->KronFlux / source->moments->KronFluxErr > 10.0) {
-	    windowRadius = PS_MAX (RADIUS, 10*source->moments->Mrf);
-	}
-
-	// re-allocate image, weight, mask arrays for each peak with box big enough to fit BIG_RADIUS
-	pmSourceRedefinePixels (source, readout, source->peak->x, source->peak->y, windowRadius + 2);
-
-	pmMoments psfMoments = *source->moments;
-
-        // replace object in image
-        if (source->tmpFlags & PM_SOURCE_TMPF_SUBTRACTED) {
-            pmSourceAdd (source, PM_MODEL_OP_FULL, maskVal);
-        }
-
-	psphotKronMag (source, windowRadius, MIN_KRON_RADIUS, maskVal);
-
-        // re-subtract the object, leave local sky
-        pmSourceSub (source, PM_MODEL_OP_FULL, maskVal);
-
-	continue;
-
-	// XXX skip this code
-	// generate the pixel masks
-	// int Xo = source->moments->Mx;
-	// int Yo = source->moments->My;
-	float dX = source->moments->Mx - source->peak->xf;
-	float dY = source->moments->My - source->peak->yf;
-	float dR = hypot(dX, dY);
-	
-	float Xo = (dR < 2.0) ? source->moments->Mx : source->peak->xf;
-	float Yo = (dR < 2.0) ? source->moments->My : source->peak->yf;
-	
-	int Kr = 2.5*source->moments->Mrf;
-	int Kr2 = Kr*Kr;
-	
-	for (int iy = Yo - Kr; iy < Yo + Kr + 1; iy++) {
-	    if (iy < 0) continue;
-	    if (iy >= Ny) continue;
-	    for (int ix = Xo - Kr; ix < Xo + Kr + 1; ix++) {
-		if (ix < 0) continue;
-		if (ix >= Nx) continue;
-
-		if (PS_SQR(ix - Xo) + PS_SQR(iy - Yo) > Kr2) continue;
-		if (kronMask->data.U8[iy][ix] < 0xff) {
-		    kronMask->data.U8[iy][ix] ++;
-		}
-	    }
-	}
-    }		
-
+	for (int j = 0; j < 4; j++) {
+	    // XXX use some S/N criterion to limit?
+	    // if (source->moments->KronFlux / source->moments->KronFluxErr > 10.0) {
+	    // windowRadius = PS_MAX (RADIUS, 10*source->moments->Mrf);
+	    // }
+
+	    // re-allocate image, weight, mask arrays for each peak with box big enough to fit BIG_RADIUS
+	    pmSourceRedefinePixels (source, readout, source->peak->x, source->peak->y, windowRadius + 2);
+
+	    // mask the pixels not contained by the footprint
+	    // if (psphotMaskFootprint (readout, source, markVal)) {
+	    // 	// psphotVisualShowMask (kapa, source->maskObj, "mask", 2);
+	    // 	// psphotVisualRangeImage (kapa, source->pixels, "image", 1, -50, 300);
+	    // 	// fprintf (stderr, "masked\n");
+	    // }	    
+
+	    // mask the pixels associated with near neighbors
+	    if (psphotKronMask (source, kronMask, markVal)) {
+		// psphotVisualShowMask (kapa, source->maskObj, "mask", 2);
+		// psphotVisualRangeImage (kapa, source->pixels, "image", 1, -50, 300);
+	     	// fprintf (stderr, "masked\n");
+	    }	    
+	    // this function populates moments->Mrf,KronFlux,KronFluxErr
+	    // XXX what about KronInner, KronOuter, etc?
+	    psphotKronMag (source, windowRadius, MIN_KRON_RADIUS, maskVal);
+	    windowRadius = PS_MIN(PS_MAX(RADIUS, 4.0*source->moments->Mrf), EXT_FIT_MAX_RADIUS);
+	    psImageMaskPixels (source->maskObj, "AND", PS_NOT_IMAGE_MASK(markVal));
+	}
+        pmSourceMagnitudes (source, psf, photMode, maskVal, markVal, source->apRadius);
+	float kmag = -2.5*log10(source->moments->KronFlux);
+	if (source->psfMag - kmag > 0.25) {
+	    // psphotVisualShowMask (kapa, source->maskObj, "mask", 1);
+	    // psphotVisualRangeImage (kapa, source->pixels, "image", 0, -50, 300);
+	    // fprintf (stderr, "psf: %f, kron: %f, dmag: %f\n", source->psfMag, kmag, source->psfMag - kmag);
+	    // fprintf (stderr, "continue\n");
+	}
+	// re-subtract the object, leave local sky
+	pmSourceSub (source, PM_MODEL_OP_FULL, maskVal);
+    }
+
+    // psphotSaveImage (NULL, readout->mask, "kron.masked.fits");
     // psphotSaveImage (NULL, kronMask, "kronmask.fits");
+
     psLogMsg ("psphot.kron", PS_LOG_DETAIL, "measure masked kron magnitudes : %f sec for %ld objects\n", psTimerMark ("psphot.kron"), sources->n);
 
@@ -153,4 +171,6 @@
     return true;
 }
+
+# define WEIGHTED 0
 
 bool psphotKronMag (pmSource *source, float radius, float minKronRadius, psImageMaskType maskVal) {
@@ -161,7 +181,9 @@
     PS_ASSERT_FLOAT_LARGER_THAN(radius, 0.0, false);
 
-    psF32 Sum = 0.0;
-    psF32 Var = 0.0;
     psF32 R2 = PS_SQR(radius);
+
+# if (WEIGHTED)
+    float rsigma2 = 0.5 / PS_SQR(radius / 4.0); // use a Gaussian window with sigma = R_window / 2
+# endif
 
     // a note about coordinates: coordinates of objects throughout psphot refer to the primary
@@ -179,5 +201,5 @@
     psF32 RS = 0.0;
 
-# if (1)
+    // the peak position is less accurate but less subject to extreme deviations
     float dX = source->moments->Mx - source->peak->xf;
     float dY = source->moments->My - source->peak->yf;
@@ -185,8 +207,4 @@
     float Xo = (dR < 2.0) ? source->moments->Mx : source->peak->xf;
     float Yo = (dR < 2.0) ? source->moments->My : source->peak->yf;
-# else
-    float Xo = source->moments->Mx;
-    float Yo = source->moments->My;
-# endif
 
     // center of mass in subimage.  Note: the calculation below uses pixel index, so we correct
@@ -221,9 +239,14 @@
 	    if (r2 > R2) continue;
 
+# if (WEIGHTED)
+	    float z = r2 * rsigma2;
+	    assert (z >= 0.0);
+	    float weight  = exp(-z);
+	    psF32 pDiff = *vPix * weight;
+# else
 	    psF32 pDiff = *vPix;
-
-	    Sum += pDiff;
-
-	    // Kron Flux uses the 1st radial moment (NOT Gaussian windowed)
+# endif
+
+	    // Kron Flux uses the 1st radial moment (maybe Gaussian windowed?)
 	    psF32 rf = pDiff * sqrt(r2);
 	    psF32 rs = pDiff;
@@ -245,5 +268,6 @@
 
     int nKronPix = 0;
-    Sum = Var = 0.0;
+    float Sum = 0.0;
+    float Var = 0.0;
 
     for (psS32 row = 0; row < source->pixels->numRows ; row++) {
@@ -274,6 +298,14 @@
 	    if (r2 > radKron2) continue;
 
+# if (WEIGHTED)
+	    float z = r2 * rsigma2;
+	    assert (z >= 0.0);
+	    float weight  = exp(-z);
+	    psF32 pDiff = *vPix * weight;
+	    psF32 wDiff = *vWgt * weight;
+# else
 	    psF32 pDiff = *vPix;
 	    psF32 wDiff = *vWgt;
+# endif
 
 	    Sum += pDiff;
@@ -283,13 +315,173 @@
     }
 
-    source->moments->Mrh = Mrf;
-
-    // XXX for a test, save the old values here:
-    source->moments->KronFouter    = source->moments->KronFlux;
-    source->moments->KronCoreErr = source->moments->KronFluxErr;
-
-    source->moments->KronFlux    = Sum       * M_PI * PS_SQR(radKron) / nKronPix;
-    source->moments->KronFluxErr = sqrt(Var) * M_PI * PS_SQR(radKron) / nKronPix;
+    source->moments->Mrf = Mrf;
+    source->moments->KronFlux    = Sum;
+    source->moments->KronFluxErr = sqrt(Var);
+    // source->moments->KronFlux    = Sum       * M_PI * PS_SQR(radKron) / nKronPix;
+    // source->moments->KronFluxErr = sqrt(Var) * M_PI * PS_SQR(radKron) / nKronPix;
 
     return true;
 }
+
+enum {Y_P, Y_M, X_P, X_M};
+
+bool psphotKronMask (pmSource *source, psImage *kronMask, psImageMaskType markVal) {
+
+    int X, Y, e, dXs, dYs;
+
+    if (!source->peak) return false;
+    if (!source->peak->saddlePoints) return false;
+    if (!source->peak->saddlePoints->n) return false;
+
+    // return true;
+
+    int xMax = source->maskObj->numCols;
+    int yMax = source->maskObj->numRows;
+    // int xOff = source->maskObj->col0;
+    // int yOff = source->maskObj->row0;
+
+    int Xo = ROUND(source->peak->xf) - source->maskObj->col0;
+    int Yo = ROUND(source->peak->yf) - source->maskObj->row0;
+
+    for (int i = 0; i < source->peak->saddlePoints->n; i++) {
+
+	psVector *saddlePoint = source->peak->saddlePoints->data[i];
+	int Xs = ROUND(saddlePoint->data.S32[0]) - source->maskObj->col0;
+	int Ys = ROUND(saddlePoint->data.S32[1]) - source->maskObj->row0;
+	
+	// fprintf (stderr, "mask %d,%d @ %d,%d (%d,%d)\n", Xo, Yo, Xs, Ys, saddlePoint->data.S32[0], saddlePoint->data.S32[1]);
+
+	// We want to mask the pixels between the edge of the image and the line perpendicular
+	// to the connecting line.  Call dX,dY = (Xs-Xo,Ys-Yo).  There are 4 cases: 
+	// +y : dY > 0, fabs(dY) > fabs(dX)
+	// -y : dY < 0, fabs(dY) > fabs(dX)
+	// +x : dX > 0, fabs(dX) > fabs(dY)
+	// -x : dX < 0, fabs(dX) > fabs(dY)
+
+	int dX = Xs - Xo;
+	int dY = Ys - Yo;
+
+	// +y : dY > 0, fabs(dY) > fabs(dX)
+	// -y : dY < 0, fabs(dY) > fabs(dX)
+	if (fabs(dY) > fabs(dX)) {
+	    // points to right of saddle
+	    Y = Ys;
+	    e = 0;
+	    dXs = (dY > 0) ? +dY : -dY;  // this forces dXs > 0
+	    dYs = (dY > 0) ? -dX : +dX;
+	    for (X = Xs; X < xMax; X++) {
+		if (dY > 0) {
+		    for (int Ym = Y; Ym < yMax; Ym++) {
+			source->maskObj->data.PS_TYPE_IMAGE_MASK_DATA[Ym][X] |= markVal;
+			// kronMask->data.S32[Ym+yOff][X+xOff] = source->id;
+		    }
+		} else {
+		    for (int Ym = Y; Ym >= 0; Ym--) {
+			source->maskObj->data.PS_TYPE_IMAGE_MASK_DATA[Ym][X] |= markVal;
+			// kronMask->data.S32[Ym+yOff][X+xOff] = source->id;
+		    }
+		}
+		e += dYs;
+		int e2 = 2 * e;
+		if (e2 > dXs) {
+		    Y ++;
+		    e -= dXs;
+		}
+		if (e2 < -dXs) {
+		    Y --;
+		    e += dXs;
+		}
+		if (Y >= yMax) break;
+		if (Y < 0) break;
+	    }
+	    // points to left of saddle
+	    Y = Ys;
+	    e = 0;
+	    for (X = Xs; X >= 0; X--) {
+		if (dY > 0) {
+		    for (int Ym = Y; Ym < yMax; Ym++) {
+			source->maskObj->data.PS_TYPE_IMAGE_MASK_DATA[Ym][X] |= markVal;
+			// kronMask->data.S32[Ym+yOff][X+xOff] = source->id;
+		    }
+		} else {
+		    for (int Ym = Y; Ym >= 0; Ym--) {
+			source->maskObj->data.PS_TYPE_IMAGE_MASK_DATA[Ym][X] |= markVal;
+			// kronMask->data.S32[Ym+yOff][X+xOff] = source->id;
+		    }
+		}
+		e -= dYs;
+		int e2 = 2 * e;
+		if (e2 > dXs) {
+		    Y ++;
+		    e -= dXs;
+		}
+		if (e2 < -dXs) {
+		    Y --;
+		    e += dXs;
+		}
+		if (Y >= yMax) break;
+		if (Y < 0) break;
+	    }
+	} else {
+	    // points to right of saddle
+	    X = Xs;
+	    e = 0;
+	    dXs = (dX > 0) ? -dY : +dY;
+	    dYs = (dX > 0) ? +dX : -dX;  // this forces dYs > 0
+	    for (Y = Ys; Y < yMax; Y++) {
+		if (dX > 0) {
+		    for (int Xm = X; Xm < xMax; Xm++) {
+			source->maskObj->data.PS_TYPE_IMAGE_MASK_DATA[Y][Xm] |= markVal;
+			// kronMask->data.S32[Y+yOff][Xm+xOff] = source->id;
+		    }
+		} else {
+		    for (int Xm = X; Xm >= 0; Xm--) {
+			source->maskObj->data.PS_TYPE_IMAGE_MASK_DATA[Y][Xm] |= markVal;
+			// kronMask->data.S32[Y+yOff][Xm+xOff] = source->id;
+		    }
+		}
+		e += dXs;
+		int e2 = 2 * e;
+		if (e2 > dYs) {
+		    X ++;
+		    e -= dYs;
+		}
+		if (e2 < -dYs) {
+		    X --;
+		    e += dYs;
+		}
+		if (X >= yMax) break;
+		if (X < 0) break;
+	    }
+	    // points to left of saddle
+	    X = Xs;
+	    e = 0;
+	    for (Y = Ys; Y >= 0; Y--) {
+		if (dX > 0) {
+		    for (int Xm = X; Xm < xMax; Xm++) {
+			source->maskObj->data.PS_TYPE_IMAGE_MASK_DATA[Y][Xm] |= markVal;
+			// kronMask->data.S32[Y+yOff][Xm+xOff] = source->id;
+		    }
+		} else {
+		    for (int Xm = X; Xm >= 0; Xm--) {
+			source->maskObj->data.PS_TYPE_IMAGE_MASK_DATA[Y][Xm] |= markVal;
+			// kronMask->data.S32[Y+yOff][Xm+xOff] = source->id;
+		    }
+		}
+		e -= dXs;
+		int e2 = 2 * e;
+		if (e2 > dYs) {
+		    X ++;
+		    e -= dYs;
+		}
+		if (e2 < -dYs) {
+		    X --;
+		    e += dYs;
+		}
+		if (X >= yMax) break;
+		if (X < 0) break;
+	    }
+	}
+    }
+    return true;
+}
Index: branches/eam_branches/ipp-20110710/psphot/src/psphotModelTestReadout.c
===================================================================
--- branches/eam_branches/ipp-20110710/psphot/src/psphotModelTestReadout.c	(revision 32157)
+++ branches/eam_branches/ipp-20110710/psphot/src/psphotModelTestReadout.c	(revision 32205)
@@ -84,5 +84,5 @@
 
     // define the source of interest
-    float xObj     = psMetadataLookupF32 (&status, recipe, "TEST_FIT_X");
+   float xObj     = psMetadataLookupF32 (&status, recipe, "TEST_FIT_X");
     float yObj     = psMetadataLookupF32 (&status, recipe, "TEST_FIT_Y");
     if (!isfinite(xObj) || !isfinite(yObj)) psAbort ("object position is not defined");
@@ -102,4 +102,21 @@
     status = pmSourceLocalSky (source, PS_STAT_SAMPLE_MEDIAN, INNER, maskVal, markVal);
     if (!status) psAbort("pmSourceLocalSky error");
+
+    { 
+	// XXX I want to test an iterative aperture for brighter sources
+	float radius = mRADIUS;
+	for (int i = 0; i < 10; i++) {
+
+	    // get the source moments
+	    status = pmSourceMoments (source, radius, 0.0, 0.0, MIN_KRON_RADIUS, maskVal);
+	    if (!status) psAbort("psSourceMoments error");
+
+	    float oldRadius = radius;
+	    radius = source->moments->Mrf * RADIUS;
+	    
+	    fprintf (stderr, "%d %f  %f  %f\n", i, oldRadius, radius, source->moments->Mrf);
+	}
+	exit (0);
+    }
 
     // get the source moments
Index: branches/eam_branches/ipp-20110710/psphot/src/psphotRadiusChecks.c
===================================================================
--- branches/eam_branches/ipp-20110710/psphot/src/psphotRadiusChecks.c	(revision 32157)
+++ branches/eam_branches/ipp-20110710/psphot/src/psphotRadiusChecks.c	(revision 32205)
@@ -241,4 +241,46 @@
 }
 
+// call this function whenever you (re)-define the EXT model
+bool psphotMaskFootprint (pmReadout *readout, pmSource *source, psImageMaskType markVal) {
+
+    psAssert (source, "source not defined??");
+    psAssert (source->peak, "peak not defined??");
+
+    pmPeak *peak = source->peak;
+
+    // set the radius based on the footprint:
+    if (!peak->footprint) return false;
+    pmFootprint *footprint = peak->footprint;
+    if (!footprint->spans) return false;
+    if (footprint->spans->n < 1) return false;
+
+    int Xo = source->maskObj->col0;
+    int Yo = source->maskObj->row0;
+
+    // mark all pixels 
+    for (int j = 0; j < source->maskObj->numRows; j++) {
+	for (int i = 0; i < source->maskObj->numCols; i++) {
+	    source->maskObj->data.PS_TYPE_IMAGE_MASK_DATA[j][i] |= markVal;
+	}
+    }
+
+    psImageMaskType clearVal = PS_NOT_IMAGE_MASK(markVal);
+
+    for (int j = 0; j < footprint->spans->n; j++) {
+        pmSpan *span = footprint->spans->data[j];
+
+	// mask the rows before and after each span
+	int minX = span->x0 - Xo - 2;
+	int maxX = span->x1 - Xo + 2;
+	int myY = span->y - Yo;
+
+	// unmark pixels inside the footprint
+	for (int i = minX; i <= maxX; i++) {
+	    source->maskObj->data.PS_TYPE_IMAGE_MASK_DATA[myY][i] &= clearVal;
+	}
+    }
+    return true;
+}
+
 // alternative EXT radius based on model guess (for use without footprints)
 bool psphotSetRadiusModel (pmModel *model, pmReadout *readout, pmSource *source, psImageMaskType markVal, bool deep) {
Index: branches/eam_branches/ipp-20110710/psphot/src/psphotReadout.c
===================================================================
--- branches/eam_branches/ipp-20110710/psphot/src/psphotReadout.c	(revision 32157)
+++ branches/eam_branches/ipp-20110710/psphot/src/psphotReadout.c	(revision 32205)
@@ -205,5 +205,6 @@
     
     // but this is chosen above to be appropriate for the PSF objects (not galaxies)
-    psphotKronMasked(config, view, filerule);
+    // psphotKronMasked(config, view, filerule);
+    psphotKronIterate(config, view, filerule);
 
     // identify CRs and extended sources (only unmeasured sources are measured)
@@ -321,5 +322,5 @@
     // XXX re-measure the kron mags with models subtracted
     // Note that this uses the PSF_MOMENTS_RADIUS as a window
-    psphotKronMasked(config, view, filerule);
+    // psphotKronMasked(config, view, filerule);
 
     // measure source size for the remaining sources
Index: branches/eam_branches/ipp-20110710/psphot/src/psphotSersicModelClass.c
===================================================================
--- branches/eam_branches/ipp-20110710/psphot/src/psphotSersicModelClass.c	(revision 32157)
+++ branches/eam_branches/ipp-20110710/psphot/src/psphotSersicModelClass.c	(revision 32205)
@@ -542,5 +542,5 @@
     
     // get the PSFratio and the ASPratio : we use these to choose the model class
-    float PSFratio = momentAxes.major / (2.35 * psfAxes.major) ;
+    float PSFratio = source->moments->Mrf / (2.35 * psfAxes.major) ;
     float ASPratio = momentAxes.minor / momentAxes.major;
 
Index: branches/eam_branches/ipp-20110710/psphot/src/psphotSourceSize.c
===================================================================
--- branches/eam_branches/ipp-20110710/psphot/src/psphotSourceSize.c	(revision 32157)
+++ branches/eam_branches/ipp-20110710/psphot/src/psphotSourceSize.c	(revision 32205)
@@ -209,5 +209,5 @@
         pmSourceMagnitudes (source, psf, photMode, maskVal, markVal, source->apRadius);
 
-        float kMag = -2.5*log10(source->moments->KronFlux);
+        float kMag = -2.5*log10(source->moments->KronFluxPSF);
         float dMag = source->psfMag - kMag;
 
@@ -334,5 +334,5 @@
         psF32 Mxy = source->moments->Mxy;
 
-        float KronMag = -2.5*log10(source->moments->KronFlux);
+        float KronMag = -2.5*log10(source->moments->KronFluxPSF);
 
 	float Mminor = 0.5*(Mxx + Myy) - 0.5*sqrt(PS_SQR(Mxx - Myy) + 4.0*PS_SQR(Mxy));
@@ -509,5 +509,5 @@
         pmSourceSub (source, PM_MODEL_OP_FULL, options->maskVal);
 
-        float kMag = -2.5*log10(source->moments->KronFlux);
+        float kMag = -2.5*log10(source->moments->KronFluxPSF);
         float dMag = source->psfMag - kMag;
 
Index: branches/eam_branches/ipp-20110710/psphot/src/psphotSourceStats.c
===================================================================
--- branches/eam_branches/ipp-20110710/psphot/src/psphotSourceStats.c	(revision 32157)
+++ branches/eam_branches/ipp-20110710/psphot/src/psphotSourceStats.c	(revision 32205)
@@ -554,5 +554,5 @@
         // determine the PSF parameters from the source moment values
         pmPSFClump psfClump = pmSourcePSFClump (NULL, NULL, sources, PSF_SN_LIM, PSF_CLUMP_GRID_SCALE, MOMENTS_SX_MAX, MOMENTS_SY_MAX, MOMENTS_SX_MIN, MOMENTS_SY_MIN, MOMENTS_AR_MAX);
-        psLogMsg ("psphot", 3, "radius %.1f, nStars: %d of %d in clump, nSigma: %5.2f, X,  Y: %f, %f (%f, %f)\n", sigma[i], psfClump.nStars, psfClump.nTotal, psfClump.nSigma, psfClump.X, psfClump.Y, sqrt(psfClump.X) / sigma[i], sqrt(psfClump.Y) / sigma[i]);
+        psLogMsg ("psphot", 3, "sigma guess (pix) %.1f, nStars: %d of %d in clump, nSigma: %5.2f, X,  Y: %f, %f (%f, %f)\n", sigma[i], psfClump.nStars, psfClump.nTotal, psfClump.nSigma, psfClump.X, psfClump.Y, sqrt(psfClump.X) / sigma[i], sqrt(psfClump.Y) / sigma[i]);
 
 	Rmin[i] = pmSourceMinKronRadius(sources, PSF_SN_LIM);