Index: /branches/eam_branches/ipp-20130904/psphot/src/psphotExtendedSourceFits.c
===================================================================
--- /branches/eam_branches/ipp-20130904/psphot/src/psphotExtendedSourceFits.c	(revision 36309)
+++ /branches/eam_branches/ipp-20130904/psphot/src/psphotExtendedSourceFits.c	(revision 36310)
@@ -315,5 +315,5 @@
 
 // XXX TEST 
-	    if (false && !isInteractive) {
+	    if (!isInteractive) {
 		if (!psThreadJobAddPending(job)) {
 		    psError(PS_ERR_UNKNOWN, false, "Unable to guess model.");
Index: /branches/eam_branches/ipp-20130904/psphot/src/psphotSourceFits.c
===================================================================
--- /branches/eam_branches/ipp-20130904/psphot/src/psphotSourceFits.c	(revision 36309)
+++ /branches/eam_branches/ipp-20130904/psphot/src/psphotSourceFits.c	(revision 36310)
@@ -19,4 +19,6 @@
 static int NfitIterPCM = 0;
 static int NfitPixPCM = 0;
+
+bool psphotPCMfitCheckSize (pmPCMdata *pcm, pmSource *source, psImageMaskType maskVal, float psfSize);
 
 bool psphotFitInit (int nThreads) {
@@ -585,6 +587,6 @@
     }
 
-    float t1, t2, t4, t5;
-    t1 = t2 = t4 = t5 = 0.0;
+    float t1, t2, t3, t4, t5;
+    t1 = t2 = t3 = t4 = t5 = 0.0;
     if (TIMING) { psTimerStart ("psphotFitPCM"); }
 
@@ -642,5 +644,5 @@
     }
 
-    if (TIMING) { t4 = psTimerMark ("psphotFitPCM"); }
+    if (TIMING) { t3 = psTimerMark ("psphotFitPCM"); }
 
     // psTraceSetLevel("psLib.math.psMinimizeLMChi2", 5);
@@ -648,9 +650,12 @@
     NfitIterPCM += pcm->modelConv->nIter;
     NfitPixPCM += pcm->modelConv->nDOF;
+    if (TIMING) { t4 = psTimerMark ("psphotFitPCM"); }
+
+    psphotPCMfitCheckSize (pcm, source, maskVal, psfSize);
     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);
+    	fprintf (stderr, "psphotFitPCM : nIter: %2d, radius: %6.1f, npix: %5d of %5d, t1: %6.4f, t2: %6.4f, t3: %6.4f, t4: %6.4f, t5: %6.4f\n", model->nIter, model->fitRadius, model->nPix, nPixBig, t1, t2, t3, t4, t5);
     }
     if (EXTRA_VERBOSE && !TIMING) {
@@ -1246,3 +1251,95 @@
 }
 
- 
+// # define N_REFF_CHECK 11
+// float drefCheck[] = {-0.02, -0.04, -0.06, 0.0, 0.85, 0.90, 0.95, 1.00, 1.05, 1.10, 1.15, 1.20, 1.25};
+
+// we have an initial fit, check to see if the current size is besst
+bool psphotPCMfitCheckSize (pmPCMdata *pcm, pmSource *source, psImageMaskType maskVal, float psfSize) {
+
+    // PAR is already at my current best guess
+    psF32 *PAR = pcm->modelConv->params->data.F32;
+
+    // store best guess as a shape
+    psEllipseAxes centerAxes;
+    pmModelParamsToAxes (&centerAxes, PAR[PM_PAR_SXX], PAR[PM_PAR_SXY], PAR[PM_PAR_SYY], true);
+
+    float xMin = NAN;
+    float iMin = NAN;
+    float rMin = NAN;
+
+    // loop over Reff, keeping the ARatio and Theta constant
+    for (int j = -20; j < 21; j++) {
+
+	float dref = j * 0.02;
+
+	psEllipseAxes guessAxes;
+	guessAxes.major = centerAxes.major + dref;
+	guessAxes.minor = guessAxes.major * centerAxes.minor / centerAxes.major;
+	guessAxes.theta = centerAxes.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
+	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;
+	    rMin = dref;
+	}
+	// fprintf (stderr, "%d | %f %f %f | %f %f %f\n", j, dref, Io, Chisq, rMin, iMin, xMin);
+    }
+
+    psEllipseAxes guessAxes;
+    guessAxes.major = centerAxes.major + rMin;
+    guessAxes.minor = guessAxes.major * centerAxes.minor / centerAxes.major;
+    guessAxes.theta = centerAxes.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;
+
+    return true;
+}
+
Index: /branches/eam_branches/ipp-20130904/psphot/test/tap_psphot_galaxygrid.pro
===================================================================
--- /branches/eam_branches/ipp-20130904/psphot/test/tap_psphot_galaxygrid.pro	(revision 36309)
+++ /branches/eam_branches/ipp-20130904/psphot/test/tap_psphot_galaxygrid.pro	(revision 36310)
@@ -607,4 +607,5 @@
   $psphotConfig = $psphotConfig -Db PSPHOT:POISSON.ERRORS.PHOT.LMM F
   $psphotConfig = $psphotConfig -Db PSPHOT:EXTENDED_SOURCE_FITS_POISSON F
+   $psphotConfig = $psphotConfig -Di PSPHOT:EXT_FIT_ITER 40
 
   # ppImage / psphot on the output
@@ -614,4 +615,15 @@
   exec ppImage $ppImageConfig $psphotConfig -file $basename.fits $outname >>& $outname.log
   break -auto on
+end
+
+macro cmf.load.reset
+  $fields = X Y M T R r MT
+  foreach field $fields
+    foreach set in ot
+      delete -q $field\$set\_s
+    end
+  end
+
+  delete -q min_S Min_S
 end
 
@@ -1262,14 +1274,16 @@
   local radius Nrun 
 
+  if (1)
+    $Nrun = 0
+    foreach radius 1.0 1.5 2.0 2.5 3.0 4.0 6.0 8.0
+      # mkexp.devexp.single tests.20131120/testrad.ps1v1.exp.$Nrun EXP $radius 1.0
+      fitexp tests.20131120/testrad.ps1v1.exp.$Nrun tests.20131120/testrad.ps1v1.exp.$Nrun.t4s EXP_CONV
+      $Nrun ++
+    end
+  end
+
   $Nrun = 0
   foreach radius 1.0 1.5 2.0 2.5 3.0 4.0 6.0 8.0
-    mkexp.devexp.single tests.20131120/testrad.ps1v1.exp.$Nrun EXP $radius 1.0
-    fitexp tests.20131120/testrad.ps1v1.exp.$Nrun tests.20131120/testrad.ps1v1.exp.$Nrun.t4 EXP_CONV
-    $Nrun ++
-  end
-
-  $Nrun = 0
-  foreach radius 1.0 1.5 2.0 2.5 3.0 4.0 6.0 8.0
-    cmf.load.concat tests.20131120/testrad.ps1v1.exp.$Nrun.dat tests.20131120/testrad.ps1v1.exp.$Nrun.t4.cmf EXP
+    cmf.load.concat tests.20131120/testrad.ps1v1.exp.$Nrun.dat tests.20131120/testrad.ps1v1.exp.$Nrun.t4s.cmf EXP
     set dM_m = Mot_m - Min_m
     vstat -q dM_m
@@ -1284,2 +1298,44 @@
   grid.plots.devexp a
 end
+
+macro check.radius.loop
+  if ($0 != 2)
+    echo "USAGE: check.radius.loop (ext)"
+    break
+  end
+
+  # cmf.load.reset
+
+  local radius Nrun 
+
+  $Nrun = 0
+  foreach radius 1.0 1.5 2.0 2.5 3.0 4.0 6.0 8.0
+    cmf.load.concat tests.20131120/testrad.ps1v1.exp.$Nrun.dat tests.20131120/testrad.ps1v1.exp.$Nrun.$1.cmf EXP
+    set dM_m = Mot_m - Min_m
+    vstat -q dM_m
+    $mag_off = $MEDIAN
+    set dR_m = Rot_m - Rin_m
+    vstat -q dR_m
+    $rad_off = $MEDIAN
+    echo $radius $mag_off $rad_off
+    $Nrun ++
+  end
+
+  grid.plots.devexp a
+end
+
+macro check.fit
+  if ($0 != 3)
+    echo "USAGE: check.fit (dat) (cmf)"
+    break
+  end
+
+  cmf.load.concat $1 $2 EXP
+  set dM_m = Mot_m - Min_m
+  vstat -q dM_m
+  $mag_off = $MEDIAN
+  set dR_m = Rot_m - Rin_m
+  vstat -q dR_m
+  $rad_off = $MEDIAN
+  echo $mag_off $rad_off
+end