Index: trunk/psModules/src/objects/pmSourcePhotometry.c
===================================================================
--- trunk/psModules/src/objects/pmSourcePhotometry.c	(revision 25980)
+++ trunk/psModules/src/objects/pmSourcePhotometry.c	(revision 27531)
@@ -109,7 +109,10 @@
 	psAssert (isfinite(fluxScale), "how can the flux scale be invalid? source at %d, %d\n", source->peak->x, source->peak->y);
 	psAssert (fluxScale > 0.0, "how can the flux scale be negative? source at %d, %d\n", source->peak->x, source->peak->y);
-	source->psfMag = -2.5*log10(fluxScale * source->modelPSF->params->data.F32[PM_PAR_I0]);
+	source->psfFlux = fluxScale * source->modelPSF->params->data.F32[PM_PAR_I0];
+	source->psfFluxErr = fluxScale * source->modelPSF->dparams->data.F32[PM_PAR_I0];
+	source->psfMag = -2.5*log10(source->psfFlux);
     } else {
-        status = pmSourcePhotometryModel (&source->psfMag, source->modelPSF);
+        status = pmSourcePhotometryModel (&source->psfMag, &source->psfFlux, source->modelPSF);
+	source->psfFluxErr = source->psfFlux * (source->modelPSF->dparams->data.F32[PM_PAR_I0] / source->modelPSF->params->data.F32[PM_PAR_I0]);
     }
 
@@ -119,5 +122,5 @@
 	for (int i = 0; i < source->modelFits->n; i++) {
 	    pmModel *model = source->modelFits->data[i];
-	    status = pmSourcePhotometryModel (&model->mag, model);
+	    status = pmSourcePhotometryModel (&model->mag, NULL, model);
 	    if (model == source->modelEXT) foundEXT = true;
 	}
@@ -125,9 +128,9 @@
 	    source->extMag = source->modelEXT->mag;
 	} else {
-	    status = pmSourcePhotometryModel (&source->extMag, source->modelEXT);
+	    status = pmSourcePhotometryModel (&source->extMag, NULL, source->modelEXT);
 	}
     } else {
 	if (source->modelEXT) {
-	    status = pmSourcePhotometryModel (&source->extMag, source->modelEXT);
+	    status = pmSourcePhotometryModel (&source->extMag, NULL, source->modelEXT);
 	}
     }
@@ -143,4 +146,9 @@
     if (mode & PM_SOURCE_PHOT_WEIGHT) {
         pmSourcePixelWeight (&source->pixWeight, model, source->maskObj, maskVal);
+    }
+
+    // measure the contribution of included pixels
+    if (mode & PM_SOURCE_PHOT_DIFFSTATS) {
+        pmSourceMeasureDiffStats (source, maskVal);
     }
 
@@ -217,21 +225,26 @@
 
 // return source model magnitude
-bool pmSourcePhotometryModel (float *fitMag, pmModel *model)
-{
-    PS_ASSERT_PTR_NON_NULL(fitMag, false);
-    if (model == NULL) {
-        return false;
-    }
-
-    float fitSum = 0;
-    *fitMag = NAN;
+bool pmSourcePhotometryModel (float *fitMag, float *fitFlux, pmModel *model)
+{
+    psAssert (fitMag || fitFlux, "at least one of magnitude or flux must be requested (not NULL)");
+    if (model == NULL) return false;
+
+    float mag  = NAN;
+    float flux = NAN;
 
     // measure fitMag
-    fitSum = model->modelFlux (model->params);
-    if (fitSum <= 0)
-        return false;
-    if (!isfinite(fitSum))
-        return false;
-    *fitMag = -2.5*log10(fitSum);
+    flux = model->modelFlux (model->params);
+    if (flux > 0) {
+	mag = -2.5*log10(flux);
+    }
+    if (fitMag) {
+	*fitMag = mag;
+    }
+    if (fitFlux) {
+	*fitFlux = flux;
+    }
+
+    if (flux <= 0) return false;
+    if (!isfinite(flux)) return false;
 
     return (true);
@@ -356,4 +369,55 @@
 
     *pixWeight = validSum / modelSum;
+    return (true);
+}
+
+# define FLUX_LIMIT 3.0
+
+// return source aperture magnitude
+bool pmSourceMeasureDiffStats (pmSource *source, psImageMaskType maskVal)
+{
+    PS_ASSERT_PTR_NON_NULL(source, false);
+
+    if (source->diffStats == NULL) {
+	source->diffStats = pmSourceDiffStatsAlloc();
+    }
+
+    float fGood = 0.0;
+    float fBad  = 0.0;
+    int   nGood = 0;
+    int   nMask = 0;
+    int   nBad  = 0;
+    
+    psImage *flux     = source->pixels;
+    psImage *variance = source->variance;
+    psImage *mask     = source->maskObj;
+
+    for (int iy = 0; iy < flux->numRows; iy++) {
+	for (int ix = 0; ix < flux->numCols; ix++) {
+            if (mask->data.PS_TYPE_IMAGE_MASK_DATA[iy][ix] & maskVal) {
+		nMask ++;
+                continue;
+	    }
+
+	    float SN = flux->data.F32[iy][ix] / sqrt(variance->data.F32[iy][ix]);
+
+	    if (SN > +FLUX_LIMIT) { 
+		nGood ++;
+		fGood += fabs(flux->data.F32[iy][ix]);
+	    }
+
+	    if (SN < -FLUX_LIMIT) { 
+		nBad ++;
+		fBad += fabs(flux->data.F32[iy][ix]);
+	    }
+	}
+    }
+
+    source->diffStats->nGood      = nGood;
+    source->diffStats->fRatio     = (fGood + fBad         == 0.0) ? NAN : fGood / (fGood + fBad);	   
+    source->diffStats->nRatioBad  = (nGood + nBad         == 0)   ? NAN : nGood / (float) (nGood + nBad);	   
+    source->diffStats->nRatioMask = (nGood + nMask        == 0)   ? NAN : nGood / (float) (nGood + nMask);	   
+    source->diffStats->nRatioAll  = (nGood + nMask + nBad == 0)   ? NAN : nGood / (float) (nGood + nMask + nBad);
+
     return (true);
 }