Index: trunk/psphot/src/psphotExtendedSourceFits.c
===================================================================
--- trunk/psphot/src/psphotExtendedSourceFits.c	(revision 36108)
+++ trunk/psphot/src/psphotExtendedSourceFits.c	(revision 36375)
@@ -65,6 +65,4 @@
     psphotVisualShowImage(readout);
 
-    // psphotSaveImage (NULL, readout->image, "test.01.fits");
-
     pmDetections *detections = psMetadataLookupPtr (&status, readout->analysis, "PSPHOT.DETECTIONS");
     psAssert (detections, "missing detections?");
@@ -111,4 +109,9 @@
     if (!status || !isfinite(fitMaxTol) || fitMaxTol <= 0) {
 	fitMaxTol = 1.0;
+    }
+
+    float fitNsigmaConv = psMetadataLookupF32 (&status, recipe, "EXT_FIT_NSIGMA_CONV"); // number of sigma for the convolution
+    if (!status || !isfinite(fitNsigmaConv) || fitNsigmaConv <= 0) {
+	fitNsigmaConv = 5.0;
     }
 
@@ -137,8 +140,13 @@
     fitOptions->minTol         = fitMinTol;
     fitOptions->maxTol         = fitMaxTol;
+    fitOptions->nsigma         = fitNsigmaConv;
 
     fitOptions->gainFactorMode   = gainFactorMode;
     fitOptions->chisqConvergence = chisqConvergence;
     fitOptions->isInteractive    = isInteractive;
+
+    // use poissonian errors or local-sky errors
+    fitOptions->poissonErrors = psMetadataLookupBool (&status, recipe, "EXTENDED_SOURCE_FITS_POISSON");
+    if (!status) fitOptions->poissonErrors = true;
 
     // maskVal is used to test for rejected pixels, and must include markVal
@@ -285,4 +293,5 @@
             psArrayAdd(job->args, 1, cells->data[j]); // sources
             psArrayAdd(job->args, 1, models);
+
             // Allocate a metadata iterator here because psMetadataIteratorAlloc/Free are not thread safe
             psMetadataIterator *iter = psMetadataIteratorAlloc (models, PS_LIST_HEAD, NULL);
@@ -386,6 +395,4 @@
     psphotSersicModelClassCleanup();
 
-    // psphotSaveImage (NULL, readout->image, "test.02.fits");
-
     psphotVisualShowResidualImage (readout, false);
 
@@ -428,4 +435,12 @@
     // psTraceSetLevel ("psLib.math.psMinimizeLMChi2_Alt", 5);
 
+    pmModelStatus badModel = PM_MODEL_STATUS_NONE;
+    badModel |= PM_MODEL_STATUS_BADARGS;
+    badModel |= PM_MODEL_STATUS_OFFIMAGE;
+    badModel |= PM_MODEL_STATUS_NAN_CHISQ;
+    badModel |= PM_MODEL_SERSIC_PCM_FAIL_GUESS;
+    badModel |= PM_MODEL_SERSIC_PCM_FAIL_GRID;
+    badModel |= PM_MODEL_PCM_FAIL_GUESS;
+
     // choose the sources of interest
     for (int i = 0; i < sources->n; i++) {
@@ -433,23 +448,13 @@
         pmSource *source = sources->data[i];
 
-        // skip PSF-like and non-astronomical objects
-        if (!(source->mode & PM_SOURCE_MODE_EXT_LIMIT)) continue;
-        if (source->type == PM_SOURCE_TYPE_DEFECT) continue;
-        if (source->type == PM_SOURCE_TYPE_SATURATED) continue;
-
-	// skip saturated stars modeled with a radial profile 
-        if (source->mode2 & PM_SOURCE_MODE2_SATSTAR_PROFILE) continue;
-
-        // XXX this should use peak?
+	// rules for measuring petrosian parameters for specific objects are set in
+	// psphotChooseAnalysisOptions.c
+	if (!(source->tmpFlags & PM_SOURCE_TMPF_EXT_FIT)) continue;
+
+	// limit selection by analysis region (XXX move this into psphotChooseAnalysisOption?)
         if (source->peak->x < region->x0) continue;
         if (source->peak->y < region->y0) continue;
         if (source->peak->x > region->x1) continue;
         if (source->peak->y > region->y1) continue;
-
-
-	// XXX for a test, just do the obvious trail
-	// XXX if (source->peak->xf < 1100) continue;
-	// XXX if (source->peak->xf > 1400) continue;
-	// XXX if (source->peak->yf >  245) continue;
 
         // replace object in image
@@ -462,7 +467,4 @@
 	psphotSetRadiusMomentsExact(&fitRadius, &windowRadius, readout, source, markVal); // NOTE : 6 allocs
 
-	// XXX WATCH OUT HERE!!
-	// fitRadius = 30;
-
 	// UPDATE : we have changed the moments calculation.  There is now an iteration within 
 	// psphotKronMasked to determine moments appropriate for a larger object.  The values
@@ -475,10 +477,7 @@
         }
 
+# ifdef TEST_OBJECT
 	bool testObject = false;
-	// testObject |= ((fabs(source->peak->xf -  179) < 5) && (fabs(source->peak->yf - 1138) < 5));
-	// testObject |= ((fabs(source->peak->xf - 5047) < 5) && (fabs(source->peak->yf -  151) < 5));
-	// testObject |= ((fabs(source->peak->xf - 3929) < 5) && (fabs(source->peak->yf - 4109) < 5));
-	// testObject |= ((fabs(source->peak->xf -  915) < 5) && (fabs(source->peak->yf - 5998) < 5));
-	// testObject |= ((fabs(source->peak->xf - 5406) < 5) && (fabs(source->peak->yf -  326) < 5));
+	testObject |= ((fabs(source->peak->xf -  179) < 5) && (fabs(source->peak->yf - 1138) < 5));
 	if (testObject) {
 	    fprintf (stderr, "test object @ %f, %f\n", source->peak->xf, source->peak->yf);
@@ -486,4 +485,5 @@
 	    psTraceSetLevel ("psphot.psphotExtendedSourceFits_Threaded", 5);
 	}
+# endif 
 
         // loop here over the models chosen for each source (exclude by S/N)
@@ -501,15 +501,13 @@
           assert (status);
 
-	  // limit selection to some SN limit
-	  bool skipSource = false;
-	  if (source->mode & PM_SOURCE_MODE_EXT_LIMIT) {
-	      skipSource = (source->moments->KronFlux < FIT_SN_LIM * source->moments->KronFluxErr);
-	  } else {
-	      skipSource = (sqrt(source->peak->detValue) < FIT_SN_LIM);
-	  }
-          if (skipSource) {
+	  // limit selection to some SN limit for specific models (this value only applies if > EXTENDED_SOURCE_SN_LIM)
+	  if (isfinite(FIT_SN_LIM)) {
+	    if (source->moments->KronFlux < FIT_SN_LIM * source->moments->KronFluxErr) {
 	      Nfaint ++;
 	      continue;
+	    }
 	  }
+
+	  source->mode2 |= PM_SOURCE_MODE2_EXT_FITS_RUN;
 
           // check on the model type
@@ -529,4 +527,5 @@
                   psTrace ("psphot", 5, "failed to fit psf-conv model for object at %f, %f", source->moments->Mx, source->moments->My);
 		  Nfail ++;
+		  source->mode2 |= PM_SOURCE_MODE2_EXT_FITS_FAIL;
                   continue;
               }
@@ -534,5 +533,5 @@
 		       source->moments->Mx, source->moments->My, pmModelClassGetName (modelFit->type), modelFit->chisq, modelFit->nPix, modelFit->nIter);
               Nconvolve ++;
-              if (!(modelFit->flags & (PM_MODEL_STATUS_BADARGS | PM_MODEL_STATUS_NONCONVERGE | PM_MODEL_STATUS_OFFIMAGE))) {
+              if (!(modelFit->flags & badModel)) {
                   NconvolvePass ++;
 		  source->mode |= PM_SOURCE_MODE_EXTENDED_FIT;
@@ -548,9 +547,10 @@
 		      Nfail ++;
 		      doneFits = true;
+		      source->mode2 |= PM_SOURCE_MODE2_EXT_FITS_FAIL;
 		      continue;
 		  }
 		  psTrace ("psphot", 4, "fit plain model for %f, %f : %s chisq = %f (npix: %d, niter: %d)\n", source->moments->Mx, source->moments->My, pmModelClassGetName (modelFit->type), modelFit->chisq, modelFit->nPix, modelFit->nIter);
 		  Nplain ++;
-		  if (!(modelFit->flags & (PM_MODEL_STATUS_BADARGS | PM_MODEL_STATUS_NONCONVERGE | PM_MODEL_STATUS_OFFIMAGE))) {
+		  if (!(modelFit->flags & badModel)) {
 		      NplainPass ++;
 		      source->mode |= PM_SOURCE_MODE_EXTENDED_FIT;
@@ -576,10 +576,10 @@
 			  fprintf (stderr, "update window : %f %f : %f -> %f\n", source->peak->xf, source->peak->yf, fitRadius, 2*fitRadius);
 			  psphotSetWindowTrail (&fitRadius, &windowRadius, readout, source, markVal, fitRadius*2.0);
+			  source->mode2 |= PM_SOURCE_MODE2_EXT_FITS_RETRY;
 		      }
 		  }
 	      }
           }
-	  // XXX really need to do this in a cleaner way:
-	  if (!modelFit) continue;
+	  psAssert (modelFit, "modelFit not set?");
 
           // test for fit quality / result
@@ -590,5 +590,16 @@
         }
 
+	// we are allowed to fit both stars and non-stars here -- if we have fitted
+	// something which we think is a star, we should use that model to subtract the
+	// object from the image.
+        if (source->type == PM_SOURCE_TYPE_STAR) {
+	  // ensure the modelPSF is cached
+	  pmSourceCacheModel (source, maskVal);
+          pmSourceSub (source, PM_MODEL_OP_FULL, maskVal);
+          continue;
+        }
+
         // evaluate the relative quality of the models, choose one
+	// the PSF model might be the best fit : allow it to succeed
         float minChisq = NAN;
         int minModel = -1;
@@ -596,9 +607,9 @@
             pmModel *model = source->modelFits->data[i];
 
+	    // skip the really bad fits
             if (!(model->flags & PM_MODEL_STATUS_FITTED)) continue;
-
-            if (model->flags & (PM_MODEL_STATUS_BADARGS)) continue;
+            if (model->flags & badModel) continue;
+
             // if (model->flags & (PM_MODEL_STATUS_NONCONVERGE)) continue;
-            if (model->flags & (PM_MODEL_STATUS_OFFIMAGE)) continue;
 
             if ((minModel < 0) || (model->chisq < minChisq)) {
@@ -618,4 +629,5 @@
 	  pmSourceCacheModel (source, maskVal);
 
+# if (PS_TRACE_ON)
 	  pmModel *model = source->modelFits->data[0];
 	  int flags = 0xffffffff;
@@ -623,8 +635,9 @@
 	    flags = model->flags;
 	  }
-
           fprintf (stderr, "failed to fit extended source model to object %d @ %f, %f (%x)\n", source->id, source->moments->Mx, source->moments->My, flags);
+#endif
           pmSourceSub (source, PM_MODEL_OP_FULL, maskVal);
 
+	  source->mode2 |= PM_SOURCE_MODE2_EXT_FITS_NONE;
           continue;
         }
@@ -643,9 +656,10 @@
 	if (source->modelEXT->isPCM) {
 	    // fprintf (stderr, "subtract PCM extended source model for object %d @ %f, %f\n", source->id, source->moments->Mx, source->moments->My);
-	    pmPCMCacheModel (source, maskVal, psfSize);
+	    pmPCMCacheModel (source, maskVal, psfSize, fitOptions->nsigma);
 	} else {
 	    // fprintf (stderr, "subtract non-PCM extended source model for object %d @ %f, %f\n", source->id, source->moments->Mx, source->moments->My);
 	    pmSourceCacheModel (source, maskVal);
 	}
+        source->modelEXT->flags |= PM_MODEL_BEST_FIT;
 
         // subtract the best fit from the object, leave local sky
@@ -655,8 +669,10 @@
         psTrace ("psphot", 5, "extended source model for source at %7.1f, %7.1f", source->moments->Mx, source->moments->My);
 
+# ifdef TEST_OBJECT
 	if (testObject) {
 	    psTraceSetLevel ("psModules.objects.pmPCM_MinimizeChisq", 0);
 	    psTraceSetLevel ("psphot.psphotExtendedSourceFits_Threaded", 0);
 	}
+# endif
     }