Index: branches/simmosaic_branches/psphot/src/psphotMagnitudes.c
===================================================================
--- branches/simmosaic_branches/psphot/src/psphotMagnitudes.c	(revision 24860)
+++ branches/simmosaic_branches/psphot/src/psphotMagnitudes.c	(revision 27839)
@@ -1,14 +1,58 @@
 # include "psphotInternal.h"
 
-bool psphotMagnitudes(pmConfig *config, pmReadout *readout, const pmFPAview *view, psArray *sources, pmPSF *psf) {
+bool psphotMagnitudes (pmConfig *config, const pmFPAview *view)
+{
+    bool status = true;
+
+    // select the appropriate recipe information
+    psMetadata *recipe  = psMetadataLookupPtr (&status, config->recipes, PSPHOT_RECIPE);
+    psAssert (recipe, "missing recipe?");
+
+    int num = psMetadataLookupS32 (&status, config->arguments, "PSPHOT.INPUT.NUM");
+    psAssert (status, "programming error: must define PSPHOT.INPUT.NUM");
+
+    // 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, "PSPHOT.INPUT", 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 (!psphotMagnitudesReadout (config, recipe, view, readout, sources, psf)) {
+            psError (PSPHOT_ERR_CONFIG, false, "failed to measure magnitudes for PSPHOT.INPUT entry %d", i);
+	    return false;
+	}
+    }
+    return true;
+}
+
+bool psphotMagnitudesReadout(pmConfig *config, psMetadata *recipe, const pmFPAview *view, pmReadout *readout, psArray *sources, pmPSF *psf) {
 
     bool status = false;
     int Nap = 0;
 
+    if (!sources->n) {
+	psLogMsg ("psphot", PS_LOG_INFO, "no sources, skipping source magnitudes");
+	return true;
+    }
+	
     psTimerStart ("psphot.mags");
-
-    // select the appropriate recipe information
-    psMetadata *recipe  = psMetadataLookupPtr (&status, config->recipes, PSPHOT_RECIPE);
-    assert (recipe);
 
     // determine the number of allowed threads
@@ -41,10 +85,12 @@
     }
 
-    bool IGNORE_GROWTH = psMetadataLookupBool (&status, recipe, "IGNORE_GROWTH");
+    bool IGNORE_GROWTH  = psMetadataLookupBool (&status, recipe, "IGNORE_GROWTH");
     bool INTERPOLATE_AP = psMetadataLookupBool (&status, recipe, "INTERPOLATE_AP");
+    bool DIFF_STATS     = psMetadataLookupBool (&status, recipe, "INTERPOLATE_AP");
 
     pmSourcePhotometryMode photMode = PM_SOURCE_PHOT_APCORR | PM_SOURCE_PHOT_WEIGHT;
     if (!IGNORE_GROWTH) photMode |= PM_SOURCE_PHOT_GROWTH;
     if (INTERPOLATE_AP) photMode |= PM_SOURCE_PHOT_INTERP;
+    if (DIFF_STATS)     photMode |= PM_SOURCE_PHOT_DIFFSTATS;
 
     // choose Cx, Cy (see psphotThreadTools.c for overview of the concepts)
@@ -71,4 +117,5 @@
             PS_ARRAY_ADD_SCALAR(job->args, photMode, PS_TYPE_S32);
             PS_ARRAY_ADD_SCALAR(job->args, maskVal,  PS_TYPE_IMAGE_MASK);
+            PS_ARRAY_ADD_SCALAR(job->args, markVal,  PS_TYPE_IMAGE_MASK);
             PS_ARRAY_ADD_SCALAR(job->args, 0,        PS_TYPE_S32); // this is used as a return value for nAp
 
@@ -102,5 +149,5 @@
                 fprintf (stderr, "error with job\n");
             } else {
-                psScalar *scalar = job->args->data[7];
+                psScalar *scalar = job->args->data[8];
                 Nap += scalar->data.S32;
             }
@@ -127,9 +174,26 @@
     pmSourcePhotometryMode photMode = PS_SCALAR_VALUE(job->args->data[5],S32);
     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);
 
     for (int i = 0; i < sources->n; i++) {
         pmSource *source = (pmSource *) sources->data[i];
-        status = pmSourceMagnitudes (source, psf, photMode, maskVal);
+
+        // replace object in image
+        if (source->tmpFlags & PM_SOURCE_TMPF_SUBTRACTED) {
+            pmSourceAdd (source, PM_MODEL_OP_FULL, maskVal);
+        }
+
+	// 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);
+
+        status = pmSourceMagnitudes (source, psf, photMode, maskVal); // maskVal includes markVal
         if (status && isfinite(source->apMag)) Nap ++;
+
+	// clear the mask bit 
+	psImageMaskPixels (source->maskObj, "AND", PS_NOT_IMAGE_MASK(markVal));
+
+        // re-subtract the object, leave local sky
+        pmSourceSub (source, PM_MODEL_OP_FULL, maskVal);
 
         if (backModel) {
@@ -155,48 +219,9 @@
 
     // change the value of a scalar on the array (wrap this and put it in psArray.h)
-    psScalar *scalar = job->args->data[7];
+    psScalar *scalar = job->args->data[8];
     scalar->data.S32 = Nap;
 
     return true;
 }
-
-# if (0)
-bool psphotMagnitudes_Unthreaded (int *nap, psArray *sources, pmPSF *psf, psImageBinning *binning, pmReadout *backModel, pmReadout *backStdev, pmSourcePhotometryMode photMode, psImageMaskType maskVal) {
-
-    bool status;
-    int Nap = 0;
-
-    for (int i = 0; i < sources->n; i++) {
-        pmSource *source = (pmSource *) sources->data[i];
-        status = pmSourceMagnitudes (source, psf, photMode, maskVal);
-        if (status && isfinite(source->apMag)) Nap ++;
-
-        if (backModel) {
-            psAssert (binning, "if backModel is defined, so should binning be");
-            source->sky = psImageUnbinPixel(source->peak->x, source->peak->y, backModel->image, binning);
-            if (isnan(source->sky) && false) {
-                psLogMsg ("psphot.magnitudes", PS_LOG_DETAIL, "error setting pmSource.sky");
-            }
-        } else {
-            source->sky = NAN;
-        }
-
-        if (backStdev) {
-            psAssert (binning, "if backStdev is defined, so should binning be");
-            source->skyErr = psImageUnbinPixel(source->peak->x, source->peak->y, backStdev->image, binning);
-            if (isnan(source->skyErr) && false) {
-                psLogMsg ("psphot.magnitudes", PS_LOG_DETAIL, "error setting pmSource.skyErr");
-            }
-        } else {
-            source->skyErr = NAN;
-        }
-    }
-
-    // change the value of a scalar on the array (wrap this and put it in psArray.h)
-    *nap = Nap;
-
-    return true;
-}
-# endif
 
 bool psphotPSFWeights(pmConfig *config, pmReadout *readout, const pmFPAview *view, psArray *sources) {