Index: /trunk/psphot/src/psphotFullForceSummary.c
===================================================================
--- /trunk/psphot/src/psphotFullForceSummary.c	(revision 36730)
+++ /trunk/psphot/src/psphotFullForceSummary.c	(revision 36731)
@@ -186,5 +186,6 @@
     psFitsClose(fits);
 
-    // Convert to an array indexed by ID
+    // Convert to a set of arrays indexed by model type + 1
+    // which contain pointers to arrays indexed by ID
 #define MAX_MODEL_TYPE 10
     psArray *sortedTables = psArrayAlloc(MAX_MODEL_TYPE+1);
Index: /trunk/psphot/src/psphotFullForceSummaryReadout.c
===================================================================
--- /trunk/psphot/src/psphotFullForceSummaryReadout.c	(revision 36730)
+++ /trunk/psphot/src/psphotFullForceSummaryReadout.c	(revision 36731)
@@ -12,7 +12,7 @@
     psVector    *fRMajor;
     psVector    *fRMinor;
-    psArray *zeroPt;
-    psArray *exptime;
-    psArray *cffTables;
+    psArray *zeroPt;    // zero points for each input used with exptime to scale flux
+    psArray *exptime;   // exposure times for each input
+    psArray *cffTables; // one for each model type index is (model_type + 1) entry 0 is no extended model (star)
 } galaxyShapeOptions;
 
@@ -384,9 +384,10 @@
                 // save the flux and dFlux values from entry with lowest chisq
                 psF64 fluxBest = fluxVec->data.F32[min_k];
-                psF64 dFlux0   = dFluxVec->data.F32[min_k];
+                psF64 dFluxBest = dFluxVec->data.F32[min_k];
+                psF64 dFlux0   = NAN;
                 psF64 flux0    = NAN;
 
                 if (useFit) {
-                    #define NUM_TRIALS_INIT 9     // this isn't used for anything execpt for the initial vector lengths
+                    #define NUM_TRIALS_INIT 9
                     psVector *major = psVectorAllocEmpty(NUM_TRIALS_INIT, PS_TYPE_F64);
                     psVector *minor = psVectorAllocEmpty(NUM_TRIALS_INIT, PS_TYPE_F64);
@@ -477,4 +478,5 @@
                     psPolynomial2D *chisqFit = psPolynomial2DAlloc(PS_POLYNOMIAL_ORD, 2, 2);
                     psPolynomial2D *fluxFit = psPolynomial2DAlloc(PS_POLYNOMIAL_ORD, 2, 2);
+                    psPolynomial2D *fluxErrorFit = psPolynomial2DAlloc(PS_POLYNOMIAL_ORD, 2, 2);
                     chisqFit->coeffMask[2][2] = PS_POLY_MASK_SET;
                     chisqFit->coeffMask[2][1] = PS_POLY_MASK_SET;
@@ -483,4 +485,7 @@
                     fluxFit->coeffMask[2][1] = PS_POLY_MASK_SET;
                     fluxFit->coeffMask[1][2] = PS_POLY_MASK_SET;
+                    fluxErrorFit->coeffMask[2][2] = PS_POLY_MASK_SET;
+                    fluxErrorFit->coeffMask[2][1] = PS_POLY_MASK_SET;
+                    fluxErrorFit->coeffMask[1][2] = PS_POLY_MASK_SET;
 
                     bool goodFit = psVectorFitPolynomial2D (chisqFit, NULL, 0xff, chisq, NULL, major, minor);
@@ -504,4 +509,12 @@
                         } else {
                             flux0 = fluxBest;
+                        }
+
+                        // .. and compute dFlux at the minimum chisq position
+                        bool goodFluxErrorFit = psVectorFitPolynomial2D(fluxErrorFit, NULL, 0xFF, dFlux, NULL, major, minor);
+                        if (goodFluxErrorFit) {
+                            dFlux0  = psPolynomial2DEval(fluxErrorFit, major0, minor0);
+                        } else {
+                            dFlux0 = dFluxBest;
                         }
 
@@ -540,4 +553,6 @@
                     }
                     psFree(chisqFit);
+                    psFree(fluxFit);
+                    psFree(fluxErrorFit);
                     psFree(major);
                     psFree(minor);
@@ -552,4 +567,5 @@
                     model->chisq = minChisq; 
                     flux0 = fluxBest;
+                    dFlux0 = dFluxBest;
 #ifdef PRINTVALS
                     fprintf(stderr, "%4d %3d %3ld | %6.3f %6.3f %6.3f %6.3f skip fit\n", 
@@ -581,5 +597,5 @@
 static bool setOptions(galaxyShapeOptions *opt, pmReadout *readout, psMetadata *recipe, bool makeVectors) {
     bool status;
-    bool useAnalysis;
+    bool useAnalysis;   // fall back to recipe if we dont' find values in analysis. Probably should no longer do this
 
     psMetadataLookupF32(&useAnalysis, readout->analysis, "GALAXY_SHAPES_FR_MAJOR_MIN");