Index: /trunk/psphot/src/Makefile.am
===================================================================
--- /trunk/psphot/src/Makefile.am	(revision 31451)
+++ /trunk/psphot/src/Makefile.am	(revision 31452)
@@ -167,4 +167,5 @@
 	psphotSourcePlots.c	       \
 	psphotRadialPlot.c	       \
+	psphotKronMasked.c             \
 	psphotDeblendSatstars.c	       \
 	psphotMosaicSubimage.c	       \
Index: /trunk/psphot/src/psphot.h
===================================================================
--- /trunk/psphot/src/psphot.h	(revision 31451)
+++ /trunk/psphot/src/psphot.h	(revision 31452)
@@ -149,4 +149,7 @@
 bool            psphotPSFstats (pmReadout *readout, pmPSF *psf);
 bool            psphotMomentsStats (pmReadout *readout, psArray *sources);
+
+bool            psphotKronMasked (pmConfig *config, const pmFPAview *view, const char *filerule);
+bool            psphotKronMaskedReadout(pmConfig *config, psMetadata *recipe, const pmFPAview *view, pmReadout *readout, psArray *sources, pmPSF *psf);
 
 // in psphotGuessModel.c
@@ -310,5 +313,5 @@
 
 bool psphotMakeFluxScale (psImage *image, psMetadata *recipe, pmPSF *psf);
-bool psphotMakeGrowthCurve (pmReadout *readout, psMetadata *recipe, pmPSF *psf);
+bool psphotMakeGrowthCurve (pmReadout *readout, psMetadata *recipe, pmPSF *psf, psArray *sources);
 bool psphotDumpPSFStars (pmReadout *readout, pmPSFtry *try, float radius, psImageMaskType maskVal, psImageMaskType markVal);
 
Index: /trunk/psphot/src/psphotApResid.c
===================================================================
--- /trunk/psphot/src/psphotApResid.c	(revision 31451)
+++ /trunk/psphot/src/psphotApResid.c	(revision 31452)
@@ -1,4 +1,4 @@
 # include "psphotInternal.h"
-// # define DEBUG
+# define DEBUG
 
 # define SKIPSTAR(MSG) { psTrace ("psphot", 3, "invalid : %s", MSG); continue; }
@@ -206,5 +206,5 @@
 
         // XXX make this user-configurable?
-        if (source->errMag > 0.01) continue;
+        if (source->psfMagErr > 0.01) continue;
 
         // aperture residual for this source
@@ -222,5 +222,5 @@
                  source->peak->xf, source->peak->yf,
                  source->modelPSF->params->data.F32[PM_PAR_XPOS], source->modelPSF->params->data.F32[PM_PAR_YPOS],
-                 source->psfMag, source->apMag, source->errMag,
+                 source->psfMag, source->apMag, source->psfMagErr,
                  source->modelPSF->params->data.F32[PM_PAR_I0],
                  source->modelPSF->params->data.F32[PM_PAR_SXX], source->modelPSF->params->data.F32[PM_PAR_SXY], source->modelPSF->params->data.F32[PM_PAR_SYY],
@@ -228,5 +228,5 @@
 # endif
         if (!isfinite(source->psfMag)) psAbort ("nan in psfMag");
-        if (!isfinite(source->errMag)) psAbort ("nan in errMag");
+        if (!isfinite(source->psfMagErr)) psAbort ("nan in psfMagErr");
         if (!isfinite(source->apMag)) psAbort ("nan in apMag");
         if (!isfinite(model->params->data.F32[PM_PAR_XPOS])) psAbort ("nan in xPos");
@@ -234,5 +234,5 @@
 
         psVectorAppend (mag, source->psfMag);
-        psVectorAppend (dMag,source->errMag);
+        psVectorAppend (dMag,source->psfMagErr);
         psVectorAppend (apResid, dap);
         psVectorAppend (xPos, model->params->data.F32[PM_PAR_XPOS]);
@@ -328,5 +328,21 @@
     psMetadataAddF32 (readout->analysis, PS_LIST_TAIL, "DAPMIFIT", PS_META_REPLACE, "ap residual scatter", psf->dApResid);
     psMetadataAddS32 (readout->analysis, PS_LIST_TAIL, "NAPMIFIT", PS_META_REPLACE, "number of apresid stars", psf->nApResid);
-    psMetadataAddF32 (readout->analysis, PS_LIST_TAIL, "APLOSS",   PS_META_REPLACE, "aperture loss (mag)", psf->growth ? psf->growth->apLoss : NAN);
+
+    // curve-of-growth offset
+    if (psf->growth) {
+	float gaussSigma = psMetadataLookupF32(&status, readout->analysis, "MOMENTS_GAUSS_SIGMA");
+	if (!status) {
+	    gaussSigma = psMetadataLookupF32(&status, recipe, "MOMENTS_GAUSS_SIGMA");
+	}
+	float apScale = psMetadataLookupF32(&status, recipe, "PSF_APERTURE_SCALE");
+	float PSF_APERTURE = (int)(apScale*gaussSigma);
+	
+	float offset = pmGrowthCurveCorrect (psf->growth, PSF_APERTURE);
+	psMetadataAddF32 (readout->analysis, PS_LIST_TAIL, "APLOSS",   PS_META_REPLACE, "aperture loss (mag)", psf->growth->apLoss);
+	psMetadataAddF32 (readout->analysis, PS_LIST_TAIL, "APREFOFF", PS_META_REPLACE, "offset to ref aperture", offset);
+    } else {
+	psMetadataAddF32 (readout->analysis, PS_LIST_TAIL, "APLOSS",   PS_META_REPLACE, "aperture loss (mag)", NAN);
+	psMetadataAddF32 (readout->analysis, PS_LIST_TAIL, "APREFOFF", PS_META_REPLACE, "offset to ref aperture", NAN);
+    }
 
     psLogMsg ("psphot.apresid", PS_LOG_DETAIL, "aperture residual: %f +/- %f\n", psf->ApResid, psf->dApResid);
Index: /trunk/psphot/src/psphotBlendFit.c
===================================================================
--- /trunk/psphot/src/psphotBlendFit.c	(revision 31451)
+++ /trunk/psphot/src/psphotBlendFit.c	(revision 31452)
@@ -242,4 +242,7 @@
         if (source->mode &  PM_SOURCE_MODE_PAIR) continue;
 
+	// do not include MOMENTS_FAILURES in the fit
+        if (source->mode & PM_SOURCE_MODE_MOMENTS_FAILURE) continue;
+
         // limit selection to some SN limit
         if (sqrt(source->peak->detValue) < FIT_SN_LIM) continue;
@@ -280,8 +283,6 @@
         if (source->mode & PM_SOURCE_MODE_EXT_LIMIT) {
             if (psphotFitBlob (readout, source, newSources, psf, fitOptions, maskVal, markVal)) {
-                source->type = PM_SOURCE_TYPE_EXTENDED;
                 psTrace ("psphot", 5, "source at %7.1f, %7.1f is ext", source->peak->xf, source->peak->yf);
                 Next ++;
-                source->mode |= PM_SOURCE_MODE_NONLINEAR_FIT;
                 continue;
             }
@@ -291,5 +292,4 @@
                 psTrace ("psphot", 5, "source at %7.1f, %7.1f is psf", source->peak->xf, source->peak->yf);
                 Npsf ++;
-                source->mode |= PM_SOURCE_MODE_NONLINEAR_FIT;
                 continue;
             }
Index: /trunk/psphot/src/psphotChoosePSF.c
===================================================================
--- /trunk/psphot/src/psphotChoosePSF.c	(revision 31451)
+++ /trunk/psphot/src/psphotChoosePSF.c	(revision 31452)
@@ -350,5 +350,5 @@
 
     // build curve-of-growth vector for this psf
-    if (!psphotMakeGrowthCurve (readout, recipe, try->psf)) {
+    if (!psphotMakeGrowthCurve (readout, recipe, try->psf, try->sources)) {
         psError(PSPHOT_ERR_PSF, false, "Unable to construct flux scale for PSF");
         psFree (models);
Index: /trunk/psphot/src/psphotEfficiency.c
===================================================================
--- /trunk/psphot/src/psphotEfficiency.c	(revision 31451)
+++ /trunk/psphot/src/psphotEfficiency.c	(revision 31452)
@@ -400,4 +400,7 @@
                 source->type = PM_SOURCE_TYPE_STAR;
 
+		source->modelPSF->fitRadius = sourceRadius;
+		source->apRadius = sourceRadius;
+
                 numFound++;
                 psArrayAdd(sources, sources->n, source);
@@ -468,8 +471,8 @@
                     source->modelPSF->params->data.F32[PM_PAR_XPOS],
                     source->modelPSF->params->data.F32[PM_PAR_YPOS],
-                    magRef, source->psfMag, source->errMag);
+                    magRef, source->psfMag, source->psfMagErr);
 #endif
             magDiff->data.F32[j] = source->psfMag - magRef;
-            magErr->data.F32[j] = source->errMag;
+            magErr->data.F32[j] = source->psfMagErr;
         }
         magDiff->n = sources->n;
Index: /trunk/psphot/src/psphotExtendedSourceFits.c
===================================================================
--- /trunk/psphot/src/psphotExtendedSourceFits.c	(revision 31451)
+++ /trunk/psphot/src/psphotExtendedSourceFits.c	(revision 31452)
@@ -333,5 +333,5 @@
 	// this uses the footprint to judge both radius and aperture?
 	// XXX save the psf-based moments for output
-	if (!pmSourceMoments (source, radius, 0.0, 0.0, maskVal)) {
+	if (!pmSourceMoments (source, radius, 0.0, 0.0, 0.0, maskVal)) {
 	    fprintf (stderr, "skipping (2) %f, %f\n", source->peak->xf, source->peak->yf);
 	    // subtract the best fit from the object, leave local sky
@@ -379,8 +379,6 @@
           assert (status);
 
-	  // fprintf (stderr, "xfit: %f, %f : %f\n", source->peak->xf, source->peak->yf, source->peak->SN); 
-
           // limit selection to some SN limit
-          assert (source->peak); // how can a source not have a peak?
+          // assert (source->peak); // how can a source not have a peak?
           if (sqrt(source->peak->detValue) < SNlim) {
 	      Nfaint ++;
Index: /trunk/psphot/src/psphotFitSourcesLinear.c
===================================================================
--- /trunk/psphot/src/psphotFitSourcesLinear.c	(revision 31451)
+++ /trunk/psphot/src/psphotFitSourcesLinear.c	(revision 31452)
@@ -137,4 +137,7 @@
         if (source->mode & PM_SOURCE_MODE_CR_LIMIT) continue;
 
+	// do not include MOMENTS_FAILURES in the fit
+        if (source->mode & PM_SOURCE_MODE_MOMENTS_FAILURE) continue;
+
 	// XXX count saturated stars
         if (source->mode & PM_SOURCE_MODE_SATSTAR) {
@@ -164,23 +167,37 @@
 
 	// check the integral of the model : is it large enough?
+	// apply mask?
 	float modelSum = 0.0;
+	float maskedSum = 0.0;
 	for (int iy = 0; iy < source->modelFlux->numRows; iy++) {
 	    for (int ix = 0; ix < source->modelFlux->numCols; ix++) {
 		modelSum += source->modelFlux->data.F32[iy][ix];
+		if (source->maskObj->data.PS_TYPE_IMAGE_MASK_DATA[iy][ix] & maskVal) continue;
+		maskedSum += source->modelFlux->data.F32[iy][ix];
 	    }
 	}
 	if (modelSum < 0.5) continue; // skip sources with no model constraint (somewhat arbitrary limit)
-	// if (modelSum < 0.01) continue; // skip sources with no model constraint (somewhat arbitrary limit)
 	if (modelSum < 0.8) {
-	    fprintf (stderr, "low-sig model @ %f, %f (%f sum, %f peak)\n",
-		     source->peak->xf, source->peak->yf, modelSum, source->peak->rawFlux);
+	    fprintf (stderr, "low-sig model @ %f, %f (%f masked sum, %f sum, %f peak)\n",
+		     source->peak->xf, source->peak->yf, maskedSum, modelSum, source->peak->rawFlux);
 	}
-
-        pmModel *model = pmSourceGetModel (NULL, source);
+	if (maskedSum < 0.5) {
+	    fprintf (stderr, "worrying model @ %f, %f (%f masked sum, %f sum, %f peak)\n",
+		     source->peak->xf, source->peak->yf, maskedSum, modelSum, source->peak->rawFlux);
+	}
+
+	bool isPSF = false;
+        pmModel *model = pmSourceGetModel (&isPSF, source);
 
         psImageMaskPixels (source->maskObj, "AND", PS_NOT_IMAGE_MASK(markVal));
         psImageKeepCircle (source->maskObj, source->peak->x, source->peak->y, model->fitRadius, "OR", markVal);
 
-        source->mode |= PM_SOURCE_MODE_LINEAR_FIT;
+	// we call this function multiple times. for the first time, we have only PSF models for all objects
+	// the second time has extended sources.  If we ever fit the PSF model, we should raise this bit
+	source->mode |= PM_SOURCE_MODE_LINEAR_FIT;
+	if (isPSF) {
+	    source->mode |= PM_SOURCE_MODE_PSFMODEL;
+	}	    
+	
         psArrayAdd (fitSources, 100, source);
     }
@@ -260,4 +277,33 @@
         }
     }
+
+# if (0)
+    static int npass = 0;
+    char name[128];
+    FILE *f1 = NULL;
+    int fd = -1;
+
+    snprintf (name, 128, "sparse.Aij.%02d.dat", npass);
+    f1 = fopen (name, "w");
+    psAssert (f1, "failed to open file\n");
+    fd = fileno (f1);
+    p_psVectorPrint (fd, sparse->Aij, "Aij");
+    fclose (f1);
+
+    snprintf (name, 128, "sparse.Bfj.%02d.dat", npass);
+    f1 = fopen (name, "w");
+    psAssert (f1, "failed to open file\n");
+    fd = fileno (f1);
+    p_psVectorPrint (fd, sparse->Bfj, "Bfj");
+    fclose (f1);
+
+    snprintf (name, 128, "sparse.Qii.%02d.dat", npass);
+    f1 = fopen (name, "w");
+    psAssert (f1, "failed to open file\n");
+    fd = fileno (f1);
+    p_psVectorPrint (fd, sparse->Qii, "Qii");
+    fclose (f1);
+    npass ++;
+# endif
 
     psSparseResort (sparse);
Index: /trunk/psphot/src/psphotFitSourcesLinearStack.c
===================================================================
--- /trunk/psphot/src/psphotFitSourcesLinearStack.c	(revision 31451)
+++ /trunk/psphot/src/psphotFitSourcesLinearStack.c	(revision 31452)
@@ -69,4 +69,17 @@
             if (!pmSourceCacheModel (source, maskVal)) continue;
           }
+
+	  // check the integral of the model : is it large enough?
+	  float modelSum = 0.0;
+	  for (int iy = 0; iy < source->modelFlux->numRows; iy++) {
+	      for (int ix = 0; ix < source->modelFlux->numCols; ix++) {
+		  modelSum += source->modelFlux->data.F32[iy][ix];
+	      }
+	  }
+	  if (modelSum < 0.5) continue; // skip sources with no model constraint (somewhat arbitrary limit)
+	  if (modelSum < 0.8) {
+	      fprintf (stderr, "low-sig model @ %f, %f (%f sum, %f peak)\n",
+		       source->peak->xf, source->peak->yf, modelSum, source->peak->rawFlux);
+	  }
 
           source->mode |= PM_SOURCE_MODE_LINEAR_FIT;
Index: /trunk/psphot/src/psphotGuessModels.c
===================================================================
--- /trunk/psphot/src/psphotGuessModels.c	(revision 31451)
+++ /trunk/psphot/src/psphotGuessModels.c	(revision 31452)
@@ -160,8 +160,4 @@
         pmSource *source = sources->data[i];
 
-        if (source->mode & PM_SOURCE_MODE_MOMENTS_FAILURE) {
-	    fprintf (stderr, "moment failure\n");
-	}
-
         // this is used to mark sources for which the model is measured. We check later that
         // all are used.
@@ -169,4 +165,5 @@
 
         // skip non-astronomical objects (very likely defects)
+        if (source->mode & PM_SOURCE_MODE_MOMENTS_FAILURE) continue;
         if (source->type == PM_SOURCE_TYPE_DEFECT) continue;
         if (source->type == PM_SOURCE_TYPE_SATURATED) continue;
@@ -196,4 +193,5 @@
         }
 
+# if (0)
 	if (source->mode & PM_SOURCE_MODE_SATSTAR) {
 	    fprintf (stderr, "satstar: %f,%f vs %f,%f : %c\n", 
@@ -202,4 +200,5 @@
 		     (useMoments ? 'T' : 'F'));
 	}
+# endif
 
         // set PSF parameters for this model (apply 2D shape model to coordinates Xo, Yo)
Index: /trunk/psphot/src/psphotKronMasked.c
===================================================================
--- /trunk/psphot/src/psphotKronMasked.c	(revision 31452)
+++ /trunk/psphot/src/psphotKronMasked.c	(revision 31452)
@@ -0,0 +1,273 @@
+# include "psphotInternal.h"
+
+bool psphotKronMag (pmSource *source, float radius, psImageMaskType maskVal);
+
+bool psphotKronMasked (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 (!psphotKronMaskedReadout (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;
+}
+
+bool psphotKronMaskedReadout(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");
+    }
+
+    // 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;
+
+    // generate the mask image: increment counter for every source overlapping the pixel
+    psImage *kronMask = psImageAlloc (readout->image->numCols, readout->image->numRows, PS_TYPE_U8);
+    psImageInit (kronMask, 0);
+    int Nx = kronMask->numCols;
+    int Ny = kronMask->numRows;
+    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);
+        }
+
+	psphotKronMag (source, 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] ++;
+		}
+	    }
+	}
+    }		
+    // 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);
+
+    psFree (kronMask);
+    return true;
+}
+
+bool psphotKronMag (pmSource *source, float radius, 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 Sum = 0.0;
+    psF32 Var = 0.0;
+    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;
+
+# if (1)
+    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;
+# 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
+    // 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
+
+    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];
+
+	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++) {
+	    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;
+
+	    psF32 pDiff = *vPix;
+
+	    Sum += pDiff;
+
+	    // Kron Flux uses the 1st radial moment (NOT Gaussian windowed)
+	    psF32 rf = pDiff * sqrt(r2);
+	    psF32 rs = pDiff;
+
+	    RF  += rf;
+	    RS  += rs;
+	}
+    }
+
+    // Saturate the 1st radial moment
+    float Mrf = MIN(radius, MAX(0.25*radius, RF/RS));
+
+    // Calculate the Kron magnitude (make this block optional?)
+    float radKron  = 2.5*Mrf;
+    float radKron2 = radKron*radKron;
+
+    int nKronPix = 0;
+    Sum = Var = 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];
+
+	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++) {
+	    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;
+
+	    psF32 pDiff = *vPix;
+	    psF32 wDiff = *vWgt;
+
+	    Sum += pDiff;
+	    Var += wDiff;
+	    nKronPix ++;
+	}
+    }
+    // XXX for a test, save the old values here:
+    source->moments->KronCore    = 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;
+
+    return true;
+}
Index: /trunk/psphot/src/psphotLoadSRCTEXT.c
===================================================================
--- /trunk/psphot/src/psphotLoadSRCTEXT.c	(revision 31451)
+++ /trunk/psphot/src/psphotLoadSRCTEXT.c	(revision 31452)
@@ -52,5 +52,5 @@
 
 	    // NOTE: most of these values are irrelevant for loaded source positions
-	    source->seq       = 0;
+	    source->seq       = source->id;
 	    PAR[PM_PAR_XPOS]  = X;
 	    PAR[PM_PAR_YPOS]  = Y;
@@ -67,5 +67,5 @@
 
 	    source->psfMag    = 0.0;
-	    source->errMag    = 0.0;
+	    source->psfMagErr    = 0.0;
 	    source->apMag     = 0.0;
 
Index: /trunk/psphot/src/psphotMagnitudes.c
===================================================================
--- /trunk/psphot/src/psphotMagnitudes.c	(revision 31451)
+++ /trunk/psphot/src/psphotMagnitudes.c	(revision 31452)
@@ -167,4 +167,18 @@
         pmSource *source = (pmSource *) sources->data[i];
 
+	bool saveTest = false;
+	psImage *testImage = NULL;
+# if (0)
+	if ((fabs(source->peak->xf-3518) < 5) && (fabs(source->peak->yf-3178) < 5)) {
+	    saveTest = true;
+	    psRegion subregion = psRegionSet (source->peak->xf - 200, source->peak->xf + 200, source->peak->yf - 200, source->peak->yf + 200);
+	    testImage = psImageSubset ((psImage *) source->pixels->parent, subregion);
+	}
+# endif
+
+	if (saveTest) {
+	    psphotSaveImage(NULL, testImage, "test.image.1.fits");
+	}
+
         // replace object in image
         if (source->tmpFlags & PM_SOURCE_TMPF_SUBTRACTED) {
@@ -172,10 +186,33 @@
         }
 
+	if (saveTest) {
+	    psphotSaveImage(NULL, testImage, "test.image.2.fits");
+	}
+
+	float xPos = source->peak->xf;
+	float yPos = source->peak->yf;
+
+	pmModel *model = source->modelPSF;
+	if (model) {
+	    xPos = model->params->data.F32[PM_PAR_XPOS];
+	    yPos = model->params->data.F32[PM_PAR_YPOS];
+	} else {
+	    bool useMoments = pmSourcePositionUseMoments(source);
+	    if (useMoments) {
+		xPos = source->moments->Mx;
+		yPos = source->moments->My;
+	    }
+	}
+
         // clear the mask bit and set the circular mask pixels
         psImageMaskPixels (source->maskObj, "AND", PS_NOT_IMAGE_MASK(markVal));
-        psImageKeepCircle (source->maskObj, source->peak->x, source->peak->y, source->apRadius, "OR", markVal);
+        psImageKeepCircle (source->maskObj, xPos, yPos, source->apRadius, "OR", markVal);
 
         status = pmSourceMagnitudes (source, psf, photMode, maskVal, markVal, source->apRadius);
-        if (status && isfinite(source->apMag)) Nap ++;
+        if (status && isfinite(source->apFlux)) {
+	    Nap ++;
+	} else {
+	    fprintf (stderr, "failed to measure mag for source @ %f,%f\n", source->peak->xf, source->peak->yf);
+	}
 
         // clear the mask bit
@@ -185,7 +222,11 @@
         pmSourceSub (source, PM_MODEL_OP_FULL, maskVal);
 
+	if (saveTest) {
+	    psphotSaveImage(NULL, testImage, "test.image.3.fits");
+	}
+
         if (backModel) {
             psAssert (binning, "if backModel is defined, so should binning be");
-            source->sky = psImageUnbinPixel(source->peak->x, source->peak->y, backModel->image, binning);
+            source->sky = psImageUnbinPixel(xPos, yPos, backModel->image, binning);
             if (isnan(source->sky) && false) {
                 psLogMsg ("psphot.magnitudes", PS_LOG_DETAIL, "error setting pmSource.sky");
@@ -197,5 +238,5 @@
         if (backStdev) {
             psAssert (binning, "if backStdev is defined, so should binning be");
-            source->skyErr = psImageUnbinPixel(source->peak->x, source->peak->y, backStdev->image, binning);
+            source->skyErr = psImageUnbinPixel(xPos, yPos, backStdev->image, binning);
             if (isnan(source->skyErr) && false) {
                 psLogMsg ("psphot.magnitudes", PS_LOG_DETAIL, "error setting pmSource.skyErr");
Index: /trunk/psphot/src/psphotMakeGrowthCurve.c
===================================================================
--- /trunk/psphot/src/psphotMakeGrowthCurve.c	(revision 31451)
+++ /trunk/psphot/src/psphotMakeGrowthCurve.c	(revision 31452)
@@ -1,5 +1,5 @@
 # include "psphotInternal.h"
 
-bool psphotMakeGrowthCurve (pmReadout *readout, psMetadata *recipe, pmPSF *psf) {
+bool psphotMakeGrowthCurve (pmReadout *readout, psMetadata *recipe, pmPSF *psf, psArray *sources) {
 
     bool status;
@@ -22,18 +22,41 @@
 
     // measure the aperture loss as a function of radius for PSF
-    bool IGNORE_GROWTH = psMetadataLookupBool (&status, recipe, "IGNORE_GROWTH");
     float REF_RADIUS = psMetadataLookupF32 (&status, recipe, "PSF_REF_RADIUS");
     float PSF_FIT_PAD   = psMetadataLookupF32 (&status, recipe, "PSF_FIT_PADDING");
+    
+    float gaussSigma = psMetadataLookupF32(&status, readout->analysis, "MOMENTS_GAUSS_SIGMA");
+    if (!status) {
+	gaussSigma = psMetadataLookupF32(&status, recipe, "MOMENTS_GAUSS_SIGMA");
+    }
+    float apScale = psMetadataLookupF32(&status, recipe, "PSF_APERTURE_SCALE");
+    float PSF_APERTURE = (int)(apScale*gaussSigma);
 
     psf->growth = pmGrowthCurveAlloc (PSF_FIT_PAD, 100.0, REF_RADIUS);
 
-    if (!pmGrowthCurveGenerate (readout, psf, IGNORE_GROWTH, maskVal, markVal)) {
-        psError(PSPHOT_ERR_APERTURE, false, "Fitting aperture corrections");
-        psFree(psf->growth); psf->growth = NULL;
-        return false;
+    bool GROWTH_FROM_SOURCES = psMetadataLookupBool (&status, recipe, "GROWTH_FROM_SOURCES");
+    if (!status) GROWTH_FROM_SOURCES = false;
+
+    if (GROWTH_FROM_SOURCES) {
+	bool INTERPOLATE_AP = psMetadataLookupBool (&status, recipe, "INTERPOLATE_AP");
+	if (!pmGrowthCurveGenerateFromSources (readout, psf, sources, INTERPOLATE_AP, maskVal, markVal)) {
+	    psError(PSPHOT_ERR_APERTURE, false, "Fitting aperture corrections");
+	    psFree(psf->growth); psf->growth = NULL;
+	    return false;
+	}
+
+    } else {
+	bool IGNORE_GROWTH = psMetadataLookupBool (&status, recipe, "IGNORE_GROWTH");
+	if (!pmGrowthCurveGenerate (readout, psf, IGNORE_GROWTH, maskVal, markVal)) {
+	    psError(PSPHOT_ERR_APERTURE, false, "Fitting aperture corrections");
+	    psFree(psf->growth); psf->growth = NULL;
+	    return false;
+	}
     }
 
     psLogMsg ("psphot", PS_LOG_MINUTIA, "built growth curve: %f sec\n", psTimerMark ("psphot.growth"));
 
+    float offset = pmGrowthCurveCorrect (psf->growth, PSF_APERTURE);
+    psLogMsg ("psphot", PS_LOG_DETAIL, "correction from %f to %f: %f mags\n", PSF_APERTURE, REF_RADIUS, offset);
+
     return true;
 }
Index: /trunk/psphot/src/psphotMergeSources.c
===================================================================
--- /trunk/psphot/src/psphotMergeSources.c	(revision 31451)
+++ /trunk/psphot/src/psphotMergeSources.c	(revision 31452)
@@ -137,4 +137,9 @@
                 pmSource *source = extSourcesTXT->data[i];
                 source->mode |= PM_SOURCE_MODE_EXTERNAL;
+
+                // the supplied peak flux needs to be re-normalized
+                source->peak->rawFlux = 1.0;
+                source->peak->smoothFlux = 1.0;
+                source->peak->detValue = 1.0;
 
                 // drop the loaded source modelPSF
Index: /trunk/psphot/src/psphotModelTest.c
===================================================================
--- /trunk/psphot/src/psphotModelTest.c	(revision 31451)
+++ /trunk/psphot/src/psphotModelTest.c	(revision 31452)
@@ -226,5 +226,5 @@
     // measure the source mags
     pmSourcePhotometryModel (&fitMag, model);
-    pmSourcePhotometryAper  (&obsMag, NULL, NULL, model, source->pixels, source->variance, source->maskObj, maskVal);
+    pmSourcePhotometryAper  (NULL, &obsMag, NULL, NULL, model, source->pixels, source->variance, source->maskObj, maskVal);
     fprintf (stderr, "ap: %f, fit: %f, apmifit: %f, nIter: %d\n", obsMag, fitMag, obsMag - fitMag, model->nIter);
 
Index: /trunk/psphot/src/psphotOutput.c
===================================================================
--- /trunk/psphot/src/psphotOutput.c	(revision 31451)
+++ /trunk/psphot/src/psphotOutput.c	(revision 31452)
@@ -75,5 +75,5 @@
 		 source->peak->xf, source->peak->yf, 
 		 model->params->data.F32[PM_PAR_XPOS], model->params->data.F32[PM_PAR_YPOS], 
-		 source->psfMag, source->apMag, source->errMag,
+		 source->psfMag, source->apMag, source->psfMagErr,
 		 model->params->data.F32[PM_PAR_I0], 
 		 model->params->data.F32[PM_PAR_SXX], model->params->data.F32[PM_PAR_SXY], model->params->data.F32[PM_PAR_SYY], 
@@ -251,4 +251,5 @@
     psMetadataItemSupplement (&status, header, analysis, "NPSFSTAR");
     psMetadataItemSupplement (&status, header, analysis, "APLOSS");
+    psMetadataItemSupplement (&status, header, analysis, "APREFOFF");
 
     psMetadataItemSupplement (&status, header, analysis, "FWHM_MAJ");
Index: /trunk/psphot/src/psphotPetrosianStats.c
===================================================================
--- /trunk/psphot/src/psphotPetrosianStats.c	(revision 31451)
+++ /trunk/psphot/src/psphotPetrosianStats.c	(revision 31452)
@@ -163,5 +163,5 @@
 	if (refRadius->data.F32[i] > apRadius) {
 	    if (i == 0) {
-		psWarning ("does this case make any sense? (refRadius[0] > apRadius)");
+		psWarning ("does this case make any sense? (refRadius[0] %f > apRadius %f)", refRadius->data.F32[i], apRadius);
 		continue;
 	    } else {
@@ -188,5 +188,5 @@
 	if (!found50 && (fluxSum->data.F32[i] > flux50)) {
 	    if (i == 0) {
-		psWarning ("does this case make any sense? (fluxSum[0] > flux50)");
+		psWarning ("does this case make any sense? (fluxSum[0] %f > flux50 %f)", fluxSum->data.F32[i], flux50);
 		continue;
 	    } else {
@@ -198,5 +198,5 @@
 	if (!found90 && (fluxSum->data.F32[i] > flux90)) {
 	    if (i == 0) {
-		psWarning ("does this case make any sense? (fluxSum[0] > flux90)");
+		psWarning ("does this case make any sense? (fluxSum[0] %f > flux90 %f)", fluxSum->data.F32[i], flux90);
 		continue;
 	    } else {
Index: /trunk/psphot/src/psphotRadialApertures.c
===================================================================
--- /trunk/psphot/src/psphotRadialApertures.c	(revision 31451)
+++ /trunk/psphot/src/psphotRadialApertures.c	(revision 31452)
@@ -95,5 +95,7 @@
 	if (source->type == PM_SOURCE_TYPE_SATURATED) continue;
 	if (source->mode & PM_SOURCE_MODE_DEFECT) continue;
-	if (source->mode & PM_SOURCE_MODE_SATSTAR) continue;
+
+	// XXX measure radial apertures even for saturated stars
+	// if (source->mode & PM_SOURCE_MODE_SATSTAR) continue;
 
 	// limit selection to some SN limit
@@ -257,9 +259,12 @@
 	float SBstdv = sqrt((fluxStd->data.F32[i] / nPix) - PS_SQR(SBmean));
 
+	// XXX report the total flux or the mask-corrected flux?
 	// flux->data.F32[i]    = SBmean * Area;
+	// fluxErr->data.F32[i] = sqrt(fluxErr->data.F32[i]) * Area / nPix;
+
+	fluxErr->data.F32[i] = sqrt(fluxErr->data.F32[i]);
 	fluxStd->data.F32[i] = SBstdv * Area;
-	fluxErr->data.F32[i] = sqrt(fluxErr->data.F32[i]) * Area / nPix;
-
-	// fill->data.F32[i] /= Area;
+	fill->data.F32[i] /= Area;
+
 	psTrace ("psphot", 5, "radial bins: %3d  %5.1f : %8.1f +/- %7.2f : %8.1f +/- %8.1f : %4.2f %6.1f\n", 
 		 i, aperRadii->data.F32[i], flux->data.F32[i], fluxErr->data.F32[i], SBmean, SBstdv, fill->data.F32[i], Area);
Index: /trunk/psphot/src/psphotRadiusChecks.c
===================================================================
--- /trunk/psphot/src/psphotRadiusChecks.c	(revision 31451)
+++ /trunk/psphot/src/psphotRadiusChecks.c	(revision 31452)
@@ -184,8 +184,5 @@
 bool psphotSetRadiusModel (pmModel *model, pmReadout *readout, pmSource *source, psImageMaskType markVal, bool deep) {
 
-    psF32 *PAR = model->params->data.F32;
-
-    pmMoments *moments = source->moments;
-    if (moments == NULL) return false;
+    pmPeak *peak = source->peak;
 
     // set the fit radius based on the object flux limit and the model
@@ -199,8 +196,8 @@
 
     // redefine the pixels if needed
-    pmSourceRedefinePixels (source, readout, PAR[PM_PAR_XPOS], PAR[PM_PAR_YPOS], model->fitRadius);
-
-    // set the mask to flag the excluded pixels
-    psImageKeepCircle (source->maskObj, PAR[PM_PAR_XPOS], PAR[PM_PAR_YPOS], model->fitRadius, "OR", markVal);
-    return true;
-}
+    pmSourceRedefinePixels (source, readout, peak->xf, peak->yf, model->fitRadius);
+
+    // set the mask to flag the excluded pixels
+    psImageKeepCircle (source->maskObj, peak->xf, peak->yf, model->fitRadius, "OR", markVal);
+    return true;
+}
Index: /trunk/psphot/src/psphotReadout.c
===================================================================
--- /trunk/psphot/src/psphotReadout.c	(revision 31451)
+++ /trunk/psphot/src/psphotReadout.c	(revision 31452)
@@ -200,4 +200,7 @@
     psphotDumpChisqs (config, view, filerule);
 
+    // XXX re-measure the kron mags with models subtracted
+    psphotKronMasked(config, view, filerule);
+
     // identify CRs and extended sources (only unmeasured sources are measured)
     psphotSourceSize (config, view, filerule, true); // pass 1 (detections->allSources)
@@ -309,4 +312,7 @@
 pass1finish:
 
+    // XXX re-measure the kron mags with models subtracted
+    psphotKronMasked(config, view, filerule);
+
     // measure source size for the remaining sources
     // NOTE: applies only to NEW (unmeasured) sources
@@ -345,16 +351,12 @@
       psErrorStackPrint(stderr, "Unable to replace sky");
       psErrorClear();
-
-/*       psLogMsg("psphot", 3, "failed on psphotSkyReplace"); */
-/*       return(psphotReadoutCleanup(config, view, filerule)); */
-    }
+    }
+
     // drop the references to the image pixels held by each source
     if (!psphotSourceFreePixels (config, view, filerule)) { // pass 1
       psErrorStackPrint(stderr, "Unable to free source pixels");
       psErrorClear();
-
-/*       psLogMsg ("psphot", 3, "failed on psphotSourceFreePixels"); */
-/*       return(psphotReadoutCleanup(config, view, filerule)); */
-    }
+    }
+
     // create the exported-metadata and free local data
     return psphotReadoutCleanup(config, view, filerule);
Index: /trunk/psphot/src/psphotReplaceUnfit.c
===================================================================
--- /trunk/psphot/src/psphotReplaceUnfit.c	(revision 31451)
+++ /trunk/psphot/src/psphotReplaceUnfit.c	(revision 31452)
@@ -67,4 +67,11 @@
     psAssert (maskVal, "missing mask value?");
 
+    // 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;
+
     for (int i = 0; i < sources->n; i++) {
       source = sources->data[i];
Index: /trunk/psphot/src/psphotRoughClass.c
===================================================================
--- /trunk/psphot/src/psphotRoughClass.c	(revision 31451)
+++ /trunk/psphot/src/psphotRoughClass.c	(revision 31452)
@@ -207,5 +207,9 @@
         return false;
     }
-    psLogMsg ("psphot", 3, "psf clump  X,  Y: %f, %f : DX, DY: %f, %f : nStars %d of %d\n", psfClump.X, psfClump.Y, psfClump.dX, psfClump.dY, psfClump.nStars, psfClump.nTotal);
+    if (!havePSF) {
+	psLogMsg ("psphot", 3, "psf clump  X,  Y: %f, %f : DX, DY: %f, %f : nStars %d of %d\n", psfClump.X, psfClump.Y, psfClump.dX, psfClump.dY, psfClump.nStars, psfClump.nTotal);
+    } else {
+	psLogMsg ("psphot", 3, "psf clump  X,  Y: %f, %f : DX, DY: %f, %f : loaded from metadata\n", psfClump.X, psfClump.Y, psfClump.dX, psfClump.dY);
+    }
 
     // get basic parameters, or set defaults
Index: /trunk/psphot/src/psphotSetThreads.c
===================================================================
--- /trunk/psphot/src/psphotSetThreads.c	(revision 31451)
+++ /trunk/psphot/src/psphotSetThreads.c	(revision 31452)
@@ -25,5 +25,5 @@
     psFree(task);
 
-    task = psThreadTaskAlloc("PSPHOT_SOURCE_STATS", 10);
+    task = psThreadTaskAlloc("PSPHOT_SOURCE_STATS", 11);
     task->function = &psphotSourceStats_Threaded;
     psThreadTaskAdd(task);
Index: /trunk/psphot/src/psphotSourceFits.c
===================================================================
--- /trunk/psphot/src/psphotSourceFits.c	(revision 31451)
+++ /trunk/psphot/src/psphotSourceFits.c	(revision 31452)
@@ -146,4 +146,5 @@
     pmSourceSub (source, PM_MODEL_OP_FULL, maskVal);
     source->mode |=  PM_SOURCE_MODE_BLEND_FIT;
+    source->mode |= PM_SOURCE_MODE_NONLINEAR_FIT;
     return true;
 }
@@ -186,4 +187,5 @@
 
     // build cached model and subtract
+    source->mode |= PM_SOURCE_MODE_NONLINEAR_FIT;
     pmSourceCacheModel (source, maskVal);
     pmSourceSub (source, PM_MODEL_OP_FULL, maskVal);
@@ -233,5 +235,5 @@
     // XXX save the psf-based moments for output
     psfMoments = *source->moments;
-    if (!pmSourceMoments (source, radius, 0.0, 0.5, maskVal)) {
+    if (!pmSourceMoments (source, radius, 0.0, 0.5, 0.0, maskVal)) {
       *source->moments = psfMoments;
       return false;
@@ -319,4 +321,5 @@
     source->type = PM_SOURCE_TYPE_EXTENDED;
     source->mode |= PM_SOURCE_MODE_EXTMODEL;
+    source->mode |= PM_SOURCE_MODE_NONLINEAR_FIT;
 
     // build cached model and subtract
@@ -346,4 +349,5 @@
     source->modelPSF = psMemIncrRefCounter (DBL->data[0]);
     source->mode     |= PM_SOURCE_MODE_PAIR;
+    source->mode |= PM_SOURCE_MODE_NONLINEAR_FIT;
 
     // copy most data from the primary source (modelEXT, blends stay NULL)
@@ -489,5 +493,5 @@
     // psTraceSetLevel("psLib.math.psMinimizeLMChi2", 5);
     pmSourceFitModel (source, model, &options, maskVal);
-    // fprintf (stderr, "chisq: %f, nIter: %d, radius: %f, npix: %d\n\n", model->chisqNorm, model->nIter, model->fitRadius, model->nPix);
+    // fprintf (stderr, "chisq: %f, nIter: %d, radius: %f, npix: %d\n", model->chisqNorm, model->nIter, model->fitRadius, model->nPix);
 
     // psTraceSetLevel("psLib.math.psMinimizeLMChi2", 0);
@@ -554,5 +558,5 @@
     // psTraceSetLevel("psLib.math.psMinimizeLMChi2", 5);
     pmSourceFitPCM (pcm, source, &options, maskVal, markVal, psfSize);
-    // fprintf (stderr, "chisq: %f, nIter: %d, radius: %f, npix: %d\n\n", model->chisqNorm, model->nIter, model->fitRadius, model->nPix);
+    // fprintf (stderr, "chisq: %f, nIter: %d, radius: %f, npix: %d\n", model->chisqNorm, model->nIter, model->fitRadius, model->nPix);
 
     // psTraceSetLevel("psLib.math.psMinimizeLMChi2", 0);
Index: /trunk/psphot/src/psphotSourceSize.c
===================================================================
--- /trunk/psphot/src/psphotSourceSize.c	(revision 31451)
+++ /trunk/psphot/src/psphotSourceSize.c	(revision 31452)
@@ -213,5 +213,5 @@
 
         psVectorAppend (ApOff, dMag);
-        psVectorAppend (ApErr, source->errMag);
+        psVectorAppend (ApErr, source->psfMagErr);
     }
     if (num == 0) {
@@ -520,5 +520,5 @@
         // set nSigmaMAG to include both systematic and poisson error terms.  we include a hard
 	// floor on the Ap Sys Err (to be a bit generous).  XXX put the floor in the recipe...
-        float nSigmaMAG = (dMag - options->ApResid) / hypot(source->errMag, hypot(options->ApSysErr, 0.02));
+        float nSigmaMAG = (dMag - options->ApResid) / hypot(source->psfMagErr, hypot(options->ApSysErr, 0.02));
         source->extNsigma = nSigmaMAG;
 
Index: /trunk/psphot/src/psphotSourceStats.c
===================================================================
--- /trunk/psphot/src/psphotSourceStats.c	(revision 31451)
+++ /trunk/psphot/src/psphotSourceStats.c	(revision 31452)
@@ -98,6 +98,8 @@
 
     // generate the array of sources, define the associated pixel
+    bool firstPass = false;
     if (!detections->newSources) {
         detections->newSources = psArrayAllocEmpty (peaks->n);
+	firstPass = true;
     }
     sources = detections->newSources;
@@ -125,4 +127,8 @@
         if (source->mode & PM_SOURCE_MODE_MOMENTS_FAILURE) {
 	    fprintf (stderr, "moment failure\n");
+	}
+
+	if (firstPass) {
+            source->mode2 |= PM_SOURCE_MODE2_PASS1_SRC;
 	}
 
@@ -159,4 +165,8 @@
         RADIUS = psMetadataLookupF32 (&status, recipe, "PSF_MOMENTS_RADIUS");
     }
+    float MIN_KRON_RADIUS = psMetadataLookupF32 (&status, readout->analysis, "MOMENTS_MIN_KRON");
+    if (!status) {
+        MIN_KRON_RADIUS = 0.75*SIGMA;
+    }
 
     // threaded measurement of the source magnitudes
@@ -186,4 +196,5 @@
             PS_ARRAY_ADD_SCALAR(job->args, RADIUS,  PS_TYPE_F32);
             PS_ARRAY_ADD_SCALAR(job->args, SIGMA,   PS_TYPE_F32);
+            PS_ARRAY_ADD_SCALAR(job->args, MIN_KRON_RADIUS, PS_TYPE_F32);
 
             PS_ARRAY_ADD_SCALAR(job->args, maskVal, PS_TYPE_IMAGE_MASK);
@@ -215,9 +226,9 @@
             } else {
                 psScalar *scalar = NULL;
-                scalar = job->args->data[7];
+                scalar = job->args->data[8];
                 Nmoments += scalar->data.S32;
-                scalar = job->args->data[8];
+                scalar = job->args->data[9];
                 Nfail += scalar->data.S32;
-                scalar = job->args->data[9];
+                scalar = job->args->data[10];
                 Nfaint += scalar->data.S32;
             }
@@ -231,4 +242,7 @@
 
     psphotVisualShowMoments (sources);
+
+    // clear the mark bits
+    // psImageMaskPixels (readout->mask, "AND", PS_NOT_IMAGE_MASK(markVal));
 
     if (detections->allSources) {
@@ -373,7 +387,8 @@
     float RADIUS                    = PS_SCALAR_VALUE(job->args->data[3],F32);
     float SIGMA                     = PS_SCALAR_VALUE(job->args->data[4],F32);
-
-    psImageMaskType maskVal         = PS_SCALAR_VALUE(job->args->data[5],PS_TYPE_IMAGE_MASK_DATA);
-    psImageMaskType markVal         = PS_SCALAR_VALUE(job->args->data[6],PS_TYPE_IMAGE_MASK_DATA);
+    float MIN_KRON_RADIUS           = PS_SCALAR_VALUE(job->args->data[5],F32);
+
+    psImageMaskType maskVal         = PS_SCALAR_VALUE(job->args->data[6],PS_TYPE_IMAGE_MASK_DATA);
+    psImageMaskType markVal         = PS_SCALAR_VALUE(job->args->data[7],PS_TYPE_IMAGE_MASK_DATA);
 
     // if no valid pixels, massive swing or very large Mrf, likely saturated source,
@@ -431,5 +446,5 @@
 	if (!(source->peak->type == PM_PEAK_SUSPECT_SATURATION)) {
 	    // measure basic source moments (no S/N clipping on input pixels)
-	    status = pmSourceMoments (source, RADIUS, SIGMA, 0.0, maskVal);
+	    status = pmSourceMoments (source, RADIUS, SIGMA, 0.0, MIN_KRON_RADIUS, maskVal);
 	} else { 
 	    // For saturated stars, choose a much larger box NOTE this is slightly sleazy, but
@@ -447,5 +462,5 @@
 
 	    psTrace ("psphot", 4, "retrying moments for %d, %d\n", source->peak->x, source->peak->y);
-	    status = pmSourceMoments (source, BIG_RADIUS, BIG_SIGMA, 0.0, maskVal);
+	    status = pmSourceMoments (source, BIG_RADIUS, BIG_SIGMA, 0.0, MIN_KRON_RADIUS, maskVal);
 	    source->mode |= PM_SOURCE_MODE_BIG_RADIUS;
 	}
@@ -456,4 +471,10 @@
 	    continue;
 	}
+
+	// XXX test of masking neighbors when measureing moments (does this fail?)
+	// psEllipseMoments moments = {source->moments->Mxx, source->moments->Myy, source->moments->Mxy};
+	// psEllipseAxes axes = psEllipseMomentsToAxes(moments, 20.0);
+        // psImageMaskCircle (source->maskView, source->peak->x, source->peak->y, 3.0*axes.major, "OR", markVal);
+
 	Nmoments ++;
 
@@ -492,4 +513,5 @@
     float sigma[NSIGMA] = {1.0, 2.0, 3.0, 4.5, 6.0, 9.0, 12.0, 18.0};
     float Sout[NSIGMA];
+    float Rmin[NSIGMA];
     int   Nout[NSIGMA]; // number of stars found in clump : use this to control the number of regions measured by psphotRoughClass
 
@@ -513,5 +535,5 @@
 
             // measure basic source moments (no S/N clipping on input pixels)
-            status = pmSourceMoments (source, 4*sigma[i], sigma[i], 0.0, maskVal);
+            status = pmSourceMoments (source, 4*sigma[i], sigma[i], 0.0, 0.0, maskVal);
         }
 
@@ -524,4 +546,6 @@
         pmPSFClump psfClump = pmSourcePSFClump (NULL, NULL, sources, PSF_SN_LIM, PSF_CLUMP_GRID_SCALE, MOMENTS_SX_MAX, MOMENTS_SY_MAX, 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]);
+
+	Rmin[i] = pmSourceMinKronRadius(sources, PSF_SN_LIM);
 
 #if 0
@@ -549,4 +573,5 @@
     // we are looking for sigma for which Sout = 0.65 (or some other value)
     int Nstars = 0;
+    float minKronRadius = NAN;
     float Sigma = NAN;
     float minS = Sout[0];
@@ -565,4 +590,5 @@
         Sigma = sigma[i] + (0.65 - Sout[i])*(sigma[i+1] - sigma[i])/(Sout[i+1] - Sout[i]);
 	Nstars = 0.5*(Nout[i] + Nout[i+1]);
+        minKronRadius = Rmin[i] + (0.65 - Sout[i])*(Rmin[i+1] - Rmin[i])/(Sout[i+1] - Sout[i]);
     }
     psAssert (isfinite(Sigma), "did we miss a case?");
@@ -575,6 +601,7 @@
     psMetadataAddF32(analysis, PS_LIST_TAIL, "PSF_MOMENTS_RADIUS", PS_META_REPLACE, "moments limit", 4.0*Sigma);
     psMetadataAddF32(analysis, PS_LIST_TAIL, "PSF_CLUMP_NSTARS", PS_META_REPLACE, "number of stars in clump", Nstars);
-
-    psLogMsg ("psphot", 3, "using window function with sigma = %f\n", Sigma);
+    psMetadataAddF32(analysis, PS_LIST_TAIL, "MOMENTS_MIN_KRON", PS_META_REPLACE, "minimum Kron Radius", minKronRadius);
+
+    psLogMsg ("psphot", 3, "using window function with sigma = %f (Min Kron Radius = %f, Nstars = %d)\n", Sigma, minKronRadius, Nstars);
     return true;
 }
Index: /trunk/psphot/src/psphotVisual.c
===================================================================
--- /trunk/psphot/src/psphotVisual.c	(revision 31451)
+++ /trunk/psphot/src/psphotVisual.c	(revision 31452)
@@ -1867,5 +1867,5 @@
 	    continue;
 	}
-	if (source->errMag > 0.1) {
+	if (source->psfMagErr > 0.1) {
 	    xLOW->data.F32[nLOW] = source->moments->Mxx;
 	    yLOW->data.F32[nLOW] = source->moments->Myy;
@@ -2591,5 +2591,5 @@
 	x->data.F32[n] = source->psfMag;
 	y->data.F32[n] = dMag;
-	dy->data.F32[n] = source->errMag;
+	dy->data.F32[n] = source->psfMagErr;
 	graphdata.xmin = PS_MIN(graphdata.xmin, x->data.F32[n]);
 	graphdata.xmax = PS_MAX(graphdata.xmax, x->data.F32[n]);
Index: /trunk/psphot/test/tap_psphot_stackphot.pro
===================================================================
--- /trunk/psphot/test/tap_psphot_stackphot.pro	(revision 31451)
+++ /trunk/psphot/test/tap_psphot_stackphot.pro	(revision 31452)
@@ -1,4 +1,6 @@
 #!/usr/bin/env mana
 # -*-sh-*-
+
+$KAPA = kapa -noX
 
 # config for ppImage to generate chip, mask, weight