Index: branches/eam_branches/ipp-20130904/psphot/src/psphotGalaxyShape.c
===================================================================
--- branches/eam_branches/ipp-20130904/psphot/src/psphotGalaxyShape.c	(revision 36303)
+++ branches/eam_branches/ipp-20130904/psphot/src/psphotGalaxyShape.c	(revision 36304)
@@ -253,4 +253,7 @@
     }
 
+    float fRmajorBest = NAN;
+    float fRminorBest = NAN;
+    float chisqBest = NAN;
     for (float fRmajor = opt->fRmajorMin; fRmajor < opt->fRmajorMax + 0.5*opt->fRmajorDel; fRmajor += opt->fRmajorDel) {
 	for (float fRminor = opt->fRminorMin; fRminor < opt->fRminorMax + 0.5*opt->fRminorDel; fRminor += opt->fRminorDel) {
@@ -262,7 +265,38 @@
 	    pmPSF_AxesToModel (PAR, testAxes, modelType);
 	    
-	    // where do the results go??
-	    psphotGalaxyShapeSource (pcm, source, maskVal, psfSize);
+	    psphotGalaxyShapeSource (pcm, source, maskVal, psfSize, true);
+
+            int i = source->galaxyFits->chisq->n - 1;
+            float thisChisq = source->galaxyFits->chisq->data.F32[i];
+            if (isfinite(thisChisq) && (!isfinite(chisqBest) || thisChisq < chisqBest)) {
+                chisqBest = thisChisq;
+                fRmajorBest = fRmajor;
+                fRminorBest = fRminor;
+            }
 	}
+    }
+
+    if (isfinite(chisqBest)) {
+        // now save the best fitting model as the source's extended model
+        psEllipseAxes testAxes = guessAxes;
+        testAxes.major = guessAxes.major * fRmajorBest;
+        testAxes.minor = guessAxes.minor * fRminorBest;
+        
+        pmPSF_AxesToModel (PAR, testAxes, modelType);
+        
+        psphotGalaxyShapeSource (pcm, source, maskVal, psfSize, false);
+
+        psFree (source->modelEXT);
+
+        source->modelEXT = psMemIncrRefCounter (pcm->modelConv);
+        source->type = PM_SOURCE_TYPE_EXTENDED;
+        source->mode |= PM_SOURCE_MODE_EXTMODEL;
+        source->mode |= PM_SOURCE_MODE_NONLINEAR_FIT;
+
+        // adjust the window so the subtraction covers the faint wings
+        // psphotSetRadiusMoments(&fitRadius, &windowRadius, readout, source, markVal);
+
+        // cache the model flux
+        pmPCMCacheModel (source, maskVal, psfSize);
     }
 
@@ -273,5 +307,5 @@
 // fit the given model to the source and find chisq & normalization
 // XXX is this a single-component model? sersic with a supplied index, Reff, axis ratio, and theta?
-bool psphotGalaxyShapeSource (pmPCMdata *pcm, pmSource *source, psImageMaskType maskVal, int psfSize) {
+bool psphotGalaxyShapeSource (pmPCMdata *pcm, pmSource *source, psImageMaskType maskVal, int psfSize, bool saveResults) {
 
     PS_ASSERT_PTR_NON_NULL(source, false);
@@ -326,25 +360,11 @@
     float dflux = flux * (dIo / Io);
 
-// Set this to be the model of choice.
-// XXX: Maybe: figure out which is the best and save that one.
-{
-    psFree (source->modelEXT);
-
-    source->modelEXT = psMemIncrRefCounter (pcm->modelConv);
-    source->type = PM_SOURCE_TYPE_EXTENDED;
-    source->mode |= PM_SOURCE_MODE_EXTMODEL;
-    source->mode |= PM_SOURCE_MODE_NONLINEAR_FIT;
-
-    // adjust the window so the subtraction covers the faint wings
-    // psphotSetRadiusMoments(&fitRadius, &windowRadius, readout, source, markVal);
-
-    // cache the model flux
-    pmPCMCacheModel (source, maskVal, psfSize);
-}
-
-    psVectorAppend (source->galaxyFits->Flux, flux);
-    psVectorAppend (source->galaxyFits->dFlux, dflux);
-    psVectorAppend (source->galaxyFits->chisq, Chisq);
-    source->galaxyFits->nPix = nPix;
+
+    if (saveResults) {
+        psVectorAppend (source->galaxyFits->Flux, flux);
+        psVectorAppend (source->galaxyFits->dFlux, dflux);
+        psVectorAppend (source->galaxyFits->chisq, Chisq);
+        source->galaxyFits->nPix = nPix;
+    }
 
     return true;