Index: trunk/psModules/src/objects/pmSourcePhotometry.c
===================================================================
--- trunk/psModules/src/objects/pmSourcePhotometry.c	(revision 21511)
+++ trunk/psModules/src/objects/pmSourcePhotometry.c	(revision 25754)
@@ -84,4 +84,5 @@
 
     // we must have a valid model
+    // XXX allow aperture magnitudes for sources without a model
     model = pmSourceGetModel (&isPSF, source);
     if (model == NULL) {
@@ -90,4 +91,5 @@
     }
 
+    // XXX handle negative flux, low-significance
     if (model->dparams->data.F32[PM_PAR_I0] > 0) {
         SN = model->params->data.F32[PM_PAR_I0] / model->dparams->data.F32[PM_PAR_I0];
@@ -101,37 +103,32 @@
 
     // measure PSF model photometry
-    if (psf->FluxScale) {
+    // XXX TEST: do not use flux scale
+    if (0 && psf->FluxScale) {
         // the source peak pixel is guaranteed to be on the image, and only minimally different from the source center
         double fluxScale = pmTrend2DEval (psf->FluxScale, (float)source->peak->x, (float)source->peak->y);
-        if (!isfinite(fluxScale)) {
-            // XXX objects on the edge can be slightly outside -- if we get an
-            // error, use the full fit.
-            psErrorClear();
-            status = pmSourcePhotometryModel (&source->psfMag, source->modelPSF);
-        } else {
-            if (isfinite(fluxScale) && (fluxScale > 0.0)) {
-                source->psfMag = -2.5*log10(fluxScale * source->modelPSF->params->data.F32[PM_PAR_I0]);
-            } else {
-                source->psfMag = NAN;
-            }
-        }
+	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]);
     } else {
         status = pmSourcePhotometryModel (&source->psfMag, source->modelPSF);
     }
 
-    // if we have a collection of model fits, one of them is a pointer to modelEXT?
-    // XXX not guaranteed
+    // if we have a collection of model fits, check if one of them is a pointer to modelEXT
     if (source->modelFits) {
-      for (int i = 0; i < source->modelFits->n; i++) {
-        pmModel *model = source->modelFits->data[i];
-        status = pmSourcePhotometryModel (&model->mag, model);
-      }
-      if (source->modelEXT) {
-        source->extMag = source->modelEXT->mag;
-      }
+        bool foundEXT = false;
+	for (int i = 0; i < source->modelFits->n; i++) {
+	    pmModel *model = source->modelFits->data[i];
+	    status = pmSourcePhotometryModel (&model->mag, model);
+	    if (model == source->modelEXT) foundEXT = true;
+	}
+	if (foundEXT) {
+	    source->extMag = source->modelEXT->mag;
+	} else {
+	    status = pmSourcePhotometryModel (&source->extMag, source->modelEXT);
+	}
     } else {
-      if (source->modelEXT) {
-        status = pmSourcePhotometryModel (&source->extMag, source->modelEXT);
-      }
+	if (source->modelEXT) {
+	    status = pmSourcePhotometryModel (&source->extMag, source->modelEXT);
+	}
     }
 
@@ -160,11 +157,6 @@
     }
 
-    // replace source flux
-    // XXX need to be certain which model and size of mask for prior subtraction
-    // XXX full model or just analytical?
-    // XXX use pmSourceAdd instead?
-    if (source->tmpFlags & PM_SOURCE_TMPF_SUBTRACTED) {
-        pmModelAdd (source->pixels, source->maskObj, model, PM_MODEL_OP_FULL, maskVal);
-    }
+    // if we measure aperture magnitudes, the source must not currently be subtracted!
+    psAssert (!(source->tmpFlags & PM_SOURCE_TMPF_SUBTRACTED), "cannot measure ap mags if source is subtracted!");
 
     // if we are measuring aperture photometry and applying the growth correction,
@@ -201,9 +193,10 @@
     if (isfinite (source->apMag) && isPSF && psf) {
         if (psf->growth && (mode & PM_SOURCE_PHOT_GROWTH)) {
-            source->apMag += pmGrowthCurveCorrect (psf->growth, model->radiusFit);
+            source->apMag += pmGrowthCurveCorrect (psf->growth, source->apRadius);
         }
         if (mode & PM_SOURCE_PHOT_APCORR) {
+	    // XXX this should be removed -- we no longer fit for the 'sky bias'
             rflux   = pow (10.0, 0.4*source->psfMag);
-            source->apMag -= PS_SQR(model->radiusFit)*rflux * psf->skyBias + psf->skySat / rflux;
+            source->apMag -= PS_SQR(source->apRadius)*rflux * psf->skyBias + psf->skySat / rflux;
         }
     }
@@ -211,10 +204,4 @@
         psFree(flux);
         psFree(mask);
-    }
-
-    // if source was originally subtracted, re-subtract object, leave local sky
-    // XXX replace with pmSourceSub...
-    if (source->tmpFlags & PM_SOURCE_TMPF_SUBTRACTED) {
-        pmModelSub (source->pixels, source->maskObj, model, PM_MODEL_OP_FULL, maskVal);
     }
 
@@ -763,55 +750,2 @@
     return flux;
 }
-
-// XXX this is test code to verify the shift is doing the right thing (seems to be)
-# if (0)
-// measure centroid of unshifted gaussian (should be 16.0,16.0)
-        {
-	  psImage *image = source->pixels;
-	  float xo = 0.0;
-	  float yo = 0.0;
-	  float xo2 = 0.0;
-	  float yo2 = 0.0;
-	  float no = 0.0;
-	  for (int j = 0; j < image->numRows; j++)
-	  {
-	    for (int i = 0; i < image->numCols; i++) {
-	      xo += i*image->data.F32[j][i];
-	      yo += j*image->data.F32[j][i];
-	      xo2 += i*i*image->data.F32[j][i];
-	      yo2 += j*j*image->data.F32[j][i];
-	      no += image->data.F32[j][i];
-	    }
-	  }
-	  xo /= no;
-	  yo /= no;
-	  xo2 = sqrt (xo2/no - xo*xo);
-	  yo2 = sqrt (yo2/no - yo*yo);
-	  fprintf (stderr, "pre-shift centroid: %f,%f, sigma: %f,%f: flux: %f\n", xo, yo, xo2, yo2, no);
-        }
-
-// measure centroid of unshifted gaussian (should be 16.0,16.0)
-        {
-	  psImage *image = flux;
-	  float xo = 0.0;
-	  float yo = 0.0;
-	  float xo2 = 0.0;
-	  float yo2 = 0.0;
-	  float no = 0.0;
-	  for (int j = 0; j < image->numRows; j++)
-	  {
-	    for (int i = 0; i < image->numCols; i++) {
-	      xo += i*image->data.F32[j][i];
-	      yo += j*image->data.F32[j][i];
-	      xo2 += i*i*image->data.F32[j][i];
-	      yo2 += j*j*image->data.F32[j][i];
-	      no += image->data.F32[j][i];
-	    }
-	  }
-	  xo /= no;
-	  yo /= no;
-	  xo2 = sqrt (xo2/no - xo*xo);
-	  yo2 = sqrt (yo2/no - yo*yo);
-	  fprintf (stderr, "pre-shift centroid: %f,%f, sigma: %f,%f: flux: %f\n", xo, yo, xo2, yo2, no);
-        }
-# endif