Index: /branches/eam_branches/ipp-20130711/psphot/src/psphotSourceFits.c
===================================================================
--- /branches/eam_branches/ipp-20130711/psphot/src/psphotSourceFits.c	(revision 36024)
+++ /branches/eam_branches/ipp-20130711/psphot/src/psphotSourceFits.c	(revision 36025)
@@ -558,5 +558,7 @@
 }
 
-# define TIMING 0
+# define TIMING 1
+
+bool psphotSersicModelGuessPCM (pmPCMdata *pcm, pmSource *source, psImageMaskType maskVal, float psfSize);
 
 pmModel *psphotFitPCM (pmReadout *readout, pmSource *source, pmSourceFitOptions *fitOptions, pmModelType modelType, psImageMaskType maskVal, psImageMaskType markVal, int psfSize) {
@@ -580,5 +582,12 @@
 
     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;
+    } else {
+	options.mode = PM_SOURCE_FIT_EXT_AND_SKY; // XXX note that there may be a conflict with psphotExtendedSourceFits.c:133
+    }
 
     pmPCMdata *pcm = pmPCMinit (source, &options, model, maskVal, psfSize);
@@ -594,5 +603,5 @@
     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;
@@ -613,12 +622,7 @@
     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_AND_SKY; // XXX note that there may be a conflict with psphotExtendedSourceFits.c:133
-    }
     // 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"); }
+    // pmPCMupdate(pcm, source, &options, model);
+    // if (TIMING) { t4 = psTimerMark ("psphotFitPCM"); }
 
     // psTraceSetLevel("psLib.math.psMinimizeLMChi2", 5);
@@ -751,2 +755,115 @@
     return true;
 }
+
+// float indexGuessInv[] = {0.5, 0.33, 0.25, 0.167, 0.125, 0.083};
+float indexGuessInv[] = {0.5, 0.25, 0.167, 0.125};
+# define N_INDEX_GUESS_INV 4
+
+// 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);
+
+    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;
+
+    // loop over index and Reff, keeping the ARatio and Theta constant?
+    // loop over index guesses and find the best fit
+    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 ((i == 0) || (Chisq < xMin)) {
+	    xMin = Chisq;
+	    iMin = Io;
+	    sMin = indexGuessInv[i];
+	}
+	fprintf (stderr, "%d | %f %f %f | %f %f %f", i, indexGuessInv[i], Io, Chisq, sMin, iMin, xMin);
+	fprintf (stderr, "\n");
+    }
+
+    PAR[PM_PAR_I0] = iMin;
+    PAR[PM_PAR_7] = sMin;
+
+    return true;
+}