Index: /branches/eam_branches/ipp-20130904/psphot/src/psphotExtendedSourceFits.c
===================================================================
--- /branches/eam_branches/ipp-20130904/psphot/src/psphotExtendedSourceFits.c	(revision 36317)
+++ /branches/eam_branches/ipp-20130904/psphot/src/psphotExtendedSourceFits.c	(revision 36318)
@@ -120,4 +120,9 @@
     }
 
+    float fitNsigmaConv = psMetadataLookupF32 (&status, recipe, "EXT_FIT_NSIGMA_CONV"); // number of sigma for the convolutio
+    if (!status || !isfinite(fitNsigmaConv) || fitNsigmaConv <= 0) {
+	fitNsigmaConv = 5.0;
+    }
+
     bool chisqConvergence = psMetadataLookupBool (&status, recipe, "LMM_FIT_CHISQ_CONVERGENCE"); // Fit tolerance
     if (!status) {
@@ -144,4 +149,5 @@
     fitOptions->minTol         = fitMinTol;
     fitOptions->maxTol         = fitMaxTol;
+    fitOptions->nsigma         = fitNsigmaConv;
 
     fitOptions->gainFactorMode   = gainFactorMode;
@@ -656,5 +662,5 @@
 	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);
Index: /branches/eam_branches/ipp-20130904/psphot/src/psphotGalaxyShape.c
===================================================================
--- /branches/eam_branches/ipp-20130904/psphot/src/psphotGalaxyShape.c	(revision 36317)
+++ /branches/eam_branches/ipp-20130904/psphot/src/psphotGalaxyShape.c	(revision 36318)
@@ -97,4 +97,9 @@
     assert (status);
 
+    float fitNsigmaConv = psMetadataLookupF32 (&status, recipe, "EXT_FIT_NSIGMA_CONV"); // number of sigma for the convolutio
+    if (!status || !isfinite(fitNsigmaConv) || fitNsigmaConv <= 0) {
+	fitNsigmaConv = 5.0;
+    }
+
     // Define source fitting parameters for extended source fits
     // we are not doing LMM fitting, so most options are irrelevant
@@ -102,4 +107,5 @@
     fitOptions->mode           = PM_SOURCE_FIT_EXT_AND_SKY;
     fitOptions->covarFactor    = psImageCovarianceFactorForAperture(readout->covariance, 10.0); // Covariance matrix
+    fitOptions->nsigma         = fitNsigmaConv;
 
     // Poisson or Constant weight for chisq tests?
@@ -298,5 +304,5 @@
 
         // cache the model flux
-        pmPCMCacheModel (source, maskVal, psfSize);
+        pmPCMCacheModel (source, maskVal, psfSize, fitOptions->nsigma);
     }
 
@@ -316,5 +322,5 @@
     // need to reset here each time since we assign below to calculate the flux
     // XXX note that this does not add sky to model
-    pmPCMMakeModel (source, pcm->modelConv, maskVal, psfSize);
+    pmPCMMakeModel (source, pcm->modelConv, pcm->nsigma, maskVal, psfSize);
     pcm->modelConv->params->data.F32[PM_PAR_I0] = 1.0;
 	
Index: /branches/eam_branches/ipp-20130904/psphot/src/psphotSourceFits.c
===================================================================
--- /branches/eam_branches/ipp-20130904/psphot/src/psphotSourceFits.c	(revision 36317)
+++ /branches/eam_branches/ipp-20130904/psphot/src/psphotSourceFits.c	(revision 36318)
@@ -653,8 +653,8 @@
     if (TIMING) { t4 = psTimerMark ("psphotFitPCM"); }
 
-    psphotPCMfitCheckSize (pcm, source, maskVal, psfSize);
-    if (pcm->modelConv->nIter == fitOptions->nIter) {
-	psphotPCMfitRetry (pcm, source, &options, maskVal, markVal, psfSize);
-    }
+    // psphotPCMfitCheckSize (pcm, source, maskVal, psfSize);
+    // if (pcm->modelConv->nIter == fitOptions->nIter) {
+    // 	psphotPCMfitRetry (pcm, source, &options, maskVal, markVal, psfSize);
+    // }
     if (TIMING) { t5 = psTimerMark ("psphotFitPCM"); }
 
@@ -787,10 +787,14 @@
 }
 
-// float indexGuessInv[] = {0.5, 0.33, 0.25, 0.167, 0.125, 0.083};
-float indexGuessInv[] = {0.5, 0.4, 0.3, 0.25, 0.20, 0.15, 0.125};
-# define N_INDEX_GUESS_INV 7
-
-// float reffGuess[] = {3.0, 10.0, 20.0, 30.0, 40.0};
-float reffGuess[] = {0.5, 0.75, 1.0, 1.4, 2.0};
+// 0.5 / n for (1.0, 1.25, 1.66, 2.0, 3.33, 4.0)
+// float indexGuessInv[] = {0.5, 0.4, 0.3, 0.25, 0.20, 0.15, 0.125};
+
+// 0.5 / n for (0.5, 1.0, 1.5, 2.0, 3.0, 4.0, 5.0, 6.0)
+float indexGuessInv[] = {1.00, 0.50, 0.333, 0.25, 1.666, 0.125, 0.10, 0.0833};
+float indexGuessR1q[] = {1.06, 1.19, 1.335, 1.48, 1.840, 2.290, 2.84, 3.5300};
+# define N_INDEX_GUESS_INV 8
+
+// we are going to guess in fractions about the R1-based guess
+float reffGuess[] = {0.8, 0.9, 1.0, 1.12, 1.25};
 # define N_REFF_GUESS 5
 
@@ -800,5 +804,5 @@
 
     // we get a reasonable guess from:
-    // * Reff = Kron R1
+    // * Reff = Kron R1 / Q(index) -- Q comes from Graham & Driver 
     // * Rmajor / Rminor & Theta from moments
     // * Io from total Kron flux
@@ -822,18 +826,4 @@
     psEllipseAxes momentAxes = psEllipseMomentsToAxes (moments, 20.0);
 
-    if (0) {
-	psEllipseAxes guessAxes;
-	guessAxes.major = source->moments->Mrf;
-	guessAxes.minor = (momentAxes.minor / momentAxes.minor) * guessAxes.major;
-	guessAxes.theta = momentAxes.theta;
-
-	if (!isfinite(guessAxes.major)) return false;
-	if (!isfinite(guessAxes.minor)) return false;
-	if (!isfinite(guessAxes.theta)) return false;
-
-	// convert the major,minor,theta to shape parameters for an Reff-like model
-	pmModelAxesToParams (&PAR[PM_PAR_SXX], &PAR[PM_PAR_SXY], &PAR[PM_PAR_SYY], guessAxes, true);
-    }
-
     // set the model position
     if (!pmModelSetPosition(&PAR[PM_PAR_XPOS], &PAR[PM_PAR_YPOS], source)) {
@@ -855,23 +845,22 @@
     // loop over index guesses and find the best fit
     for (int j = 0; j < N_REFF_GUESS; j++) {
-
-	psEllipseAxes guessAxes;
-	guessAxes.major = reffGuess[j] * source->moments->Mrf;
-	guessAxes.minor = guessAxes.major * (momentAxes.minor / momentAxes.major);
-	guessAxes.theta = momentAxes.theta;
-
-	if (!isfinite(guessAxes.major)) return false;
-	if (!isfinite(guessAxes.minor)) return false;
-	if (!isfinite(guessAxes.theta)) return false;
-
-	// convert the major,minor,theta to shape parameters for an Reff-like model
-	pmModelAxesToParams (&PAR[PM_PAR_SXX], &PAR[PM_PAR_SXY], &PAR[PM_PAR_SYY], guessAxes, true);
-
 	for (int i = 0; i < N_INDEX_GUESS_INV; i++) {
 	    PAR[PM_PAR_7] = indexGuessInv[i];
 
+	    psEllipseAxes guessAxes;
+	    guessAxes.major = reffGuess[j] * source->moments->Mrf / indexGuessR1q[i];
+	    guessAxes.minor = guessAxes.major * (momentAxes.minor / momentAxes.major);
+	    guessAxes.theta = momentAxes.theta;
+
+	    if (!isfinite(guessAxes.major)) return false;
+	    if (!isfinite(guessAxes.minor)) return false;
+	    if (!isfinite(guessAxes.theta)) return false;
+
+	    // convert the major,minor,theta to shape parameters for an Reff-like model
+	    pmModelAxesToParams (&PAR[PM_PAR_SXX], &PAR[PM_PAR_SXY], &PAR[PM_PAR_SYY], guessAxes, true);
+
 	    // generated the modelFlux
 	    // XXX note that this does not add sky to model
-	    pmPCMMakeModel (source, pcm->modelConv, maskVal, psfSize);
+	    pmPCMMakeModel (source, pcm->modelConv, pcm->nsigma, maskVal, psfSize);
 	
 	    float YY = 0.0;
@@ -911,5 +900,5 @@
 		iMin = Io;
 		sMin = indexGuessInv[i];
-		rMin = reffGuess[j];
+		rMin = reffGuess[j] / indexGuessR1q[i];
 	    }
 	    // fprintf (stderr, "%d | %f %f %f %f | %f %f %f %f", i, indexGuessInv[i], reffGuess[j], Io, Chisq, sMin, rMin, iMin, xMin);
@@ -1023,5 +1012,5 @@
     psVectorAppend (chi2, model->chisqNorm);
 
-    PAR[PM_PAR_7] = (nStart < N_INDEX_GUESS_INV - 1) ? 0.5*(indexGuessInv[nStart + 1] + indexGuessInv[nStart]) : 0.1;
+    PAR[PM_PAR_7] = (nStart < N_INDEX_GUESS_INV - 1) ? 0.5*(indexGuessInv[nStart + 1] + indexGuessInv[nStart]) : indexGuessInv[nStart - 1];
     pmSourceFitPCM (pcm, source, &options, maskVal, markVal, psfSize);
     if (EXTRA_VERBOSE) fprintf (stderr, "%d >>> %d %f : %f - %f %f - %f %f %f - %f\n", source->id, model->nIter, model->chisqNorm, PAR[7], PAR[2], PAR[3], PAR[4], PAR[5], PAR[6], PAR[1]);
@@ -1029,23 +1018,9 @@
     psVectorAppend (chi2, model->chisqNorm);
 
-    PAR[PM_PAR_7] = (nStart > 0) ? 0.5*(indexGuessInv[nStart - 1] + indexGuessInv[nStart]) : 0.55;
+    PAR[PM_PAR_7] = (nStart > 0) ? 0.5*(indexGuessInv[nStart - 1] + indexGuessInv[nStart]) : indexGuessInv[nStart + 1];
     pmSourceFitPCM (pcm, source, &options, maskVal, markVal, psfSize);
     if (EXTRA_VERBOSE) fprintf (stderr, "%d >>> %d %f : %f - %f %f - %f %f %f - %f\n", source->id, model->nIter, model->chisqNorm, PAR[7], PAR[2], PAR[3], PAR[4], PAR[5], PAR[6], PAR[1]);
     psVectorAppend (Sidx, 100*PAR[PM_PAR_7]);
     psVectorAppend (chi2, model->chisqNorm);
-
-    if (chi2->data.F32[1] < chi2->data.F32[2]) {
-      if (nStart == N_INDEX_GUESS_INV - 1) {
-	PAR[PM_PAR_7] = 0.11;
-      } else {
-	PAR[PM_PAR_7] = indexGuessInv[nStart + 1];
-      }
-    } else {
-      if (nStart == 0) {
-	PAR[PM_PAR_7] = 0.52;
-      } else {
-	PAR[PM_PAR_7] = indexGuessInv[nStart - 1];
-      }
-    }
 
     psPolynomial1D *poly = psPolynomial1DAlloc (PS_POLYNOMIAL_ORD, 2);
@@ -1292,5 +1267,5 @@
 	// generated the modelFlux
 	    // XXX note that this does not add sky to model
-	pmPCMMakeModel (source, pcm->modelConv, maskVal, psfSize);
+	pmPCMMakeModel (source, pcm->modelConv, pcm->nsigma, maskVal, psfSize);
 	
 	float YY = 0.0;
@@ -1331,5 +1306,5 @@
 	    rMin = dref;
 	}
-	// fprintf (stderr, "%d | %f %f %f | %f %f %f\n", j, dref, Io, Chisq, rMin, iMin, xMin);
+	fprintf (stderr, "%d | %f %f %f | %f %f %f\n", j, dref, Io, Chisq, rMin, iMin, xMin);
     }
 
@@ -1375,5 +1350,5 @@
     // generated the modelFlux
 	    // XXX note that this does not add sky to model
-    pmPCMMakeModel (source, pcm->modelConv, maskVal, psfSize);
+    pmPCMMakeModel (source, pcm->modelConv, pcm->nsigma, maskVal, psfSize);
 	
     float YY = 0.0;