Index: trunk/psphot/src/Makefile.am
===================================================================
--- trunk/psphot/src/Makefile.am	(revision 33994)
+++ trunk/psphot/src/Makefile.am	(revision 34086)
@@ -183,5 +183,4 @@
 	psphotSourcePlots.c	       \
 	psphotRadialPlot.c	       \
-	psphotKronMasked.c             \
 	psphotKronIterate.c            \
 	psphotRadialProfileWings.c     \
@@ -209,8 +208,9 @@
 	psphotSetNFrames.c
 
-# re-instate these
+# not currently used
 #	psphotIsophotal.c	       \
 #	psphotAnnuli.c		       \
 #	psphotKron.c		       \
+#	psphotKronMasked.c             \
 #
 
Index: trunk/psphot/src/psphot.h
===================================================================
--- trunk/psphot/src/psphot.h	(revision 33994)
+++ trunk/psphot/src/psphot.h	(revision 34086)
@@ -61,4 +61,7 @@
 bool            psphotSetMaskAndVarianceReadout (pmConfig *config, const pmFPAview *view, const char *filerule, int index, psMetadata *recipe);
 
+bool            psphotUpdateVariance (pmConfig *config, const pmFPAview *view, const char *filerule);
+bool            psphotUpdateVarianceReadout (pmConfig *config, const pmFPAview *view, const char *filerule, int index, psMetadata *recipe);
+
 bool            psphotModelBackground (pmConfig *config, const pmFPAview *view, const char *filerule);
 bool            psphotModelBackgroundReadoutFileIndex (pmConfig *config, const pmFPAview *view, const char *filerule, int index);
@@ -99,10 +102,5 @@
 
 bool            psphotFitSourcesLinear (pmConfig *config, const pmFPAview *view, const char *filerule, bool final);
-
-# if (HAVE_MODEL_VAR)
 bool            psphotFitSourcesLinearReadout (psMetadata *recipe, pmReadout *readout, psArray *sources, pmPSF *psf, bool final, pmSourceFitVarMode fitVarMode);
-# else
-bool            psphotFitSourcesLinearReadout (psMetadata *recipe, pmReadout *readout, psArray *sources, pmPSF *psf, bool final);
-# endif
 
 bool            psphotSourceSize (pmConfig *config, const pmFPAview *view, const char *filerule, bool getPSFsize);
@@ -489,3 +487,16 @@
 bool psphotSetNFramesReadout (pmConfig *config, const pmFPAview *view, const char *filerule, int index);
 
+bool psphotGenerateModelVariance (pmConfig *config, const pmFPAview *view, pmFPAfile *file, int index, psMetadata *recipe, pmReadout *readout, psArray *sources);
+pmSourceFitVarMode psphotGetFitVarMode (psMetadata *recipe);
+bool psphotFreeModelVariance (pmReadout *readout, psArray *sources);
+
+bool psphotModelBackgroundReadout(psImage *model,  // Model image
+				  psImage *modelStdev, // Model stdev image
+				  psMetadata *analysis, // Analysis metadata for outputs
+				  pmReadout *readout, // Readout for which to generate a background model
+				  psImageBinning *binning, // Binning parameters
+				  const pmConfig *config,// Configuration
+				  bool useVarianceImage
+    );
+
 #endif
Index: trunk/psphot/src/psphotCleanup.c
===================================================================
--- trunk/psphot/src/psphotCleanup.c	(revision 33994)
+++ trunk/psphot/src/psphotCleanup.c	(revision 34086)
@@ -29,5 +29,13 @@
 psExit psphotGetExitStatus (void) {
 
-    psErrorCode err = psErrorCodeLast ();
+    // gcc -Wswitch complains here if err is declared as type psErrorCode
+    // the collection of ps*ErrorCode values are enums defined separately for 
+    // each module (psphot, pswarp, etc).  the lowest type, psErrorCode is only the base set and does
+    // not include the possible psphot values
+
+    // for now, to get around this, we just use an int for the switch
+
+    // psErrorCode err = psErrorCodeLast ();
+    int err = psErrorCodeLast ();
     switch (err) {
       case PS_ERR_NONE:
Index: trunk/psphot/src/psphotEfficiency.c
===================================================================
--- trunk/psphot/src/psphotEfficiency.c	(revision 33994)
+++ trunk/psphot/src/psphotEfficiency.c	(revision 34086)
@@ -420,9 +420,5 @@
 
     // psphotFitSourcesLinearReadout subtracts the model fits
-# if (HAVE_MODEL_VAR)
     if (!psphotFitSourcesLinearReadout(recipe, readout, fakeSourcesAll, psf, true, PM_SOURCE_PHOTFIT_CONST)) {
-# else
-    if (!psphotFitSourcesLinearReadout(recipe, readout, fakeSourcesAll, psf, true)) {
-# endif
         psError(PS_ERR_UNKNOWN, false, "Unable to perform linear fit on fake sources.");
         psFree(fakeSources);
Index: trunk/psphot/src/psphotFitSourcesLinear.c
===================================================================
--- trunk/psphot/src/psphotFitSourcesLinear.c	(revision 33994)
+++ trunk/psphot/src/psphotFitSourcesLinear.c	(revision 34086)
@@ -12,10 +12,4 @@
 static bool SetBorderMatrixElements (psSparseBorder *border, pmReadout *readout, psArray *sources, bool constant_weights, int SKY_FIT_ORDER, psImageMaskType markVal);
 
-# if (HAVE_MODEL_VAR)
-bool psphotGenerateModelVariance (pmConfig *config, const pmFPAview *view, pmFPAfile *file, int index, psMetadata *recipe, pmReadout *readout, psArray *sources);
-pmSourceFitVarMode psphotGetFitVarMode (psMetadata *recipe);
-bool psphotFreeModelVariance (pmReadout *readout, psArray *sources);
-# endif
-
 // for now, let's store the detections on the readout->analysis for each readout
 bool psphotFitSourcesLinear (pmConfig *config, const pmFPAview *view, const char *filerule, bool final)
@@ -30,5 +24,4 @@
     assert (recipe);
 
-# if (HAVE_MODEL_VAR)
     pmSourceFitVarMode fitVarMode = psphotGetFitVarMode (recipe);
     if (!fitVarMode) {
@@ -40,5 +33,4 @@
     // do a single pass.
     pmSourceFitVarMode fitVarModePass1 = (fitVarMode == PM_SOURCE_PHOTFIT_MODEL_VAR) ? PM_SOURCE_PHOTFIT_CONST : fitVarMode;
-# endif
 
     int num = psphotFileruleCount(config, filerule);
@@ -68,15 +60,9 @@
         psAssert (psf, "missing psf?");
 
-# if (HAVE_MODEL_VAR)
-        if (!psphotFitSourcesLinearReadout (recipe, readout, sources, psf, final, fitVarModePass1))
-# else
-        if (!psphotFitSourcesLinearReadout (recipe, readout, sources, psf, final))
-# endif
-	{
+        if (!psphotFitSourcesLinearReadout (recipe, readout, sources, psf, final, fitVarModePass1)) {
             psError (PSPHOT_ERR_CONFIG, false, "failed to fit sources (linear) for %s entry %d", filerule, i);
             return false;
         }
 
-# if (HAVE_MODEL_VAR)
 	// the MODEL_VAR weighting scheme requires knowledge of the model fluxes to generate the variance
 	// after we have determined the initial set of fits, then we can generate the variance image and 
@@ -104,5 +90,4 @@
 	    }
 	}
-# endif
 
 	psphotVisualShowResidualImage (readout, (num > 0)); 
@@ -113,5 +98,4 @@
 }
 
-# if (HAVE_MODEL_VAR)
 // look up the fit variance mode from the recipe; older recipes do not have the value
 // 'LINEAR_FIT_VARIANCE_MODE'; in those cases, look for 'CONSTANT_PHOTOMETRIC_WEIGHTS' as a boolean and
@@ -136,4 +120,7 @@
 	return PM_SOURCE_PHOTFIT_IMAGE_VAR;
     }
+    if (!strcasecmp(fitVarModeString, "SKY")   || !strcasecmp(fitVarModeString, "MODEL_SKY")) {
+	return PM_SOURCE_PHOTFIT_MODEL_SKY;
+    }
     if (!strcasecmp(fitVarModeString, "MODEL") || !strcasecmp(fitVarModeString, "MODEL_VAR")) {
 	return PM_SOURCE_PHOTFIT_MODEL_VAR;
@@ -142,12 +129,6 @@
     return PM_SOURCE_PHOTFIT_NONE;
 }
-# endif
-
-# if (HAVE_MODEL_VAR)
-bool psphotFitSourcesLinearReadout (psMetadata *recipe, pmReadout *readout, psArray *sources, pmPSF *psf, bool final, pmSourceFitVarMode fitVarMode)
-# else
-bool psphotFitSourcesLinearReadout (psMetadata *recipe, pmReadout *readout, psArray *sources, pmPSF *psf, bool final)
-# endif
-{
+
+bool psphotFitSourcesLinearReadout (psMetadata *recipe, pmReadout *readout, psArray *sources, pmPSF *psf, bool final, pmSourceFitVarMode fitVarMode) {
     bool status;
     float x;
@@ -185,11 +166,4 @@
     if (psRegionIsNaN (AnalysisRegion)) psAbort("analysis region mis-defined");
 
-# if (!HAVE_MODEL_VAR)
-    bool CONSTANT_PHOTOMETRIC_WEIGHTS =
-        psMetadataLookupBool(&status, recipe, "CONSTANT_PHOTOMETRIC_WEIGHTS");
-    if (!status) {
-        psAbort("You must provide a value for the BOOL recipe CONSTANT_PHOTOMETRIC_WEIGHTS");
-    }
-# endif
     int SKY_FIT_ORDER = psMetadataLookupS32(&status, recipe, "SKY_FIT_ORDER");
     if (!status) {
@@ -299,5 +273,5 @@
 	    source->mode |= PM_SOURCE_MODE_PSFMODEL;
 	}	    
-	
+
         psArrayAdd (fitSources, 100, source);
     }
@@ -319,5 +293,4 @@
     psSparseBorder *border = psSparseBorderAlloc (sparse, nBorder);
 
-# if (HAVE_MODEL_VAR)
     // if fitVarMode is MODEL_VAR, then we need to generate the model image variance
     // XXX we have two possibilities here: 
@@ -327,5 +300,4 @@
 
     // 2) do a single pass, and use the model guess to define the model variance (but do I trust the Model Guess?)
-# endif
 
     // fill out the sparse matrix elements and border elements (B)
@@ -336,58 +308,29 @@
 
         // diagonal elements of the sparse matrix (auto-cross-product)
-# if (HAVE_MODEL_VAR)
-        f = pmSourceModelDotModel (SRCi, SRCi, fitVarMode, covarFactor, maskVal);
-# else
-        f = pmSourceModelDotModel (SRCi, SRCi, CONSTANT_PHOTOMETRIC_WEIGHTS, covarFactor, maskVal);
-# endif
-        psSparseMatrixElement (sparse, i, i, f);
-
-# if (HAVE_MODEL_VAR)
-        // if we have used CONSTANT errors, then we need to calculate the value of the parameter error
+        float MM = pmSourceModelDotModel (SRCi, SRCi, fitVarMode, covarFactor, maskVal);
+        psSparseMatrixElement (sparse, i, i, MM);
+
+        // if we have used CONSTANT errors, then we need to re-calculate the value of the parameter error
         if (fitVarMode != PM_SOURCE_PHOTFIT_IMAGE_VAR) {
             float var = pmSourceModelDotModel (SRCi, SRCi, PM_SOURCE_PHOTFIT_IMAGE_VAR, covarFactor, maskVal);
             errors->data.F32[i] = 1.0 / sqrt(var);
         } else {
-            errors->data.F32[i] = 1.0 / sqrt(f);
-        }
-# else
-        // the formal error depends on the weighting scheme
-        if (CONSTANT_PHOTOMETRIC_WEIGHTS) {
-            float var = pmSourceModelDotModel (SRCi, SRCi, false, covarFactor, maskVal);
-            errors->data.F32[i] = 1.0 / sqrt(var);
-        } else {
-            errors->data.F32[i] = 1.0 / sqrt(f);
-        }
-# endif
+            errors->data.F32[i] = 1.0 / sqrt(MM);
+        }
 
         // find the image x model value
-# if (HAVE_MODEL_VAR)
-        f = pmSourceDataDotModel (SRCi, SRCi, fitVarMode, covarFactor, maskVal);
-# else
-        f = pmSourceDataDotModel (SRCi, SRCi, CONSTANT_PHOTOMETRIC_WEIGHTS, covarFactor, maskVal);
-# endif
-        psSparseVectorElement (sparse, i, f);
+        float FM = pmSourceDataDotModel (SRCi, SRCi, fitVarMode, covarFactor, maskVal);
+        psSparseVectorElement (sparse, i, FM);
 
         // add the per-source variances (border region)
         switch (SKY_FIT_ORDER) {
           case 1:
-# if (HAVE_MODEL_VAR)
             f = pmSourceModelWeight (SRCi, 1, fitVarMode, covarFactor, maskVal);
             psSparseBorderElementB (border, i, 1, f);
             f = pmSourceModelWeight (SRCi, 2, fitVarMode, covarFactor, maskVal);
             psSparseBorderElementB (border, i, 2, f);
-# else
-            f = pmSourceModelWeight (SRCi, 1, CONSTANT_PHOTOMETRIC_WEIGHTS, covarFactor, maskVal);
-            psSparseBorderElementB (border, i, 1, f);
-            f = pmSourceModelWeight (SRCi, 2, CONSTANT_PHOTOMETRIC_WEIGHTS, covarFactor, maskVal);
-            psSparseBorderElementB (border, i, 2, f);
-# endif
 
           case 0:
-# if (HAVE_MODEL_VAR)
             f = pmSourceModelWeight (SRCi, 0, fitVarMode, covarFactor, maskVal);
-# else
-            f = pmSourceModelWeight (SRCi, 0, CONSTANT_PHOTOMETRIC_WEIGHTS, covarFactor, maskVal);
-# endif
             psSparseBorderElementB (border, i, 0, f);
             break;
@@ -409,9 +352,5 @@
 
             // got an overlap; calculate cross-product and add to output array
-# if (HAVE_MODEL_VAR)
             f = pmSourceModelDotModel (SRCi, SRCj, fitVarMode, covarFactor, maskVal);
-# else
-            f = pmSourceModelDotModel (SRCi, SRCj, CONSTANT_PHOTOMETRIC_WEIGHTS, covarFactor, maskVal);
-# endif
             psSparseMatrixElement (sparse, j, i, f);
         }
@@ -451,9 +390,5 @@
 
     // set the sky, sky_x, sky_y components of border matrix
-# if (HAVE_MODEL_VAR)
     SetBorderMatrixElements (border, readout, fitSources, (fitVarMode == PM_SOURCE_PHOTFIT_CONST), SKY_FIT_ORDER, markVal);
-# else
-    SetBorderMatrixElements (border, readout, fitSources, CONSTANT_PHOTOMETRIC_WEIGHTS, SKY_FIT_ORDER, markVal);
-# endif
 
     psSparseConstraint constraint;
@@ -500,4 +435,9 @@
     psLogMsg ("psphot.ensemble", PS_LOG_MINUTIA, "sub models: %f sec (%d elements)\n", psTimerMark ("psphot.linear"), sparse->Nelem);
 
+    // mean stats on fit windows
+    float sumRadius = 0.0;
+    float sumPixels = 0.0;
+    float sumSource = 0.0;
+
     // measure chisq for each source
     // for (int i = 0; final && (i < fitSources->n); i++) {
@@ -505,4 +445,10 @@
         pmSource *source = fitSources->data[i];
         pmModel *model = pmSourceGetModel (NULL, source);
+
+	// accumulate fit windows statistics
+        sumRadius += model->fitRadius;
+	sumPixels += M_PI*PS_SQR(model->fitRadius);
+	sumSource += 1.0;
+
         if (!(source->mode & PM_SOURCE_MODE_NONLINEAR_FIT)) {
 	    model->nPar = 1; // LINEAR-only sources have 1 parameter; NONLINEAR sources have their original value
@@ -511,4 +457,7 @@
     }
     psLogMsg ("psphot.ensemble", PS_LOG_MINUTIA, "get chisqs: %f sec (%d elements)\n", psTimerMark ("psphot.linear"), sparse->Nelem);
+
+    float meanRadius = sumRadius / sumSource;
+    float meanPixels = sumPixels / sumSource;
 
     // psFree (index);
@@ -520,5 +469,5 @@
     psFree (border);
 
-    psLogMsg ("psphot.ensemble", PS_LOG_WARN, "measure ensemble of PSFs: %f sec\n", psTimerMark ("psphot.linear"));
+    psLogMsg ("psphot.ensemble", PS_LOG_WARN, "measure ensemble of PSFs (mean radius = %f pixels, mean area = %f pixels: %f sec\n", meanRadius, meanPixels, psTimerMark ("psphot.linear"));
 
     psphotVisualPlotChisq (sources);
@@ -531,7 +480,4 @@
     return true;
 }
-
-// XXX do we need this?
-// XXX disallow the simultaneous sky fit and remove this code...
 
 // Calculate the weight terms for the sky fit component of the matrix.  This function operates
@@ -614,5 +560,4 @@
 }
 
-# if (HAVE_MODEL_VAR)
 bool psphotModelBackgroundReadout(psImage *model,  // Model image
 				  psImage *modelStdev, // Model stdev image
@@ -624,5 +569,5 @@
     );
 
-bool psphotGenerateModelVariance (pmConfig *config, const pmFPAview *view, pmFPAfile *file, int index, psMetadata *recipe, pmReadout *readout, psArray *sources) {
+ bool psphotGenerateModelVariance (pmConfig *config, const pmFPAview *view, pmFPAfile *file, int index, psMetadata *recipe, pmReadout *readout, psArray *sources) {
 
     bool status = false;
@@ -650,4 +595,5 @@
     psImage *varModelStdev = psImageAlloc(backBinning->nXruff, backBinning->nYruff, PS_TYPE_F32); // Background model
 
+    // generate an image of the mean variance image in DN
     if (!psphotModelBackgroundReadout(varModel, varModelStdev, NULL, readout, backBinning, config, true)) {
 	psError(PS_ERR_UNKNOWN, false, "Unable to generate background model");
@@ -664,10 +610,17 @@
 	return false;
     }
-
     psFree (varModel);
     psFree (varModelStdev);
 
-    // XXX for a test:
-    psphotSaveImage (NULL, modelVar, "model.bck.fits");
+    float gain = 1.0;  // accept 1.0 as a default since it is not critical to the analysis
+    pmCell *cell = readout->parent; // The parent cell
+    if (cell) {
+      gain = psMetadataLookupF32(&status, cell->concepts, "CELL.GAIN"); // Cell gain
+      if (!status) {
+	gain = 1.0;	      // set note above
+      }
+    }
+    if (gain > 2.0) { /* warn? */ }
+    // XXX we are not actually using the gain, but need to test it to avoid gcc pedantic warnings
 
     // insert all of the source models
@@ -690,10 +643,9 @@
 
 	// add the source model to the model variance image
+	// XXX note that this should be added with gain applied
+	// var_DN = flux_DN / gain [e/DN]
+	// to do this requires an API upgrade...
 	pmSourceAdd (source, PM_MODEL_OP_MODELVAR, maskVal);
     }
-
-    // XXX for a test:
-    psphotSaveImage (NULL, modelVar, "model.var.fits");
-    psphotSaveImage (NULL, readout->variance, "image.var.fits");
 
     // we save the model variance for future reference
@@ -724,3 +676,2 @@
     return true;
 }
-# endif
Index: trunk/psphot/src/psphotImageLoop.c
===================================================================
--- trunk/psphot/src/psphotImageLoop.c	(revision 33994)
+++ trunk/psphot/src/psphotImageLoop.c	(revision 34086)
@@ -164,4 +164,5 @@
     pmFPAfileActivate (config->files, true, "PSPHOT.EXPNUM");
     while ((chip = pmFPAviewNextChip (view, input->fpa, 1)) != NULL) {
+        if (! chip->process || ! chip->file_exists) { continue; }
         if (!pmFPAfileIOChecks (config, view, PM_FPA_BEFORE)) ESCAPE ("failed attempting to load EXPNUM input for Chip in psphot.");
         while ((cell = pmFPAviewNextCell (view, input->fpa, 1)) != NULL) {
Index: trunk/psphot/src/psphotMakeResiduals.c
===================================================================
--- trunk/psphot/src/psphotMakeResiduals.c	(revision 33994)
+++ trunk/psphot/src/psphotMakeResiduals.c	(revision 34086)
@@ -177,8 +177,6 @@
 
                 mflux = 0;
-                bool offImage = false;
                 if (psImageInterpolate (&flux, &dflux, &mflux, ix, iy, interp) == PS_INTERPOLATE_STATUS_OFF) {
                     // This pixel is off the image
-                    offImage = true;
                     fmasks->data.PS_TYPE_VECTOR_MASK_DATA[i] = badMask;
                 } else {
Index: trunk/psphot/src/psphotMaskReadout.c
===================================================================
--- trunk/psphot/src/psphotMaskReadout.c	(revision 33994)
+++ trunk/psphot/src/psphotMaskReadout.c	(revision 34086)
@@ -94,5 +94,5 @@
 
     // test output of files at this stage
-    if (psTraceGetLevel("psphot") >= 5) {
+    if (psTraceGetLevel("psphot.imsave") >= 5) {
         psphotSaveImage (NULL, readout->image,  "image.fits");
         psphotSaveImage (NULL, readout->mask,   "mask.fits");
@@ -105,2 +105,99 @@
     return true;
 }
+
+// XXX this function and support below was created to test the theory that the faint-end
+// bias results from the Poisson variation of the background pixels.  This is NOT the
+// case.  Using the code below maintains the faint-end bias.
+bool psphotUpdateVariance (pmConfig *config, const pmFPAview *view, const char *filerule) {
+
+    bool status = false;
+
+    // select the appropriate recipe information
+    psMetadata *recipe  = psMetadataLookupPtr (&status, config->recipes, PSPHOT_RECIPE);
+    psAssert (recipe, "missing recipe?");
+
+    int num = psphotFileruleCount(config, filerule);
+
+    // loop over the available readouts
+    for (int i = 0; i < num; i++) {
+
+	// Generate the mask and weight images, including the user-defined analysis region of interest
+	if (!psphotUpdateVarianceReadout (config, view, filerule, i, recipe)) {
+	    psError (PSPHOT_ERR_CONFIG, false, "failed to generate mask for %s entry %d", filerule, i);
+	    return false;
+	}
+    }
+    return true;
+}
+
+// determine the mean variance image (equivalent to the background model, but for the variance image)
+// set the variance image to the MAX(input, mean)
+bool psphotUpdateVarianceReadout (pmConfig *config, const pmFPAview *view, const char *filerule, int index, psMetadata *recipe) {
+
+    pmFPAfile *file = pmFPAfileSelectSingle(config->files, filerule, index); // File of interest
+    psAssert (file, "missing file?");
+
+    // find the currently selected readout
+    pmReadout  *readout = pmFPAviewThisReadout (view, file->fpa);
+    psAssert (readout, "missing readout?");
+
+    pmSourceFitVarMode varMode = psphotGetFitVarMode (recipe);
+    if (varMode == PM_SOURCE_PHOTFIT_NONE) {
+      psError (PSPHOT_ERR_CONFIG, false, "need valid LINEAR_FIT_VARIANCE_MODE");
+      return false;
+    }
+
+    // make this an option via the recipe
+    if (varMode != PM_SOURCE_PHOTFIT_MODEL_SKY) return true;
+
+    // create a model variance image (full-scale image to take result of psImageUnbin below)
+    psImage *modelVar = psImageCopy (NULL, readout->variance, PS_TYPE_F32);
+
+    // find the binning information
+    psImageBinning *backBinning = psphotBackgroundBinning (modelVar, config);
+    assert (backBinning);
+    
+    psImage *varModel = psImageAlloc(backBinning->nXruff, backBinning->nYruff, PS_TYPE_F32); // Background model
+    psImage *varModelStdev = psImageAlloc(backBinning->nXruff, backBinning->nYruff, PS_TYPE_F32); // Background model
+
+    if (!psphotModelBackgroundReadout(varModel, varModelStdev, NULL, readout, backBinning, config, true)) {
+	psError(PS_ERR_UNKNOWN, false, "Unable to generate background model");
+	psFree (varModel);
+	psFree (varModelStdev);
+	return false;
+    }
+
+    // linear interpolation to full-scale
+    if (!psImageUnbin (modelVar, varModel, backBinning)) {
+	psError (PSPHOT_ERR_PROG, true, "inconsistent sizes for unbinning");
+	psFree (varModel);
+	psFree (varModelStdev);
+	return false;
+    }
+
+    // XXX save these?
+    psFree (varModel);
+    psFree (varModelStdev);
+
+    psImage *im = readout->image;
+    psImage *wt = readout->variance;
+    for (int j = 0; j < im->numRows; j++) {
+      for (int i = 0; i < im->numCols; i++) {
+	if (!isfinite(im->data.F32[j][i])) continue;
+	if (!isfinite(wt->data.F32[j][i])) continue;
+	// XXX for a test, make variance constant wt->data.F32[j][i] = PS_MAX(wt->data.F32[j][i], modelVar->data.F32[j][i]);
+	wt->data.F32[j][i] = modelVar->data.F32[j][i];
+      }
+    }
+
+    // test output of files at this stage
+    if (psTraceGetLevel("psphot.imsave") >= 5) {
+        psphotSaveImage (NULL, readout->image,  "image.varsky.fits");
+        psphotSaveImage (NULL, readout->mask,   "mask.varsky.fits");
+        psphotSaveImage (NULL, readout->variance, "variance.varsky.fits");
+    }
+
+    psFree (modelVar);
+
+    return true;
+}
Index: trunk/psphot/src/psphotPetrosianStats.c
===================================================================
--- trunk/psphot/src/psphotPetrosianStats.c	(revision 33994)
+++ trunk/psphot/src/psphotPetrosianStats.c	(revision 34086)
@@ -46,5 +46,5 @@
 
     bool anyPetro = false;
-    bool manyPetro = false;
+    // bool manyPetro = false;  XXX not used
     bool above = true;
     float Asum = 0.0;
@@ -122,5 +122,5 @@
 	    }
 	    above = false;
-	    if (anyPetro) manyPetro = true;
+	    // if (anyPetro) manyPetro = true;
 	    anyPetro = true;
 	}
@@ -212,5 +212,5 @@
     source->extpars->petrosianR50    = R50;
     source->extpars->petrosianR90    = R90;
-    source->extpars->petrosianFill      = petApix / petArea;
+    source->extpars->petrosianFill   = petApix / petArea;
     
     // XXX add the errors
Index: trunk/psphot/src/psphotReadout.c
===================================================================
--- trunk/psphot/src/psphotReadout.c	(revision 33994)
+++ trunk/psphot/src/psphotReadout.c	(revision 34086)
@@ -9,4 +9,49 @@
 }
 
+# if (0)
+// TEST CODE, can be removed
+bool psphotDumpFlux (pmConfig *config, const pmFPAview *view, const char *filerule) {
+
+    bool status = false;
+
+    // find the currently selected readout
+    pmFPAfile *file = pmFPAfileSelectSingle(config->files, filerule, 0); // 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?");
+
+    static int npass;
+    char filename[64];
+    snprintf (filename, 64, "mags.%d.dat", npass);
+    FILE *ftest = fopen (filename, "w");
+    for (int j = 0; j < sources->n; j++) {
+	pmSource *source = sources->data[j];
+
+	float psfMag;
+	status = pmSourcePhotometryModel (&psfMag, NULL, source->modelPSF);
+
+	float psfMagNorm;
+	float Io = source->modelPSF->params->data.F32[PM_PAR_I0];
+	source->modelPSF->params->data.F32[PM_PAR_I0] = 1.0;
+	status = pmSourcePhotometryModel (&psfMagNorm, NULL, source->modelPSF);
+	source->modelPSF->params->data.F32[PM_PAR_I0] = Io;
+
+	// double apTrend = pmTrend2DEval (psf->ApTrend, (float)source->peak->x, (float)source->peak->y);
+	fprintf (ftest, "%d %d %d  %f %f %f %f  %f %f\n", j, source->peak->x, source->peak->y, source->modelPSF->params->data.F32[PM_PAR_I0], source->modelPSF->params->data.F32[PM_PAR_SXX], source->modelPSF->params->data.F32[PM_PAR_SYY], source->modelPSF->params->data.F32[PM_PAR_SXY], psfMag, psfMagNorm);
+    }
+    fclose (ftest);
+    npass++;
+
+    return true;
+}
+# endif
+
 bool psphotReadout(pmConfig *config, const pmFPAview *view, const char *filerule) {
 
Index: trunk/psphot/src/psphotReadoutCleanup.c
===================================================================
--- trunk/psphot/src/psphotReadoutCleanup.c	(revision 33994)
+++ trunk/psphot/src/psphotReadoutCleanup.c	(revision 34086)
@@ -7,5 +7,5 @@
 
     // remove internal pmFPAfiles, if created
-    if (psErrorCodeLast() == PSPHOT_ERR_DATA) {
+    if (psErrorCodeLast() == (psErrorCode) PSPHOT_ERR_DATA) {
         psErrorStackPrint(stderr, "Error in the psphot readout analysis");
         psErrorClear();
Index: trunk/psphot/src/psphotSignificanceImage.c
===================================================================
--- trunk/psphot/src/psphotSignificanceImage.c	(revision 33994)
+++ trunk/psphot/src/psphotSignificanceImage.c	(revision 34086)
@@ -8,5 +8,4 @@
     float SIGMA_SMTH, NSIGMA_SMTH;
     bool status = false;
-    bool guess = false;
 
     // smooth the image and variance map
@@ -35,5 +34,4 @@
         SIGMA_SMTH  = 0.5*(fwhmMajor + fwhmMinor) / (2.0*sqrt(2.0*log(2.0)));
         NSIGMA_SMTH = psMetadataLookupF32 (&status, recipe, "PEAKS_SMOOTH_NSIGMA");
-        guess = false;
     } else {
         // if we do not know the FWHM, use the guess smoothing kernel supplied.
@@ -43,5 +41,4 @@
         NSIGMA_SMTH = psMetadataLookupF32 (&status, recipe, "PEAKS_SMOOTH_NSIGMA");
         PS_ASSERT (status, NULL);
-        guess = true;
     }
     // record the actual smoothing sigma
Index: trunk/psphot/src/psphotSourceFits.c
===================================================================
--- trunk/psphot/src/psphotSourceFits.c	(revision 33994)
+++ trunk/psphot/src/psphotSourceFits.c	(revision 34086)
@@ -224,14 +224,4 @@
     if (source->type == PM_SOURCE_TYPE_DEFECT) return false;
     if (source->type == PM_SOURCE_TYPE_SATURATED) return false;
-
-# define TEST_X -420.0
-# define TEST_Y 300.0
-   
-    if ((fabs(source->peak->xf - TEST_X) < 5) && (fabs(source->peak->yf - TEST_Y) < 5)) {
-	fprintf (stderr, "test galaxy\n");
-    }
-
-# undef TEST_X
-# undef TEST_Y
 
     // set the radius based on the footprint (also sets the mask pixels)
@@ -513,11 +503,4 @@
     }
 
-# define TEST_X -540.0
-# define TEST_Y 540.0
-   
-    if ((fabs(source->peak->xf - TEST_X) < 5) && (fabs(source->peak->yf - TEST_Y) < 5)) {
-	psTraceSetLevel("psModules.objects.pmPCM_MinimizeChisq", 5);
-    }
-
     float t1, t2, t4, t5;
     if (TIMING) { psTimerStart ("psphotFitPCM"); }
@@ -569,8 +552,4 @@
     }
 
-    if ((fabs(source->peak->xf - TEST_X) < 5) && (fabs(source->peak->yf - TEST_Y) < 5)) {
-	psTraceSetLevel("psModules.objects.pmPCM_MinimizeChisq", 0);
-    }
-
     // psTraceSetLevel("psLib.math.psMinimizeLMChi2", 0);
     psFree (pcm);
@@ -578,7 +557,4 @@
     return model;
 }
-
-# undef TEST_X
-# undef TEST_Y
 
 // note that these should be 1/2n of the standard sersic index
@@ -603,11 +579,4 @@
     float xMin = NAN;
     float chiSquare[N_INDEX_GUESS];
-
-# define TEST_X -540.0
-# define TEST_Y 540.0
-   
-    if ((fabs(source->peak->xf - TEST_X) < 5) && (fabs(source->peak->yf - TEST_Y) < 5)) {
-	psTraceSetLevel("psModules.objects.pmPCM_MinimizeChisq", 5);
-    }
 
     for (int i = 0; i < N_INDEX_GUESS; i++) {
@@ -635,8 +604,4 @@
     assert (iMin >= 0);
 
-    if ((fabs(source->peak->xf - TEST_X) < 5) && (fabs(source->peak->yf - TEST_Y) < 5)) {
-	psTraceSetLevel("psModules.objects.pmPCM_MinimizeChisq", 0);
-    }
-
     model->flags = PM_MODEL_STATUS_NONE; // do not attempt to handle failures here, let the next iteration deal with it
     model->params->data.F32[PM_PAR_7] = 0.5/indexGuess[iMin];
@@ -666,8 +631,4 @@
     float xMin = NAN;
     float chiSquare[N_INDEX_GUESS];
-
-    if ((fabs(source->peak->xf - TEST_X) < 5) && (fabs(source->peak->yf - TEST_Y) < 5)) {
-	psTraceSetLevel("psModules.objects.pmPCM_MinimizeChisq", 5);
-    }
 
     for (int i = 0; i < N_INDEX_GUESS; i++) {
@@ -701,8 +662,4 @@
     }
     assert (iMin >= 0);
-    
-    if ((fabs(source->peak->xf - TEST_X) < 5) && (fabs(source->peak->yf - TEST_Y) < 5)) {
-	psTraceSetLevel("psModules.objects.pmPCM_MinimizeChisq", 0);
-    }
 
     model->flags = PM_MODEL_STATUS_NONE; // do not attempt to handle failures here, let the next iteration deal with it