Index: trunk/psphot/src/psphotSourceFits.c
===================================================================
--- trunk/psphot/src/psphotSourceFits.c	(revision 31452)
+++ trunk/psphot/src/psphotSourceFits.c	(revision 32348)
@@ -211,5 +211,5 @@
 bool psphotFitBlob (pmReadout *readout, pmSource *source, psArray *newSources, pmPSF *psf, pmSourceFitOptions *fitOptions, psImageMaskType maskVal, psImageMaskType markVal) {
 
-    float radius;
+    float fitRadius, windowRadius;
     bool okEXT, okDBL;
     pmModel *ONE = NULL;
@@ -217,5 +217,5 @@
     pmModel *EXT = NULL;
     psArray *DBL = NULL;
-    pmMoments psfMoments;
+    // pmMoments psfMoments;
 
     // skip the source if we don't think it is extended
@@ -225,18 +225,17 @@
     if (source->type == PM_SOURCE_TYPE_SATURATED) return false;
 
+# define TEST_X -420.0
+# define TEST_Y 300.0
+   
+    if ((fabs(source->peak->xf - TEST_X) < 5) && (fabs(source->peak->yf - TEST_Y) < 5)) {
+	fprintf (stderr, "test galaxy\n");
+    }
+
+# undef TEST_X
+# undef TEST_Y
+
     // set the radius based on the footprint (also sets the mask pixels)
-    if (!psphotSetRadiusFootprint(&radius, readout, source, markVal, 1.0)) return false;
-
-    // XXX note that this changes the source moments that are published...
-    // XXX all published moments should use the same measurement
-    // recalculate the source moments using the larger extended-source moments radius
-    // at this stage, skip Gaussian windowing, and do not clip pixels by S/N
-    // this uses the footprint to judge both radius and aperture?
-    // XXX save the psf-based moments for output
-    psfMoments = *source->moments;
-    if (!pmSourceMoments (source, radius, 0.0, 0.5, 0.0, maskVal)) {
-      *source->moments = psfMoments;
-      return false;
-    }
+    if (!psphotSetRadiusMoments(&fitRadius, &windowRadius, readout, source, markVal)) return false;
+    // fprintf (stderr, "rad: %6.1f %6.1f  | %5.2f %5.2f %5.2f  ", source->peak->xf, source->peak->yf, source->moments->Mrf, fitRadius, windowRadius);
 
     psTrace ("psphot", 5, "trying blob...\n");
@@ -263,7 +262,7 @@
 	ONE = DBL->data[0];
 	if (ONE) {
-	    if (!isfinite(ONE->params->data.F32[PM_PAR_I0])) psAbort("nan in fit");
+	    psAssert (isfinite(ONE->params->data.F32[PM_PAR_I0]), "nan in fit");
 	    chiDBL = ONE->chisqNorm; // save chisq for double-star/galaxy comparison
-	    ONE->fitRadius = radius;
+	    ONE->fitRadius = fitRadius;
 	}
 
@@ -271,6 +270,6 @@
 	ONE = DBL->data[1];
 	if (ONE) {
-	    if (!isfinite(ONE->params->data.F32[PM_PAR_I0])) psAbort("nan in fit");
-	    ONE->fitRadius = radius;
+	    psAssert (isfinite(ONE->params->data.F32[PM_PAR_I0]), "nan in fit");
+	    ONE->fitRadius = fitRadius;
 	}
     }
@@ -284,4 +283,5 @@
 	okEXT = psphotEvalEXT (tmpSrc, EXT);
 	chiEXT = EXT ? EXT->chisqNorm : NAN;
+	EXT->fitRadius = fitRadius;
     }
 
@@ -294,5 +294,4 @@
 
     if (okEXT && okDBL) {
-        psTrace ("psphot", 5, "blob chisq: %f vs %f\n", chiEXT, chiDBL);
         // XXX EAM : a bogus bias: need to examine this better
         if (3*chiEXT > chiDBL) goto keepDBL;
@@ -303,8 +302,10 @@
     if (!okEXT && okDBL) goto keepDBL;
 
+    psTrace ("psphot", 4, "both failed: blob chisq: %f vs %f for %f,%f\n", chiEXT, chiDBL, source->peak->xf, source->peak->yf);
+
     // both models failed; reject them both
     // XXX -- change type flags to psf in this case, and make sure we subtract it?
     // reset the psf moments
-    *source->moments = psfMoments;
+    // XXX *source->moments = psfMoments;
 
     psFree (EXT);
@@ -313,4 +314,5 @@
 
 keepEXT:
+    psTrace ("psphot", 4, "goto EXT : blob chisq: %f vs %f for %f,%f\n", chiEXT, chiDBL, source->peak->xf, source->peak->yf);
     // sub EXT
     psFree (DBL);
@@ -338,8 +340,9 @@
 
     // reset the psf moments
-    *source->moments = psfMoments;
+    // XXX *source->moments = psfMoments;
     return true;
 
 keepDBL:
+    psTrace ("psphot", 4, "goto DBL : blob chisq: %f vs %f for %f,%f\n", chiEXT, chiDBL, source->peak->xf, source->peak->yf);
     // sub DLB
     psFree (EXT);
@@ -372,11 +375,6 @@
 	    psLogMsg ("psphot", 1, "PAR %d : %f +/- %f\n", i, source->modelPSF->params->data.F32[i], source->modelPSF->dparams->data.F32[i]);
 	}
-	psphotVisualShowResidualImage (readout, false);
     }
 # endif
-
-    // reset the (original) psf moments
-    *source->moments = psfMoments;
-    *newSrc->moments = psfMoments;
 
     psArrayAdd (newSources, 100, newSrc);
@@ -386,6 +384,4 @@
 
 escape:
-    // reset the psf moments
-    *source->moments = psfMoments;
     psImageMaskPixels (source->maskObj, "AND", PS_NOT_IMAGE_MASK(markVal)); 
     psFree (tmpSrc);
@@ -476,11 +472,5 @@
     // for sersic models, use a grid search to choose an index, then float the params there
     if (modelType == pmModelClassGetType("PS_MODEL_SERSIC")) {
-	// for the test fits, use a somewhat smaller radius
-	psphotSetRadiusFootprint(&model->fitRadius, readout, source, markVal, 0.5);
     	psphotFitSersicIndex (model, readout, source, fitOptions, maskVal, markVal);
-    }
-
-    if (!psphotSetRadiusModel (model, readout, source, markVal, true)) {
-	psphotSetRadiusFootprint(&model->fitRadius, readout, source, markVal, 1.0);
     }
 
@@ -494,8 +484,10 @@
     pmSourceFitModel (source, model, &options, maskVal);
     // fprintf (stderr, "chisq: %f, nIter: %d, radius: %f, npix: %d\n", model->chisqNorm, model->nIter, model->fitRadius, model->nPix);
-
     // psTraceSetLevel("psLib.math.psMinimizeLMChi2", 0);
+
     return (model);
 }
+
+# define TIMING 0
 
 pmModel *psphotFitPCM (pmReadout *readout, pmSource *source, pmSourceFitOptions *fitOptions, pmModelType modelType, psImageMaskType maskVal, psImageMaskType markVal, int psfSize) {
@@ -517,4 +509,14 @@
 	return NULL;
     }
+
+# define TEST_X -540.0
+# define TEST_Y 540.0
+   
+    if ((fabs(source->peak->xf - TEST_X) < 5) && (fabs(source->peak->yf - TEST_Y) < 5)) {
+	psTraceSetLevel("psModules.objects.pmPCM_MinimizeChisq", 5);
+    }
+
+    float t1, t2, t4, t5;
+    if (TIMING) { psTimerStart ("psphotFitPCM"); }
 
     pmPCMdata *pcm = pmPCMinit (source, &options, model, maskVal, psfSize);
@@ -524,27 +526,24 @@
 	return model;
     }
-
-    // use the source moments, etc to guess basic model parameters
-    if (!pmSourceModelGuessPCM (pcm, source, maskVal, markVal)) {
-	psFree (pcm);
-        model->flags |= PM_MODEL_STATUS_BADARGS;
-	return model;
-    }
-
-    // for sersic models, use a grid search to choose an index, then float the params there
+    if (TIMING) { t1 = psTimerMark ("psphotFitPCM"); }
+
+    // get the guess for sersic models 
     if (modelType == pmModelClassGetType("PS_MODEL_SERSIC")) {
-	// for the test fits, use a somewhat smaller radius
-	psphotSetRadiusFootprint(&model->fitRadius, readout, source, markVal, 0.5);
-
-    	if (!psphotFitSersicIndexPCM (pcm, readout, source, fitOptions, maskVal, markVal, psfSize)) {
+	// use the source moments, etc to guess basic model parameters
+	if (!psphotSersicModelClassGuessPCM (pcm, source)) {
 	    psFree (pcm);
 	    model->flags |= PM_MODEL_STATUS_BADARGS;
 	    return model;
 	}
-    }
-
-    if (!psphotSetRadiusModel (model, readout, source, markVal, true)) {
-	psphotSetRadiusFootprint(&model->fitRadius, readout, source, markVal, 1.0);
-    }
+    } else {
+	// use the source moments, etc to guess basic model parameters
+	if (!pmSourceModelGuessPCM (pcm, source, maskVal, markVal)) {
+	    psFree (pcm);
+	    model->flags |= PM_MODEL_STATUS_BADARGS;
+	    return model;
+	}
+    }
+
+    if (TIMING) { t2 = psTimerMark ("psphotFitPCM"); }
 
     if (modelType == pmModelClassGetType("PS_MODEL_SERSIC")) {
@@ -555,8 +554,18 @@
     // update the pcm elements if we have changed the circumstance (options.mode or source->pixels)
     pmPCMupdate(pcm, source, &options, model);
+    if (TIMING) { t4 = psTimerMark ("psphotFitPCM"); }
 
     // psTraceSetLevel("psLib.math.psMinimizeLMChi2", 5);
     pmSourceFitPCM (pcm, source, &options, maskVal, markVal, psfSize);
-    // fprintf (stderr, "chisq: %f, nIter: %d, radius: %f, npix: %d\n", model->chisqNorm, model->nIter, model->fitRadius, model->nPix);
+    if (TIMING) { t5 = psTimerMark ("psphotFitPCM"); }
+
+    if (TIMING) {
+	int nPixBig = source->pixels->numCols * source->pixels->numRows;
+    	fprintf (stderr, "psphotFitPCM : nIter: %2d, radius: %6.1f, npix: %5d of %5d, t1: %6.4f, t2: %6.4f, t4: %6.4f, t5: %6.4f\n", model->nIter, model->fitRadius, model->nPix, nPixBig, t1, t2, t4, t5);
+    }
+
+    if ((fabs(source->peak->xf - TEST_X) < 5) && (fabs(source->peak->yf - TEST_Y) < 5)) {
+	psTraceSetLevel("psModules.objects.pmPCM_MinimizeChisq", 0);
+    }
 
     // psTraceSetLevel("psLib.math.psMinimizeLMChi2", 0);
@@ -566,7 +575,12 @@
 }
 
+# undef TEST_X
+# undef TEST_Y
+
 // note that these should be 1/2n of the standard sersic index
-float indexGuess[] = {0.5, 0.33, 0.25, 0.167, 0.125, 0.083};
-# define N_INDEX_GUESS 6
+// float indexGuess[] = {0.8, 1.0, 1.5, 2.0, 2.5, 3.0, 3.5, 4.0, 4.5, 5.0, 5.5, 6.0};
+// float indexGuess[] = {0.5, 0.33, 0.25, 0.167, 0.125, 0.083};
+float indexGuess[] = {1.0, 2.0, 3.0, 4.0};
+# define N_INDEX_GUESS 4
 
 // A sersic model is very sensitive to the index.  attempt to find the index first by grid search in just the index
@@ -586,6 +600,13 @@
     float chiSquare[N_INDEX_GUESS];
 
+# define TEST_X -540.0
+# define TEST_Y 540.0
+   
+    if ((fabs(source->peak->xf - TEST_X) < 5) && (fabs(source->peak->yf - TEST_Y) < 5)) {
+	psTraceSetLevel("psModules.objects.pmPCM_MinimizeChisq", 5);
+    }
+
     for (int i = 0; i < N_INDEX_GUESS; i++) {
-	model->params->data.F32[PM_PAR_7] = indexGuess[i];
+	model->params->data.F32[PM_PAR_7] = 0.5/indexGuess[i];
 
 	if (!model->modelGuess(model, source)) {
@@ -594,12 +615,6 @@
 	}
 
-	// each time we change the model guess, we need to adjust the radius
-	// XXX this did not work : we do not need such a large radius -- just uses moments-based radius
-	if (false && !psphotSetRadiusModel (model, readout, source, markVal, false)) {
-	    psphotSetRadiusFootprint(&model->fitRadius, readout, source, markVal, 0.5);
-	}
-	
 	pmSourceFitModel (source, model, &options, maskVal);
-	// fprintf (stderr, "chisq: %f, nIter: %d, radius: %f, npix: %d\n", model->chisqNorm, model->nIter, model->fitRadius, model->nPix);
+	// fprintf (stderr, "index: %f, chisq: %f, nIter: %d, radius: %f, npix: %d\n", indexGuess[i], model->chisqNorm, model->nIter, model->fitRadius, model->nPix);
 
 	chiSquare[i] = model->chisqNorm;
@@ -616,12 +631,11 @@
     assert (iMin >= 0);
 
+    if ((fabs(source->peak->xf - TEST_X) < 5) && (fabs(source->peak->yf - TEST_Y) < 5)) {
+	psTraceSetLevel("psModules.objects.pmPCM_MinimizeChisq", 0);
+    }
+
     model->flags = PM_MODEL_STATUS_NONE; // do not attempt to handle failures here, let the next iteration deal with it
-    model->params->data.F32[PM_PAR_7] = indexGuess[iMin];
+    model->params->data.F32[PM_PAR_7] = 0.5/indexGuess[iMin];
     model->modelGuess(model, source);
-
-    // each time we change the model guess, we need to adjust the radius
-    // if (!psphotSetRadiusModel (model, readout, source, markVal, true)) {
-    // 	psphotSetRadiusFootprint(&model->fitRadius, readout, source, markVal);
-    // }
 
     return true;
@@ -648,4 +662,8 @@
     float xMin = NAN;
     float chiSquare[N_INDEX_GUESS];
+
+    if ((fabs(source->peak->xf - TEST_X) < 5) && (fabs(source->peak->yf - TEST_Y) < 5)) {
+	psTraceSetLevel("psModules.objects.pmPCM_MinimizeChisq", 5);
+    }
 
     for (int i = 0; i < N_INDEX_GUESS; i++) {
@@ -657,11 +675,6 @@
 	}
 
-	// each time we change the model guess, we need to adjust the radius
-	// XXX this did not work : we do not need such a large radius -- just uses moments-based radius
-	if (false && !psphotSetRadiusModel (model, readout, source, markVal, false)) {
-	    psphotSetRadiusFootprint(&model->fitRadius, readout, source, markVal, 0.5);
-	}
-	
 # if (0)
+	// this block is to test the relative speed of straight and PCM fits
 	pmSourceFitModel (source, model, &options, maskVal);
 # else
@@ -669,4 +682,5 @@
 	pmSourceFitPCM (pcm, source, &options, maskVal, markVal, psfSize);
 # endif
+	fprintf (stderr, "index: %f, chisq: %f, nIter: %d, radius: %f, npix: %d\n", indexGuess[i], model->chisqNorm, model->nIter, model->fitRadius, model->nPix);
 	// fprintf (stderr, "chisq: %f, nIter: %d, radius: %f, npix: %d\n", model->chisqNorm, model->nIter, model->fitRadius, model->nPix);
 
@@ -684,4 +698,8 @@
     assert (iMin >= 0);
     
+    if ((fabs(source->peak->xf - TEST_X) < 5) && (fabs(source->peak->yf - TEST_Y) < 5)) {
+	psTraceSetLevel("psModules.objects.pmPCM_MinimizeChisq", 0);
+    }
+
     model->flags = PM_MODEL_STATUS_NONE; // do not attempt to handle failures here, let the next iteration deal with it
     model->params->data.F32[PM_PAR_7] = indexGuess[iMin];