Index: /branches/eam_branches/ipp-20110213/ppStack/src/ppStackConvolve.c
===================================================================
--- /branches/eam_branches/ipp-20110213/ppStack/src/ppStackConvolve.c	(revision 30800)
+++ /branches/eam_branches/ipp-20110213/ppStack/src/ppStackConvolve.c	(revision 30801)
@@ -4,10 +4,10 @@
 
 // Update the value of a concept
-#define UPDATE_CONCEPT(SOURCE, NAME, VALUE) { \
-    psMetadataItem *item = psMetadataLookup(SOURCE->concepts, NAME); \
-    psAssert(item, "Concept should be present"); \
-    psAssert(item->type == PS_TYPE_F32, "Concept should be F32"); \
-    item->data.F32 = VALUE; \
-}
+#define UPDATE_CONCEPT(SOURCE, NAME, VALUE) {				\
+	psMetadataItem *item = psMetadataLookup(SOURCE->concepts, NAME); \
+	psAssert(item, "Concept should be present");			\
+	psAssert(item->type == PS_TYPE_F32, "Concept should be F32");	\
+	item->data.F32 = VALUE;						\
+    }
 
 bool ppStackConvolve(ppStackOptions *options, pmConfig *config)
@@ -45,4 +45,6 @@
     psVectorInit(renorms, NAN);
 
+    psVector *satValues = psVectorAllocEmpty(num, PS_TYPE_F32); // Renormalisation values for variances
+
     psList *fpaList = psListAlloc(NULL); // List of input FPAs, for concept averaging
     psList *cellList = psListAlloc(NULL); // List of input cells, for concept averaging
@@ -68,4 +70,6 @@
             psFree(cellList);
             psFree(target);
+	    psFree(renorms);
+	    psFree(satValues);
             return false;
         }
@@ -85,4 +89,6 @@
             psFree(cellList);
             psFree(target);
+	    psFree(renorms);
+	    psFree(satValues);
             return false;
         }
@@ -104,4 +110,6 @@
                 psFree(cellList);
                 psFree(target);
+		psFree(renorms);
+		psFree(satValues);
                 return false;
                 // Non-fatal errors
@@ -118,4 +126,6 @@
                     psFree(cellList);
                     psFree(target);
+		    psFree(renorms);
+		    psFree(satValues);
                     return false;
                 }
@@ -164,4 +174,6 @@
             psFree(rng);
             psFree(target);
+	    psFree(renorms);
+	    psFree(satValues);
             return false;
         }
@@ -175,4 +187,6 @@
             psFree(maskHeader);
             psFree(target);
+	    psFree(renorms);
+	    psFree(satValues);
             return false;
         }
@@ -184,4 +198,6 @@
             psFree(rng);
             psFree(target);
+	    psFree(renorms);
+	    psFree(satValues);
             return false;
         }
@@ -219,4 +235,6 @@
         if (options->matchZPs) {
             // I think I need to take off the exposure time because we're going to set the new exposure time
+	    // Clarification: the zero point (ZP) in the header should be set such that:
+	    // M_app = m_inst + ZP + 2.5*log(exptime), where exptime in the output is sumExposure
             psMetadataItem *zpItem = psMetadataLookup(inCell->parent->parent->concepts, "FPA.ZP");
             zpItem->data.F32 += options->norm->data.F32[i] + 2.5*log10(options->sumExposure);
@@ -227,4 +245,9 @@
             expItem = psMetadataLookup(inCell->concepts, "CELL.EXPOSURE");
             expItem->data.F32 = options->sumExposure;
+
+	    // flux_out = flux_in * ten(-0.4*norm) -- save the individual saturation values
+            psMetadataItem *satItem = psMetadataLookup(inCell->concepts, "CELL.EXPOSURE");
+	    satItem->data.F32 *= pow(10.0, -0.4*options->norm->data.F32[i]);
+            psVectorAppend (satValues, satItem->data.F32);
         }
 
@@ -235,4 +258,6 @@
             psFree(rng);
             psFree(target);
+	    psFree(renorms);
+	    psFree(satValues);
             return false;
         }
@@ -255,4 +280,6 @@
         psFree(fpaList);
         psFree(cellList);
+	psFree(renorms);
+	psFree(satValues);
         return true;
     }
@@ -279,8 +306,18 @@
         psFree(fpaList);
         psFree(cellList);
-
+	
+	// The best guess for an output saturation value depends on the recipe.  If we have
+	// 'safe' on, the we require at least 2 pixels to generate a valid output pixel.  In
+	// this case, the best value for CELL.SATURATION is the 2nd highest value in the list.
+	// If not, it should be the higest value in the list
+	bool safe = psMetadataLookupBool(&mdok, recipe, "SAFE"); // Be safe when combining small numbers of pixels
+	psVectorSortInPlace(satValues);
+	float satBest = safe && satValues->n > 1 ? satValues->data.F32[1] : satValues->data.F32[0];
+
+	// UPDATE CELL.SATURATION here
         UPDATE_CONCEPT(outFPA,  "FPA.EXPOSURE",  options->sumExposure);
         UPDATE_CONCEPT(outCell, "CELL.EXPOSURE", options->sumExposure);
         UPDATE_CONCEPT(outCell, "CELL.DARKTIME", NAN);
+        UPDATE_CONCEPT(outCell, "CELL.SATURATION", satBest);
         UPDATE_CONCEPT(outFPA,  "FPA.ZP",        options->zp);
 
@@ -288,4 +325,5 @@
         UPDATE_CONCEPT(unconvCell, "CELL.EXPOSURE", options->sumExposure);
         UPDATE_CONCEPT(unconvCell, "CELL.DARKTIME", NAN);
+        UPDATE_CONCEPT(unconvCell, "CELL.SATURATION", satBest);
         UPDATE_CONCEPT(unconvFPA,  "FPA.ZP",        options->zp);
 
@@ -295,9 +333,14 @@
             double time = psTimeToMJD(fpaTime);
             psMetadataAddF64(options->stats, PS_LIST_TAIL, "MJD_OBS", PS_META_DUPLICATE_OK,
-                                 "Average MJD_OBS of inputs", time);
-
-        }
-    }
-
+			     "Average MJD_OBS of inputs", time);
+
+        }
+    }
+
+    // XXX EAM : this may be overly harsh -- or at least it would be if I (EAM) hadn't changed
+    // the values of matchChi2 I modified pmSubtraction.c to fit a 2nd order polynomial to the
+    // star chisq distribution (because of systematic errors in the model being matched).  This
+    // fit forces the mean value to be 0.0.  Perhaps we can / should exclude images which have
+    // an excessively high value for the rms?
 
     // Reject images out-of-hand on the basis of their match chi^2
@@ -314,4 +357,6 @@
             psError(PPSTACK_ERR_PROG, false, "Unable to sort vector.");
             psFree(values);
+	    psFree(renorms);
+	    psFree(satValues);
             return false;
         }
@@ -354,4 +399,6 @@
         psErrorClear();
         psWarning("No good images survived convolution stage.");
+	psFree(renorms);
+	psFree(satValues);
         return true;
     }
@@ -364,4 +411,5 @@
     }
     psFree(renorms);
+    psFree(satValues);
 
     return true;