Index: trunk/psphot/src/psphotSourceFits.c
===================================================================
--- trunk/psphot/src/psphotSourceFits.c	(revision 35769)
+++ trunk/psphot/src/psphotSourceFits.c	(revision 36086)
@@ -559,4 +559,10 @@
 
 # define TIMING 0
+# define EXTRA_VERBOSE 0
+
+bool psphotSersicModelGuessPCM (pmPCMdata *pcm, pmSource *source, psImageMaskType maskVal, float psfSize);
+bool psphotFitSersicShapeAndIndex (pmPCMdata *pcm, pmReadout *readout, pmSource *source, pmSourceFitOptions *fitOptions, psImageMaskType maskVal, psImageMaskType markVal, int psfSize);
+bool psphotFitSersicShapeAndIndexGrid (pmPCMdata *pcm, pmReadout *readout, pmSource *source, pmSourceFitOptions *fitOptions, psImageMaskType maskVal, psImageMaskType markVal, int psfSize);
+bool psphotFitSersicShapeAndIndexGridAuto (pmPCMdata *pcm, pmReadout *readout, pmSource *source, pmSourceFitOptions *fitOptions, psImageMaskType maskVal, psImageMaskType markVal, int psfSize);
 
 pmModel *psphotFitPCM (pmReadout *readout, pmSource *source, pmSourceFitOptions *fitOptions, pmModelType modelType, psImageMaskType maskVal, psImageMaskType markVal, int psfSize) {
@@ -580,5 +586,16 @@
 
     float t1, t2, t4, t5;
+    t1 = t2 = t4 = t5 = 0.0;
     if (TIMING) { psTimerStart ("psphotFitPCM"); }
+
+    if (modelType == pmModelClassGetType("PS_MODEL_SERSIC")) {
+	options.mode = PM_SOURCE_FIT_NO_INDEX; // XXX note that there may be a conflict with psphotExtendedSourceFits.c:133
+	options.mode = PM_SOURCE_FIT_EXT_AND_SKY;
+	// if we are ever (in a given psphot implementation) going to fit a parameter, we must set the options here to include
+	// that parameter (otherwise pmPCMupdate will fail to allocate the dmodelFlux image)
+	// thus, if the sersic analysis below uses an index fit, need to use this EXT_AND_SKY mode for init
+    } else {
+	options.mode = PM_SOURCE_FIT_EXT_AND_SKY;
+    }
 
     pmPCMdata *pcm = pmPCMinit (source, &options, model, maskVal, psfSize);
@@ -594,7 +611,18 @@
     if (modelType == pmModelClassGetType("PS_MODEL_SERSIC")) {
 	// use the source moments, etc to guess basic model parameters
-	if (!psphotSersicModelClassGuessPCM (pcm, source)) {
+	if (!psphotSersicModelGuessPCM (pcm, source, maskVal, psfSize)) {
 	    psFree (pcm);
 	    model->flags |= PM_MODEL_STATUS_BADARGS;
+	    return model;
+	}
+	if (TIMING) { t2 = psTimerMark ("psphotFitPCM"); }
+
+	// psphotFitSersicShapeAndIndex (pcm, readout, source, fitOptions, maskVal, markVal, psfSize);
+	options.mode = PM_SOURCE_FIT_NO_INDEX;
+	if (!psphotFitSersicShapeAndIndexGridAuto (pcm, readout, source, &options, maskVal, markVal, psfSize)) {
+	    psFree (pcm);
+	    model->flags |= PM_MODEL_STATUS_BADARGS;
+	    psError(PS_ERR_UNKNOWN, true, "Failed to find a index & shape");
+	    psErrorClear (); // clear the polynomial error
 	    return model;
 	}
@@ -608,13 +636,4 @@
     }
 
-    if (TIMING) { t2 = psTimerMark ("psphotFitPCM"); }
-
-    if (modelType == pmModelClassGetType("PS_MODEL_SERSIC")) {
-	options.mode = PM_SOURCE_FIT_NO_INDEX;
-    } else {
-	options.mode = PM_SOURCE_FIT_EXT;
-    }
-    // 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"); }
 
@@ -629,6 +648,10 @@
     	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);
     }
-
-    // psTraceSetLevel("psLib.math.psMinimizeLMChi2", 0);
+    if (EXTRA_VERBOSE && !TIMING) {
+	int nPixBig = source->pixels->numCols * source->pixels->numRows;
+	float *PAR = model->params->data.F32;
+	fprintf (stderr, "source %d : %f - %f %f - %f %f %f - %f | nIter: %2d, radius: %6.1f, npix: %5d of %5d, chisq %f\n", source->id, PAR[7], PAR[2], PAR[3], PAR[4], PAR[5], PAR[6], PAR[1], model->nIter, model->fitRadius, model->nPix, nPixBig, model->chisqNorm);
+    }
+
     psFree (pcm);
 
@@ -748,2 +771,472 @@
     return true;
 }
+
+// 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};
+# define N_REFF_GUESS 5
+
+// A sersic model is very sensitive to the index.  attempt to find the index first by grid search in just the index
+// for a sersic model, attempt to fit just the index and normalization with a modest number of iterations
+bool psphotSersicModelGuessPCM (pmPCMdata *pcm, pmSource *source, psImageMaskType maskVal, float psfSize) {
+
+    // we get a reasonable guess from:
+    // * Reff = Kron R1
+    // * Rmajor / Rminor & Theta from moments
+    // * Io from total Kron flux
+
+    // the guesses are used to fill in PAR:
+    psF32 *PAR = pcm->modelConv->params->data.F32;
+
+    // convert the moments to Major,Minor,Theta
+    psEllipseMoments moments;
+
+    if (!isfinite(source->moments->Mrf)) return false;
+    if (!isfinite(source->moments->Mxx)) return false;
+    if (!isfinite(source->moments->Mxy)) return false;
+    if (!isfinite(source->moments->Myy)) return false;
+
+    moments.x2 = source->moments->Mxx;
+    moments.y2 = source->moments->Myy;
+    moments.xy = source->moments->Mxy;
+    
+    // limit axis ratio < 20.0
+    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)) {
+      return false;
+    }
+
+    // sky is zero (no longer fitted, but not yet deprecated)
+    PAR[PM_PAR_SKY]  = 0.0;
+
+    // for the index loop, use Io = 1.0, use fitted values to determine Io
+    PAR[PM_PAR_I0] = 1.0;
+
+    float xMin = NAN;
+    float iMin = NAN;
+    float sMin = NAN;
+    float rMin = NAN;
+
+    // loop over index and Reff, keeping the ARatio and Theta constant?
+    // 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];
+
+	    // generated the modelFlux
+	    pmPCMMakeModel (source, pcm->modelConv, maskVal, psfSize);
+	
+	    float YY = 0.0;
+	    float YM = 0.0;
+	    float MM = 0.0;
+	    bool usePoisson = false;
+
+	    for (int iy = 0; iy < source->pixels->numRows; iy++) {
+		for (int ix = 0; ix < source->pixels->numCols; ix++) {
+		    // skip masked points
+		    if (source->maskObj->data.PS_TYPE_IMAGE_MASK_DATA[iy][ix]) {
+			continue;
+		    }
+		    // skip zero-variance points
+		    if (source->variance->data.F32[iy][ix] == 0) {
+			continue;
+		    }
+		    // skip nan value points
+		    if (!isfinite(source->pixels->data.F32[iy][ix])) {
+			continue;
+		    }
+
+		    float fy = source->pixels->data.F32[iy][ix];
+		    float fm = source->modelFlux->data.F32[iy][ix];
+		    float wt = (usePoisson) ? 1.0 / source->variance->data.F32[iy][ix] : 1.0;
+
+		    YY += PS_SQR(fy) * wt;
+		    YM += fm * fy * wt;
+		    MM += PS_SQR(fm) * wt;
+		}
+	    }
+
+	    float Io = YM / MM;
+	    float Chisq = YY - 2 * Io * YM + Io * Io * MM;
+	    if (isnan(xMin) || (Chisq < xMin)) {
+		xMin = Chisq;
+		iMin = Io;
+		sMin = indexGuessInv[i];
+		rMin = reffGuess[j];
+	    }
+	    // 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");
+	}
+    }
+
+    {
+	psEllipseAxes guessAxes;
+	guessAxes.major = rMin * 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);
+    }
+
+    PAR[PM_PAR_I0] = iMin;
+    PAR[PM_PAR_7] = sMin;
+
+    return true;
+}
+
+// we have a set of guess parameters, do a small number of iterations fitting only SHAPE then only INDEX 
+bool psphotFitSersicShapeAndIndex (pmPCMdata *pcm, pmReadout *readout, pmSource *source, pmSourceFitOptions *fitOptions, psImageMaskType maskVal, psImageMaskType markVal, int psfSize) {
+
+    pmModel *model = pcm->modelConv;
+
+    assert (model->type == pmModelClassGetType("PS_MODEL_SERSIC"));
+
+    pmSourceFitOptions options = *fitOptions;
+    
+    for (int i = 0; i < 3; i++) {
+      // fit EXT (not PSF) model (set/unset the pixel mask)
+      options.mode = PM_SOURCE_FIT_SHAPE;
+      options.nIter = 2;
+
+      // update the pcm elements if we have changed the circumstance (here, options.mode)
+      pmPCMupdate(pcm, source, &options, model);
+      
+      pmSourceFitPCM (pcm, source, &options, maskVal, markVal, psfSize);
+      if (EXTRA_VERBOSE) {
+	float *PAR = model->params->data.F32;
+	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]);
+      }
+      
+      // fit EXT (not PSF) model (set/unset the pixel mask)
+      options.mode = PM_SOURCE_FIT_INDEX;
+      // options.mode = PM_SOURCE_FIT_EXT_AND_SKY;
+      options.nIter = 30;
+      
+      // update the pcm elements if we have changed the circumstance (here, options.mode)
+      pmPCMupdate(pcm, source, &options, model);
+      
+      pmSourceFitPCM (pcm, source, &options, maskVal, markVal, psfSize);
+      if (EXTRA_VERBOSE) {
+	float *PAR = model->params->data.F32;
+	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]);
+      }
+    }
+
+    // update the pcm elements if we have changed the circumstance (here, options.mode)
+    pmPCMupdate(pcm, source, fitOptions, model);
+
+    return true;
+}
+
+// we have a set of guess parameters, do a small number of iterations fitting only SHAPE then only INDEX 
+bool psphotFitSersicShapeAndIndexGridAuto (pmPCMdata *pcm, pmReadout *readout, pmSource *source, pmSourceFitOptions *fitOptions, psImageMaskType maskVal, psImageMaskType markVal, int psfSize) {
+
+    pmModel *model = pcm->modelConv;
+
+    assert (model->type == pmModelClassGetType("PS_MODEL_SERSIC"));
+
+    pmSourceFitOptions options = *fitOptions;
+    
+    psF32 *PAR = pcm->modelConv->params->data.F32;
+
+    options.mode = PM_SOURCE_FIT_SHAPE;
+    options.nIter = 7;
+    
+    // update the pcm elements if we have changed the circumstance (here, options.mode)
+    pmPCMupdate(pcm, source, &options, model);
+    
+    // we have been provided a guess at the index (P[7]) from the list of indexGuessInv
+
+    // find the matching indexGuessInv
+    int nStart = -1;
+    for (int i = 0; i < N_INDEX_GUESS_INV; i++) {
+	if (fabs(PAR[PM_PAR_7] - indexGuessInv[i]) < 0.01) {
+	    nStart = i;
+	    break;
+	}
+    }
+    if (nStart == -1) {
+	fprintf (stderr, "WARNING: could not find start guess %f\n", PAR[PM_PAR_7]);
+	return false;
+    }
+
+    psVector *chi2 = psVectorAllocEmpty (16, PS_TYPE_F32);
+    psVector *Sidx = psVectorAllocEmpty (16, PS_TYPE_F32);
+
+    PAR[PM_PAR_7] = indexGuessInv[nStart];
+    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);
+
+    PAR[PM_PAR_7] = (nStart < N_INDEX_GUESS_INV - 1) ? 0.5*(indexGuessInv[nStart + 1] + indexGuessInv[nStart]) : 0.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);
+
+    PAR[PM_PAR_7] = (nStart > 0) ? 0.5*(indexGuessInv[nStart - 1] + indexGuessInv[nStart]) : 0.55;
+    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);
+    if (!psVectorFitPolynomial1D (poly, NULL, 0, chi2, NULL, Sidx)) {
+        psError(PS_ERR_UNKNOWN, true, "Failed to find a good chisq parabola");
+	psFree (chi2);
+	psFree (Sidx);
+	psFree (poly);
+	return false;
+    }
+
+    // where is the minimum of this polynomial fit?
+    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);
+    PAR[PM_PAR_7] = Smin;
+
+    // return to the original fitting mode (fitOptions)
+    pmPCMupdate(pcm, source, fitOptions, model);
+
+    psFree (chi2);
+    psFree (Sidx);
+    psFree (poly);
+
+    return true;
+}
+
+ 
+// we have a set of guess parameters, do a small number of iterations fitting only SHAPE then only INDEX 
+bool psphotFitSersicShapeAndIndexGridAutoScaled (pmPCMdata *pcm, pmReadout *readout, pmSource *source, pmSourceFitOptions *fitOptions, psImageMaskType maskVal, psImageMaskType markVal, int psfSize) {
+
+    pmModel *model = pcm->modelConv;
+
+    assert (model->type == pmModelClassGetType("PS_MODEL_SERSIC"));
+
+    pmSourceFitOptions options = *fitOptions;
+    
+    psF32 *PAR = pcm->modelConv->params->data.F32;
+
+    options.mode = PM_SOURCE_FIT_SHAPE;
+    options.nIter = 5;
+    
+    // update the pcm elements if we have changed the circumstance (here, options.mode)
+    pmPCMupdate(pcm, source, &options, model);
+    
+    float parStart[8];
+    for (int i = 0; i < 8; i++) parStart[i] = PAR[i];
+
+    // we start with a guess at the index (P[7])
+
+    // get chisq for P[7], P[7]*1.1, P[7]*1.25 (or *0.75 depending on the result of 1.1)
+
+    psVector *chi2 = psVectorAllocEmpty (16, PS_TYPE_F32);
+    psVector *Sidx = psVectorAllocEmpty (16, PS_TYPE_F32);
+
+    PAR[PM_PAR_7] = parStart[7];
+    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);
+
+    float fI = 1.1;
+    PAR[PM_PAR_7] = parStart[7]*fI;
+    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[0]) {
+      fI = 1.3;
+    } else {
+      fI = 1.0 / 1.3;
+    }
+    
+    PAR[PM_PAR_7] = parStart[7]*fI;
+    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);
+
+    // can we fit the 3 pts with a parabola?
+    int nTry = 0;
+    psPolynomial1D *poly = psPolynomial1DAlloc (PS_POLYNOMIAL_ORD, 2);
+    while (!psVectorFitPolynomial1D (poly, NULL, 0, chi2, NULL, Sidx)) {
+      psErrorClear (); // clear the polynomial error
+      if (nTry > 4) {
+        psError(PS_ERR_UNKNOWN, true, "Failed to find a good chisq parabola");
+	psFree (chi2);
+	psFree (Sidx);
+	psFree (poly);
+	return false;
+      }
+      fI = (fI < 1.0) ? fI / 1.3 : fI * 1.3;
+      PAR[PM_PAR_7] = parStart[7]*fI;
+      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);
+      nTry ++;
+    }
+
+    // where is the minimum of this polynomial fit?
+    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);
+    PAR[PM_PAR_7] = Smin;
+
+    // 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]);
+    
+    //// for (int i = 0; i < 8; i++) PAR[i] = parStart[i];
+    //// 
+    //// for (float fI = 0.0; fI < 0.15; fI += 0.01) {
+    ////   PAR[PM_PAR_7] = parStart[7] - fI;
+    //// 
+    ////   // fit EXT (not PSF) model (set/unset the pixel mask)
+    ////   
+    ////   pmSourceFitPCM (pcm, source, &options, maskVal, markVal, psfSize);
+    ////   if (TIMING) {
+    //// 	float *PAR = model->params->data.F32;
+    //// 	fprintf (stderr, "%d %f : %f - %f %f %f - %f\n", model->nIter, model->chisqNorm, PAR[7], PAR[4], PAR[5], PAR[6], PAR[1]);
+    ////   }
+    //// }
+
+    // return to the original fitting mode (fitOptions)
+    pmPCMupdate(pcm, source, fitOptions, model);
+
+    psFree (chi2);
+    psFree (Sidx);
+    psFree (poly);
+
+    return true;
+}
+
+ 
+// we have a set of guess parameters, do a small number of iterations fitting only SHAPE then only INDEX 
+bool psphotFitSersicShapeAndIndexGrid (pmPCMdata *pcm, pmReadout *readout, pmSource *source, pmSourceFitOptions *fitOptions, psImageMaskType maskVal, psImageMaskType markVal, int psfSize) {
+
+    pmModel *model = pcm->modelConv;
+
+    assert (model->type == pmModelClassGetType("PS_MODEL_SERSIC"));
+
+    pmSourceFitOptions options = *fitOptions;
+    
+    psF32 *PAR = pcm->modelConv->params->data.F32;
+
+    options.mode = PM_SOURCE_FIT_SHAPE;
+    options.nIter = 10;
+    
+    // update the pcm elements if we have changed the circumstance (here, options.mode)
+    pmPCMupdate(pcm, source, &options, model);
+    
+    psVector *chi2 = psVectorAllocEmpty (16, PS_TYPE_F32);
+    psVector *Sidx = psVectorAllocEmpty (16, PS_TYPE_F32);
+
+    float par7[] = {0.100, 0.125, 0.150, 0.175, 0.200, 0.225, 0.250};
+    for (int i = 0; i < 7; i++) {
+      PAR[PM_PAR_7] = par7[i];
+      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, PAR[PM_PAR_7]);
+      psVectorAppend (chi2, model->chisqNorm);
+    }
+
+    psPolynomial1D *poly = psPolynomial1DAlloc (PS_POLYNOMIAL_ORD, 2);
+    if (!psVectorFitPolynomial1D (poly, NULL, 0, chi2, NULL, Sidx)) {
+      psError(PS_ERR_UNKNOWN, true, "Failed to find a good chisq parabola");
+      psFree (chi2);
+      psFree (Sidx);
+      psFree (poly);
+      return false;
+    }
+
+    // where is the minimum of this polynomial fit?
+    fprintf (stderr, "fit1d: %f + %f x + %f x^2\n", poly->coeff[0], poly->coeff[1], poly->coeff[2]);
+    float Smin = -0.5 * poly->coeff[1] / poly->coeff[2];
+
+    // 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);
+    PAR[PM_PAR_7] = Smin;
+    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]);
+    
+    //// for (int i = 0; i < 8; i++) PAR[i] = parStart[i];
+    //// 
+    //// for (float fI = 0.0; fI < 0.15; fI += 0.01) {
+    ////   PAR[PM_PAR_7] = parStart[7] - fI;
+    //// 
+    ////   // fit EXT (not PSF) model (set/unset the pixel mask)
+    ////   
+    ////   pmSourceFitPCM (pcm, source, &options, maskVal, markVal, psfSize);
+    ////   if (TIMING) {
+    //// 	float *PAR = model->params->data.F32;
+    //// 	fprintf (stderr, "%d %f : %f - %f %f %f - %f\n", model->nIter, model->chisqNorm, PAR[7], PAR[4], PAR[5], PAR[6], PAR[1]);
+    ////   }
+    //// }
+
+    // return to the original fitting mode (fitOptions)
+    pmPCMupdate(pcm, source, fitOptions, model);
+
+    psFree (chi2);
+    psFree (Sidx);
+    psFree (poly);
+
+    return true;
+}
+
+