Index: trunk/psphot/src/psphotGalaxyShape.c
===================================================================
--- trunk/psphot/src/psphotGalaxyShape.c	(revision 36633)
+++ trunk/psphot/src/psphotGalaxyShape.c	(revision 36718)
@@ -238,89 +238,89 @@
 bool psphotGalaxyShapeGrid (pmSource *source, pmSourceFitOptions *fitOptions, psphotGalaxyShapeOptions *opt, psImageMaskType maskVal, int psfSize) {
 
-    // we use the 0th model to define the initial guess shape
-    pmModel *model = source->modelFits->data[0];
-    if (!model) return false;
-
-    pmModelType modelType = model->type;
-
-    // we are using fitOptions->mode : be sure this makes sense
-    pmPCMdata *pcm = pmPCMinit (source, fitOptions, model, maskVal, psfSize);
-    if (!pcm) return false;
-
-    // we are fitting only PM_PAR_I0; the shape elements are generated from a grid
-    psF32 *PAR = pcm->modelConv->params->data.F32;
-
-    // double check that the guess is carried along...
-
-    // I have some source guess (e0, e1, e2)
-    psEllipseAxes guessAxes = pmPSF_ModelToAxes (PAR, modelType);
-
-    if (!source->galaxyFits) {
-	source->galaxyFits = pmSourceGalaxyFitsAlloc();
-    }
-
-    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) {
-  
-	    psEllipseAxes testAxes = guessAxes;
-	    testAxes.major = guessAxes.major * fRmajor;
-	    testAxes.minor = guessAxes.minor * fRminor;
-	    
-	    pmPSF_AxesToModel (PAR, testAxes, modelType);
-	    
-	    psphotGalaxyShapeSource (pcm, source, maskVal, psfSize, true);
-
-            int i = source->galaxyFits->chisq->n - 1;
-            float flux = source->galaxyFits->Flux->data.F32[i];
-            if (isfinite(flux)) {
-                float thisChisq = source->galaxyFits->chisq->data.F32[i];
-                if (isfinite(thisChisq) && isfinite(flux) && (!isfinite(chisqBest) || thisChisq < chisqBest)) {
-                    chisqBest = thisChisq;
-                    fRmajorBest = fRmajor;
-                    fRminorBest = fRminor;
+    for (int iModel = 0 ; iModel < source->modelFits->n; iModel++) {
+        pmModel *model = source->modelFits->data[iModel];
+        if (!model) return false;
+
+        pmModelType modelType = model->type;
+
+        // we are using fitOptions->mode : be sure this makes sense
+        pmPCMdata *pcm = pmPCMinit (source, fitOptions, model, maskVal, psfSize);
+        if (!pcm) return false;
+
+        // we are fitting only PM_PAR_I0; the shape elements are generated from a grid
+        psF32 *PAR = pcm->modelConv->params->data.F32;
+
+        // double check that the guess is carried along...
+
+        // I have some source guess (e0, e1, e2)
+        psEllipseAxes guessAxes = pmPSF_ModelToAxes (PAR, modelType);
+
+        if (!source->galaxyFits) {
+            source->galaxyFits = psArrayAllocEmpty(1);
+            // source->galaxyFits = pmSourceGalaxyFitsAlloc();
+        }
+        pmSourceGalaxyFits *galaxyFits = pmSourceGalaxyFitsAlloc();
+        psArrayAdd(source->galaxyFits, 1, galaxyFits);
+        psFree(galaxyFits);
+
+        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) {
+      
+                psEllipseAxes testAxes = guessAxes;
+                testAxes.major = guessAxes.major * fRmajor;
+                testAxes.minor = guessAxes.minor * fRminor;
+                
+                pmPSF_AxesToModel (PAR, testAxes, modelType);
+                
+                psphotGalaxyShapeSource (pcm, source, galaxyFits, maskVal, psfSize, true);
+
+                int i = galaxyFits->chisq->n - 1;
+                float flux = galaxyFits->Flux->data.F32[i];
+                if (isfinite(flux)) {
+                    float thisChisq = galaxyFits->chisq->data.F32[i];
+                    if (isfinite(thisChisq) && isfinite(flux) && (!isfinite(chisqBest) || thisChisq < chisqBest)) {
+                        chisqBest = thisChisq;
+                        fRmajorBest = fRmajor;
+                        fRminorBest = fRminor;
+                    }
                 }
+                // reset I0 to avoid potential problems on the next iteration
+                PAR[PM_PAR_I0] = 1.0;
             }
-            // reset I0 to avoid potential problems on the next iteration
-            PAR[PM_PAR_I0] = 1.0;
-	}
-    }
-
-#define SAVE_BEST_MODEL
-#ifdef SAVE_BEST_MODEL
-    // Save model with smallest chisq
-    if (isfinite(chisqBest)) {
-#else 
-    // Save model with nominal parameters
-    {
-        fRmajorBest = 1;
-        fRminorBest = 1;
+        }
+
+        if (isfinite(chisqBest)) {
+            // now save the best fitting model as the source's extended model ...
+            psEllipseAxes testAxes = guessAxes;
+
+            // ... unless this macro is defined
+#ifndef SAVE_NOMINAL_MODEL
+            testAxes.major = guessAxes.major * fRmajorBest;
+            testAxes.minor = guessAxes.minor * fRminorBest;
 #endif
-        // 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);
-
-        // Replace modelEXT with this model, if the model is good
-        if (isfinite(PAR[PM_PAR_I0])) {
-            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;
-
-            // cache the model flux
-            pmPCMCacheModel (source, maskVal, psfSize, fitOptions->nsigma);
-        }
-    }
-
-    psFree (pcm);
+            
+            pmPSF_AxesToModel (PAR, testAxes, modelType);
+                
+            psphotGalaxyShapeSource (pcm, source, galaxyFits, maskVal, psfSize, false);
+
+            // Replace modelEXT with the best model from the first of the model fits if one of them is good
+            if (isfinite(PAR[PM_PAR_I0]) && iModel == 0) {
+                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;
+
+                // cache the model flux
+                pmPCMCacheModel (source, maskVal, psfSize, fitOptions->nsigma);
+            }
+        }
+
+        psFree (pcm);
+    }
     return true;
 }
@@ -328,5 +328,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 saveResults) {
+bool psphotGalaxyShapeSource (pmPCMdata *pcm, pmSource *source, pmSourceGalaxyFits *galaxyFits, psImageMaskType maskVal, int psfSize, bool saveResults) {
 
     PS_ASSERT_PTR_NON_NULL(source, false);
@@ -384,8 +384,9 @@
 
     if (saveResults) {
-        psVectorAppend (source->galaxyFits->Flux, flux);
-        psVectorAppend (source->galaxyFits->dFlux, dflux);
-        psVectorAppend (source->galaxyFits->chisq, Chisq);
-        source->galaxyFits->nPix = nPix;
+        psVectorAppend (galaxyFits->Flux, flux);
+        psVectorAppend (galaxyFits->dFlux, dflux);
+        psVectorAppend (galaxyFits->chisq, Chisq);
+        galaxyFits->nPix = nPix;
+        galaxyFits->modelType = pcm->modelConv->type;
     }