Index: trunk/psphot/src/psphot.h
===================================================================
--- trunk/psphot/src/psphot.h	(revision 36701)
+++ trunk/psphot/src/psphot.h	(revision 36718)
@@ -537,5 +537,5 @@
 bool psphotGalaxyShape_Threaded (psThreadJob *job);
 bool psphotGalaxyShapeGrid (pmSource *source, pmSourceFitOptions *fitOptions, psphotGalaxyShapeOptions *opt, psImageMaskType maskVal, int psfSize);
-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);
 psphotGalaxyShapeOptions *psphotGalaxyShapeOptionsAlloc();
 
Index: trunk/psphot/src/psphotFullForceSummary.c
===================================================================
--- trunk/psphot/src/psphotFullForceSummary.c	(revision 36701)
+++ trunk/psphot/src/psphotFullForceSummary.c	(revision 36718)
@@ -187,12 +187,25 @@
 
     // Convert to an array indexed by ID
-    psArray *sortedTable = psArrayAlloc(4*inTable->n);
+#define MAX_MODEL_TYPE 10
+    psArray *sortedTables = psArrayAlloc(MAX_MODEL_TYPE+1);
     for (int i=0; i<inTable->n; i++) {
         psMetadata *row = inTable->data[i];
         psS32 ID = psMetadataLookupS32(&status, row, "ID");
+        psS32 modelType = psMetadataLookupS32(&status, row, "MODEL_TYPE");
+        // XXX: need to use the lookup table functions to be ready for changes in the model type numbers
+        if (modelType+1 >= MAX_MODEL_TYPE) {
+            psError(PSPHOT_ERR_CONFIG, false, "found modelType %d max allowed is %d", modelType, MAX_MODEL_TYPE);
+            return false;
+        }
+        psArray *sortedTable = sortedTables->data[modelType+1];
+        if (!sortedTable) {
+            sortedTable = psArrayAlloc(4*inTable->n);
+            sortedTables->data[modelType+1] = sortedTable;
+            // dont' free sortedTable the array of tables gets our reference
+        }
         if (ID >= sortedTable->n) {
             sortedTable = psArrayRealloc(sortedTable, 2*ID);
             // Why doesn't psArrayRealloc do this?????
-            sortedTable->n = sortedTable->nalloc;
+                sortedTable->n = sortedTable->nalloc;
         }
         if (sortedTable->data[ID]) {
@@ -202,7 +215,7 @@
         sortedTable->data[ID] = psMemIncrRefCounter(row);
     }
-    psMetadataAddArray(config->arguments, PS_LIST_TAIL, "CFF_TABLE", PS_META_REPLACE, "cff table", sortedTable);
+    psMetadataAddArray(config->arguments, PS_LIST_TAIL, "CFF_TABLES", PS_META_REPLACE, "cff tables", sortedTables);
     psFree(inTable);
-    psFree(sortedTable);
+    psFree(sortedTables);
 
     return true;
Index: trunk/psphot/src/psphotFullForceSummaryReadout.c
===================================================================
--- trunk/psphot/src/psphotFullForceSummaryReadout.c	(revision 36701)
+++ trunk/psphot/src/psphotFullForceSummaryReadout.c	(revision 36718)
@@ -14,5 +14,5 @@
     psArray *zeroPt;
     psArray *exptime;
-    psArray *cffTable;
+    psArray *cffTables;
 } galaxyShapeOptions;
 
@@ -71,6 +71,6 @@
                 return false;
             }
-            options.cffTable = psMetadataLookupPtr(NULL, config->arguments, "CFF_TABLE");
-            if (!options.cffTable) {
+            options.cffTables = psMetadataLookupPtr(NULL, config->arguments, "CFF_TABLES");
+            if (!options.cffTables) {
                 psError (PS_ERR_UNKNOWN, true, "Cannot find cff table in arguments.");
                 return false;
@@ -178,9 +178,19 @@
     pmSource *outSrc = NULL;
 
-    psVector *sumWeightedFlux = NULL;
-    psVector *sumInvSig2 = NULL;
-    psVector *numerator   = NULL;
-    psF32   totalNPix = 0;
-    long    vectorLength = 0;
+
+    pmSource *source0 = obj->sources->data[0];
+    int ID = source0->seq;
+
+    // Find the parameters from the cff tables
+    int nModels = 0;
+    // skip table 0 which is modelType -1 not extended
+    for (int i=1; i<options->cffTables->n; i++) {
+        psArray *table = options->cffTables->data[i];
+        if (!table) continue;
+        if (ID >= table->n) continue;
+        if (table->data[ID]) {
+            nModels++;
+        }
+    }
 
     // Loop over sources for this object. Start output source for first input that has
@@ -210,91 +220,123 @@
                 outSrc->extpars =  psMemIncrRefCounter(source->extpars);
             }
-            if (source->modelEXT) {
-                outSrc->modelEXT =  psMemIncrRefCounter(source->modelEXT);
-            }
-            if (source->modelFits) {
-                outSrc->modelFits =  psMemIncrRefCounter(source->modelFits);
-            }
-        }
-
-        // The only parameters that we are summarizing currently are galaxy shapes. 
-        // Continue for inputs that do not have modelEXT which contains the model paramters
-        if (!source->modelEXT) continue;
-
-        if (source->galaxyFits && isfinite(source->galaxyFits->nPix) && source->galaxyFits->chisq->n) {
-            if (numerator == NULL) {
-                // first source with galaxyFits allocate accumulators
-                vectorLength = source->galaxyFits->chisq->n;
-                sumWeightedFlux = psVectorAlloc(vectorLength, PS_TYPE_F32);
-                psVectorInit(sumWeightedFlux, 0.0);
-                sumInvSig2 = psVectorAlloc(vectorLength, PS_TYPE_F32);
-                psVectorInit(sumInvSig2, 0.0);
-                numerator   = psVectorAlloc(vectorLength, PS_TYPE_F32);
-                psVectorInit(numerator, 0.0);
-            }
-
-            // Die if the lengths of the vectors is not the same in all sources
-            psAssert(vectorLength == options->numTrials, "length of chisq vector %ld does not match options %d",
-                vectorLength, options->numTrials);
-            psAssert(source->galaxyFits->chisq->n == vectorLength, "length of chisq vectors do not match %ld %ld",
-                    source->galaxyFits->chisq->n, vectorLength);
-
-            psF32 scaleFactor = fluxScaleFactor->data.F32[source->imageID];
-
-            totalNPix += source->galaxyFits->nPix;
-
-            for (int k = 0; k < vectorLength; k++) {
-
-                psF32 chisq = source->galaxyFits->chisq->data.F32[k];
-                psF32 flux  = source->galaxyFits->Flux->data.F32[k]  * scaleFactor;
-                psF32 dFlux = source->galaxyFits->dFlux->data.F32[k] * scaleFactor;
-
-                numerator->data.F32[k] += chisq * source->galaxyFits->nPix;
-
-                psF32 invSig2 = 1.0 / (dFlux * dFlux);
-                sumInvSig2->data.F32[k] += invSig2;
-
-                sumWeightedFlux->data.F32[k] += flux * invSig2;
-            }
-        }
-    }
-
-    if (outSrc) {
+            outSrc->modelEXT = psFree(outSrc->modelEXT);
+        }
+        if (nModels && !outSrc->modelFits && source->modelFits && source->modelFits->n == nModels) {
+            outSrc->modelFits =  psMemIncrRefCounter(source->modelFits);
+        }
+    }
+
+    if (!outSrc) {
+        // no good measurements for this source (how?)
+        return NULL;
+    }
+
+    // now loop over the model fits and summarize the galaxy shape measurements for each
+    for (int iModel=0; outSrc->modelFits && iModel < outSrc->modelFits->n; iModel++) {
+        pmModel *model = outSrc->modelFits->data[iModel];
+        int modelType = model->type;
+        psAssert(modelType >= 0 && modelType+1 < options->cffTables->n, "model type out of range");
+
+        psArray *cffTable = options->cffTables->data[modelType+1];
+        psAssert(cffTable, "missing cff table?");
+
+        psMetadata *cffRow = cffTable->data[outSrc->seq];
+        if (!cffRow) {
+            psError (PS_ERR_PROGRAMMING, true, "NO cff data for object %d model %d", outSrc->seq, modelType);
+            psFree(outSrc);
+            return NULL;
+        }
+        bool mdok;
+        psF32 R_MAJOR = psMetadataLookupF32(&mdok, cffRow, "R_MAJOR");
+        if (!mdok) {
+            psError (PS_ERR_PROGRAMMING, true, "can't find R_MAJOR for object %d type %d", outSrc->seq, modelType);
+            psFree(outSrc);
+            return NULL;
+        }
+        psF32 R_MINOR = psMetadataLookupF32(&mdok, cffRow, "R_MINOR");
+        if (!mdok) {
+            psError (PS_ERR_PROGRAMMING, true, "can't find R_MINOR for object %d type %d", outSrc->seq, modelType);
+            psFree(outSrc);
+            return NULL;
+        }
+
+        psVector *sumWeightedFlux = NULL;
+        psVector *sumInvSig2 = NULL;
+        psVector *numerator   = NULL;
+        psF32   totalNPix = 0;
+        long    vectorLength = 0;
+        for (int i=0; i < obj->sources->n; i++) {
+            pmSource *source = obj->sources->data[i];
+
+            if (source->pixWeightNotPoor < .9) continue;
+
+            // The only parameters that we are summarizing currently are galaxy shapes. 
+            if (!source->modelFits) continue;
+            if (!source->galaxyFits) continue;
+
+            // XXX: put this into a function
+            pmSourceGalaxyFits *galaxyFits = NULL;
+            for (int j=0; j<source->galaxyFits->n; j++) {
+                galaxyFits = source->galaxyFits->data[j];
+                if (galaxyFits && galaxyFits->modelType == modelType) break;
+                galaxyFits = NULL;
+            }
+            // end of proposed function
+
+            if (galaxyFits && isfinite(galaxyFits->nPix) && galaxyFits->chisq->n) {
+                if (numerator == NULL) {
+                    // first source with galaxyFits allocate accumulators
+                    vectorLength = galaxyFits->chisq->n;
+                    sumWeightedFlux = psVectorAlloc(vectorLength, PS_TYPE_F32);
+                    psVectorInit(sumWeightedFlux, 0.0);
+                    sumInvSig2 = psVectorAlloc(vectorLength, PS_TYPE_F32);
+                    psVectorInit(sumInvSig2, 0.0);
+                    numerator   = psVectorAlloc(vectorLength, PS_TYPE_F32);
+                    psVectorInit(numerator, 0.0);
+                }
+
+                // Die if the lengths of the vectors is not the same in all sources
+                psAssert(vectorLength == options->numTrials, "length of chisq vector %ld does not match options %d",
+                    vectorLength, options->numTrials);
+                psAssert(galaxyFits->chisq->n == vectorLength, "length of chisq vectors do not match %ld %ld",
+                             galaxyFits->chisq->n, vectorLength);
+
+                psF32 scaleFactor = fluxScaleFactor->data.F32[source->imageID];
+
+                totalNPix += galaxyFits->nPix;
+
+                for (int k = 0; k < vectorLength; k++) {
+                    psF32 chisq = galaxyFits->chisq->data.F32[k];
+                    psF32 flux  = galaxyFits->Flux->data.F32[k]  * scaleFactor;
+                    psF32 dFlux = galaxyFits->dFlux->data.F32[k] * scaleFactor;
+
+                    numerator->data.F32[k] += chisq * galaxyFits->nPix;
+
+                    psF32 invSig2 = 1.0 / (dFlux * dFlux);
+                    sumInvSig2->data.F32[k] += invSig2;
+
+                    sumWeightedFlux->data.F32[k] += flux * invSig2;
+                }
+            }
+        }
+
         if (vectorLength) {
             // allocate galaxyFits for the output source
-            outSrc->galaxyFits = pmSourceGalaxyFitsAlloc();
-            outSrc->galaxyFits->nPix  = totalNPix;
+            if (!outSrc->galaxyFits) {
+                outSrc->galaxyFits = psArrayAllocEmpty(1);
+            }
+            pmSourceGalaxyFits *galaxyFits = pmSourceGalaxyFitsAlloc();
+            psArrayAdd(outSrc->galaxyFits, 1, galaxyFits);
+            psFree(galaxyFits);
+
+            galaxyFits->nPix  = totalNPix;
+            galaxyFits->modelType = modelType;
             psVector *fluxVec = 
-                outSrc->galaxyFits->Flux  = psVectorRecycle(outSrc->galaxyFits->Flux,  vectorLength, PS_TYPE_F32);
+                galaxyFits->Flux  = psVectorRecycle(galaxyFits->Flux,  vectorLength, PS_TYPE_F32);
             psVector *dFluxVec = 
-                outSrc->galaxyFits->dFlux = psVectorRecycle(outSrc->galaxyFits->dFlux, vectorLength, PS_TYPE_F32);
+                galaxyFits->dFlux = psVectorRecycle(galaxyFits->dFlux, vectorLength, PS_TYPE_F32);
             psVector *chisqVec = 
-                outSrc->galaxyFits->chisq = psVectorRecycle(outSrc->galaxyFits->chisq, vectorLength, PS_TYPE_F32);
-
-            // Find the nominal radii paramters in the cff table
-            if (outSrc->seq >= options->cffTable->n) {
-                psError(PS_ERR_PROGRAMMING, true, "object sequence number %d is larger than cff table length %ld",
-                    outSrc->seq, options->cffTable->n);
-                return NULL;
-            }
-            psMetadata *cffRow = options->cffTable->data[outSrc->seq];
-            if (!cffRow) {
-                psError (PS_ERR_PROGRAMMING, true, "NO cff data for object %d", outSrc->seq);
-                psFree(outSrc);
-                return NULL;
-            }
-            bool mdok;
-            psF32 R_MAJOR = psMetadataLookupF32(&mdok, cffRow, "R_MAJOR");
-            if (!mdok) {
-                psError (PS_ERR_PROGRAMMING, true, "can't find R_MAJOR for object %d", outSrc->seq);
-                psFree(outSrc);
-                return NULL;
-            }
-            psF32 R_MINOR = psMetadataLookupF32(&mdok, cffRow, "R_MINOR");
-            if (!mdok) {
-                psError (PS_ERR_PROGRAMMING, true, "can't find R_MINOR for object %d", outSrc->seq);
-                psFree(outSrc);
-                return NULL;
-            }
+                galaxyFits->chisq = psVectorRecycle(galaxyFits->chisq, vectorLength, PS_TYPE_F32);
+
 
             // fill the summary galaxyFits vectors and find the trial with the minimum value for chisq
@@ -319,8 +361,8 @@
             psFree(sumWeightedFlux);
 
-            if (min_k >= 0 && isfinite(minChisq) && outSrc->modelEXT) {
+            if (min_k >= 0 && isfinite(minChisq)) {
                 // copy the best fit params to the model
                 // fractional radii with the lowest chisq
-                psEllipseAxes axes = pmPSF_ModelToAxes(outSrc->modelEXT->params->data.F32, outSrc->modelEXT->type);
+                psEllipseAxes axes = pmPSF_ModelToAxes(model->params->data.F32, model->type);
 
                 // examine the params for the trial with minimum chisq.
@@ -472,5 +514,5 @@
                         axes.major = major0;
                         axes.minor = minor0;
-                        outSrc->modelEXT->chisq = chisq0;
+                        model->chisq = chisq0;
 #ifdef DUMPVECTORS
                         char fn[80];
@@ -487,5 +529,5 @@
 #endif
                     } else {
-                        outSrc->modelEXT->flags |= PM_MODEL_STATUS_NONCONVERGE;
+                        model->flags |= PM_MODEL_STATUS_NONCONVERGE;
 
 #ifdef PRINTVALS
@@ -496,12 +538,4 @@
                         // psFree(outSrc->galaxyFits);
 
-                        // XXX drop the model since it is bogus? Perhaps we do this or use the nominal values
-                        // and add a flag
-#ifdef notnow
-                        psFree(outSrc->modelEXT);
-                        outSrc->modelEXT = 0;
-                        psFree(outSrc->modelFits);
-                        outSrc->modelFits = 0;
-#endif
                     }
                     psFree(chisqFit);
@@ -516,5 +550,5 @@
                     axes.major = R_MAJOR * fRMajorBest;
                     axes.minor = R_MINOR * fRMinorBest;
-                    outSrc->modelEXT->chisq = minChisq; 
+                    model->chisq = minChisq; 
                     flux0 = fluxBest;
 #ifdef PRINTVALS
@@ -525,10 +559,10 @@
                 // now save the model parameters in the model structure
 
-                pmPSF_AxesToModel (outSrc->modelEXT->params->data.F32, axes, outSrc->modelEXT->type);
-
-                outSrc->modelEXT->mag = -2.5 * log10(flux0);
-                outSrc->modelEXT->magErr = dFlux0 / flux0; // 1 / SN
+                pmPSF_AxesToModel (model->params->data.F32, axes, modelType);
+
+                model->mag = -2.5 * log10(flux0);
+                model->magErr = dFlux0 / flux0; // 1 / SN
                 // XXX: should there be a different flag if we didn't do a fit of the chisq
-                outSrc->modelEXT->flags |= PM_MODEL_STATUS_FITTED;
+                model->flags |= PM_MODEL_STATUS_FITTED;
             }
         }
Index: trunk/psphot/src/psphotGalaxyShape.c
===================================================================
--- trunk/psphot/src/psphotGalaxyShape.c	(revision 36701)
+++ 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;
     }