Index: /branches/eam_branches/ipp-20111122/psphot/src/psphot.h
===================================================================
--- /branches/eam_branches/ipp-20111122/psphot/src/psphot.h	(revision 33085)
+++ /branches/eam_branches/ipp-20111122/psphot/src/psphot.h	(revision 33086)
@@ -99,5 +99,5 @@
 
 bool            psphotFitSourcesLinear (pmConfig *config, const pmFPAview *view, const char *filerule, bool final);
-bool            psphotFitSourcesLinearReadout (psMetadata *recipe, pmReadout *readout, psArray *sources, pmPSF *psf, bool final);
+bool            psphotFitSourcesLinearReadout (psMetadata *recipe, pmReadout *readout, psArray *sources, pmPSF *psf, bool final, pmSourceFitVarMode fitVarMode);
 
 bool            psphotSourceSize (pmConfig *config, const pmFPAview *view, const char *filerule, bool getPSFsize);
Index: /branches/eam_branches/ipp-20111122/psphot/src/psphotEfficiency.c
===================================================================
--- /branches/eam_branches/ipp-20111122/psphot/src/psphotEfficiency.c	(revision 33085)
+++ /branches/eam_branches/ipp-20111122/psphot/src/psphotEfficiency.c	(revision 33086)
@@ -420,5 +420,5 @@
 
     // psphotFitSourcesLinearReadout subtracts the model fits
-    if (!psphotFitSourcesLinearReadout(recipe, readout, fakeSourcesAll, psf, true)) {
+    if (!psphotFitSourcesLinearReadout(recipe, readout, fakeSourcesAll, psf, true, PM_SOURCE_PHOTFIT_CONST)) {
         psError(PS_ERR_UNKNOWN, false, "Unable to perform linear fit on fake sources.");
         psFree(fakeSources);
Index: /branches/eam_branches/ipp-20111122/psphot/src/psphotFitSourcesLinear.c
===================================================================
--- /branches/eam_branches/ipp-20111122/psphot/src/psphotFitSourcesLinear.c	(revision 33085)
+++ /branches/eam_branches/ipp-20111122/psphot/src/psphotFitSourcesLinear.c	(revision 33086)
@@ -12,4 +12,8 @@
 static bool SetBorderMatrixElements (psSparseBorder *border, pmReadout *readout, psArray *sources, bool constant_weights, int SKY_FIT_ORDER, psImageMaskType markVal);
 
+bool psphotGenerateModelVariance (pmConfig *config, const pmFPAview *view, int index, psMetadata *recipe, pmReadout *readout, psArray *sources);
+pmSourceFitVarMode psphotGetFitVarMode (psMetadata *recipe);
+bool psphotFreeModelVariance (pmReadout *readout, psArray *sources);
+
 // 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)
@@ -29,4 +33,7 @@
 	return false;
     }
+    // MODEL_VAR requires 2 passes -- in the first, we get the rough fluxes; in the second, we
+    // use the flux to define the model variance before fitting the objects.  Other modes only
+    // do a single pass.
     pmSourceFitVarMode fitVarModePass1 = (fitVarMode == PM_SOURCE_PHOTFIT_MODEL_VAR) ? PM_SOURCE_PHOTFIT_CONST : fitVarMode;
 
@@ -66,9 +73,12 @@
 	// re-run the fit against that variance.
 	if (fitVarMode == PM_SOURCE_PHOTFIT_MODEL_VAR) {
-	    // generate the model variance image
-	    if (!psphotGenerateModelVariance (recipe, readout, sources, psf, final, fitVarMode)) {
+	    // generate the model variance image & source pointers
+	    if (!psphotGenerateModelVariance (config, view, i, recipe, readout, sources)) {
 		psError (PSPHOT_ERR_CONFIG, false, "failed to fit sources (linear) for %s entry %d", filerule, i);
 		return false;
 	    }
+
+	    // replace all sources (use TMPF_SUBTRACTED as test)
+	    psphotReplaceAllSourcesReadout (config, view, filerule, i, recipe, false);
 
 	    // rerun fit with correct fitVarMode
@@ -77,4 +87,10 @@
 		return false;
 	    }
+
+	    // free the model variance image & source pointers
+	    if (!psphotFreeModelVariance (readout, sources)) {
+		psError (PSPHOT_ERR_CONFIG, false, "failed to fit sources (linear) for %s entry %d", filerule, i);
+		return false;
+	    }
 	}
 
@@ -86,4 +102,7 @@
 }
 
+// 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
+// set the value to either CONST or IMAGE_VAR
 pmSourceFitVarMode psphotGetFitVarMode (psMetadata *recipe) {
 
@@ -96,5 +115,5 @@
 	    psAbort("You must provide a value for LINEAR_FIT_VARIANCE_MODE or CONSTANT_PHOTOMETRIC_WEIGHTS");
 	}
-	pmsourceFitVarMode fitVarMode = CONSTANT_PHOTOMETRIC_WEIGHTS ? PM_SOURCE_PHOTFIT_CONST : PM_SOURCE_PHOTFIT_IMAGE_VAR;
+	pmSourceFitVarMode fitVarMode = CONSTANT_PHOTOMETRIC_WEIGHTS ? PM_SOURCE_PHOTFIT_CONST : PM_SOURCE_PHOTFIT_IMAGE_VAR;
 	return fitVarMode;
     } 
@@ -109,5 +128,5 @@
     }
     psError (PSPHOT_ERR_CONFIG, false, "Invalid value for LINEAR_FIT_VARIANCE_MODE (%s)", fitVarModeString);
-    return PM_SOURCE_PHOTFIT_MODEL_NONE;
+    return PM_SOURCE_PHOTFIT_NONE;
 }
 
@@ -182,9 +201,4 @@
 	// psAssert (source->peak, "source without peak??");
 	// psAssert (source->peak->footprint, "peak without footprint??");
-
-	// XXX TEST
-        if (source->mode2 & PM_SOURCE_MODE2_MATCHED) {
-	    fprintf (stderr, "forced photometry on supplied source\n");
-	}
 
         // turn this bit off and turn it on again if we pass this test
@@ -381,5 +395,5 @@
 
     // set the sky, sky_x, sky_y components of border matrix
-    SetBorderMatrixElements (border, readout, fitSources, CONSTANT_PHOTOMETRIC_WEIGHTS, SKY_FIT_ORDER, markVal);
+    SetBorderMatrixElements (border, readout, fitSources, (fitVarMode == PM_SOURCE_PHOTFIT_CONST), SKY_FIT_ORDER, markVal);
 
     psSparseConstraint constraint;
@@ -540,11 +554,25 @@
 }
 
-bool psphotGenerateModelVariance (pmConfig *config, const pmFPAview *view, const char *filerule, int index, psMetadata *recipe) {
+bool psphotGenerateModelVariance (pmConfig *config, const pmFPAview *view, int index, psMetadata *recipe, pmReadout *readout, psArray *sources) {
+
+    bool status = false;
+    psRegion fullRegion = psRegionSet (0, 0, 0, 0);
+
+    // 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;
 
     // create a model variance image
     psImage *modelVar = psImageCopy (NULL, readout->variance, PS_TYPE_F32);
 
-    // find the pmFPAfile for the readout (background model is saved on PSPHOT.BACKMDL regardless of 'filename')
-    pmFPAfile *backModellFile = pmFPAfileSelectSingle(config->files, "PSPHOT.BACKMDL", index); // File of interest
+    // find the pmFPAfile for the readout (background model is saved on PSPHOT.BACKMDL regardless of 'filerule')
+    pmFPAfile *backModelFile = pmFPAfileSelectSingle(config->files, "PSPHOT.BACKMDL", index); // File of interest
     assert (backModelFile);
 
@@ -554,5 +582,5 @@
 
     // find the binning information
-    psImageBinning *backBinning = psMetadataLookupPtr(&status, model->analysis, "PSPHOT.BACKGROUND.BINNING");
+    psImageBinning *backBinning = psMetadataLookupPtr(&status, backModel->analysis, "PSPHOT.BACKGROUND.BINNING");
     assert (backBinning);
 
@@ -572,23 +600,17 @@
 	pmSource *source = sources->data[i];
 
-	// XXX which sources need to be injected?
-	// skip non-astronomical objects (very likely defects)
-	if (source->type == PM_SOURCE_TYPE_DEFECT) continue;
-	if (source->type == PM_SOURCE_TYPE_SATURATED) continue;
-      
-	// do not include CRs in the full ensemble fit
-	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 need to identify the region appropriate for the source
+	// skip sources which were not fitted already
+	if (!(source->mode & PM_SOURCE_MODE_LINEAR_FIT)) continue;
+
+	// pixel region appropriate for the source
+	psRegion region = psRegionForImage (source->pixels, fullRegion);
 
 	// define the source->modelVar pixels (view on modelVar image)
 	psAssert (!source->modelVar, "programming error : modelVar should be NULL here");
+	psAssert (source->modelFlux, "programming error : modelFlux should not be NULL here");
 	source->modelVar = psImageSubset(readout->image, region);
 
 	// add the source model to the model variance image
-	pmSourceAdd (source, PM_MODEL_OP_XXX, maskVal);
+	pmSourceAdd (source, PM_MODEL_OP_MODELVAR, maskVal);
     }
 
@@ -599,12 +621,19 @@
     // XXX for a test:
     psphotSaveImage (NULL, modelVar, "model.var.fits");
+    psphotSaveImage (NULL, readout->image, "image.var.fits");
+
+    // we save the model variance for future reference
+    psMetadataAddImage(readout->analysis, PS_LIST_TAIL, "PSPHOT.MODEL.VAR", PS_META_REPLACE, "model variance", modelVar);
+    psFree (modelVar);
 
     return true;
 }
 
-bool psphotFreeModelVariance (pmConfig *config, const pmFPAview *view, const char *filerule, int index, psMetadata *recipe) {
+bool psphotFreeModelVariance (pmReadout *readout, psArray *sources) {
+
+    bool status = false;
 
     // find the binning information
-    psImage *modelVar = psMetadataLookupPtr(&status, model->analysis, "PSPHOT.MODEL.VAR");
+    psImage *modelVar = psMetadataLookupPtr(&status, readout->analysis, "PSPHOT.MODEL.VAR");
     assert (modelVar);
 
Index: /branches/eam_branches/ipp-20111122/psphot/src/psphotModelBackground.c
===================================================================
--- /branches/eam_branches/ipp-20111122/psphot/src/psphotModelBackground.c	(revision 33085)
+++ /branches/eam_branches/ipp-20111122/psphot/src/psphotModelBackground.c	(revision 33086)
@@ -159,6 +159,6 @@
 
     // measure clipped median for subimages
-    psRegion ruffRegion = {0,0,0,0};
-    psRegion fineRegion = {0,0,0,0};
+    psRegion ruffRegion = psRegionSet (0,0,0,0);
+    psRegion fineRegion = psRegionSet (0,0,0,0);
     for (int iy = 0; iy < model->numRows; iy++) {
         for (int ix = 0; ix < model->numCols; ix++) {
Index: /branches/eam_branches/ipp-20111122/psphot/src/psphotTest.c
===================================================================
--- /branches/eam_branches/ipp-20111122/psphot/src/psphotTest.c	(revision 33085)
+++ /branches/eam_branches/ipp-20111122/psphot/src/psphotTest.c	(revision 33086)
@@ -103,5 +103,5 @@
 # if (0)
 
-    psRegion region = {0,0,0,0};        // a region representing the entire array
+psRegion region = psRegionSet (0,0,0,0);        // a region representing the entire array
     psphotTestArguments (&argc, argv);
 
Index: /branches/eam_branches/ipp-20111122/psphot/src/psphotTestSourceOutput.c
===================================================================
--- /branches/eam_branches/ipp-20111122/psphot/src/psphotTestSourceOutput.c	(revision 33085)
+++ /branches/eam_branches/ipp-20111122/psphot/src/psphotTestSourceOutput.c	(revision 33086)
@@ -136,5 +136,5 @@
 
 	    // generate working image for this source
-	    psRegion region = {ix - dx, ix + dx, iy - dy, iy + dy};
+	    psRegion region = psRegionSet(ix - dx, ix + dx, iy - dy, iy + dy);
 
 	    psImage *vM = psImageSubset (imMo, region);