Index: branches/eam_branches/ipp-20130904/psphot/src/psphotSourceFits.c
===================================================================
--- branches/eam_branches/ipp-20130904/psphot/src/psphotSourceFits.c	(revision 36318)
+++ branches/eam_branches/ipp-20130904/psphot/src/psphotSourceFits.c	(revision 36325)
@@ -653,5 +653,5 @@
     if (TIMING) { t4 = psTimerMark ("psphotFitPCM"); }
 
-    // psphotPCMfitCheckSize (pcm, source, maskVal, psfSize);
+    if (options.isInteractive) psphotPCMfitCheckSize (pcm, source, maskVal, psfSize);
     // if (pcm->modelConv->nIter == fitOptions->nIter) {
     // 	psphotPCMfitRetry (pcm, source, &options, maskVal, markVal, psfSize);
@@ -791,5 +791,5 @@
 
 // 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 indexGuessInv[] = {1.00, 0.50, 0.333, 0.25, 0.166, 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
@@ -902,6 +902,6 @@
 		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);
-	    // fprintf (stderr, "\n");
+	    fprintf (stderr, "%d | %f %f %f %f | %f %f %f %f", i, indexGuessInv[i], reffGuess[j], Io, Chisq, sMin, rMin, iMin, xMin);
+	    fprintf (stderr, "\n");
 	}
     }
@@ -1006,5 +1006,20 @@
     psVector *Sidx = psVectorAllocEmpty (16, PS_TYPE_F32);
 
-    PAR[PM_PAR_7] = indexGuessInv[nStart];
+    float Sm = NAN, Sp = NAN, So = NAN;
+    if (nStart == 0) {
+	Sm = indexGuessInv[nStart];
+	So = 0.5*(indexGuessInv[nStart + 1] + indexGuessInv[nStart]);
+	Sp = indexGuessInv[nStart + 1];
+    } else if (nStart == N_INDEX_GUESS_INV - 1) {
+	Sp = indexGuessInv[nStart];
+	So = 0.5*(indexGuessInv[nStart - 1] + indexGuessInv[nStart]);
+	Sm = indexGuessInv[nStart - 1];
+    } else {
+	Sm = 0.5*(indexGuessInv[nStart - 1] + indexGuessInv[nStart]);
+	So = indexGuessInv[nStart];
+	Sp = 0.5*(indexGuessInv[nStart + 1] + indexGuessInv[nStart]);
+    }
+    
+    PAR[PM_PAR_7] = Sm;
     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]);
@@ -1012,5 +1027,5 @@
     psVectorAppend (chi2, model->chisqNorm);
 
-    PAR[PM_PAR_7] = (nStart < N_INDEX_GUESS_INV - 1) ? 0.5*(indexGuessInv[nStart + 1] + indexGuessInv[nStart]) : indexGuessInv[nStart - 1];
+    PAR[PM_PAR_7] = So;
     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]);
@@ -1018,5 +1033,5 @@
     psVectorAppend (chi2, model->chisqNorm);
 
-    PAR[PM_PAR_7] = (nStart > 0) ? 0.5*(indexGuessInv[nStart - 1] + indexGuessInv[nStart]) : indexGuessInv[nStart + 1];
+    PAR[PM_PAR_7] = Sp;
     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]);
@@ -1036,6 +1051,8 @@
     float Smin = -0.5 * poly->coeff[1] / poly->coeff[2] / 100.0;
 
-    // constrain Smin to be in a valid range (1.0 - 0.1, corresponding to 0.5 (Gauss) to 5.0 (slightly peakier than Dev)
-    Smin = PS_MAX(PS_MIN(Smin, 1.0), 0.1);
+    // constrain Smin to be in a valid range: allow the fitted range to go a bit beyond the 3 trial points, but no further
+    float Smx = Sm - 0.25*(So - Sm);
+    float Spx = Sp + 0.25*(Sp - So);
+    Smin = PS_MAX(PS_MIN(Smin, Smx), Spx);
     PAR[PM_PAR_7] = Smin;
 
@@ -1249,7 +1266,7 @@
 
     // loop over Reff, keeping the ARatio and Theta constant
-    for (int j = -20; j < 21; j++) {
-
-	float dref = j * 0.02;
+    for (int j = -4; j <= 4; j++) {
+
+	float dref = j * 0.01;
 
 	psEllipseAxes guessAxes;
@@ -1306,5 +1323,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);
     }