Index: trunk/ppStack/src/ppStackConvolve.c
===================================================================
--- trunk/ppStack/src/ppStackConvolve.c	(revision 31054)
+++ trunk/ppStack/src/ppStackConvolve.c	(revision 31158)
@@ -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)
@@ -47,4 +47,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
@@ -70,4 +72,6 @@
             psFree(cellList);
             psFree(target);
+	    psFree(renorms);
+	    psFree(satValues);
             return false;
         }
@@ -87,4 +91,6 @@
             psFree(cellList);
             psFree(target);
+	    psFree(renorms);
+	    psFree(satValues);
             return false;
         }
@@ -106,4 +112,6 @@
                 psFree(cellList);
                 psFree(target);
+		psFree(renorms);
+		psFree(satValues);
                 return false;
                 // Non-fatal errors
@@ -120,4 +128,6 @@
                     psFree(cellList);
                     psFree(target);
+		    psFree(renorms);
+		    psFree(satValues);
                     return false;
                 }
@@ -166,4 +176,6 @@
             psFree(rng);
             psFree(target);
+	    psFree(renorms);
+	    psFree(satValues);
             return false;
         }
@@ -177,4 +189,6 @@
             psFree(maskHeader);
             psFree(target);
+	    psFree(renorms);
+	    psFree(satValues);
             return false;
         }
@@ -186,4 +200,6 @@
             psFree(rng);
             psFree(target);
+	    psFree(renorms);
+	    psFree(satValues);
             return false;
         }
@@ -221,5 +237,8 @@
         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");
+	    float inZP = zpItem->data.F32;
             zpItem->data.F32 += options->norm->data.F32[i] + 2.5*log10(options->sumExposure);
 
@@ -228,5 +247,15 @@
 
             expItem = psMetadataLookup(inCell->concepts, "CELL.EXPOSURE");
+	    float inExptime = expItem->data.F32;
             expItem->data.F32 = options->sumExposure;
+
+	    // flux_out = flux_in * ten(-0.4*norm) -- save the individual saturation values
+            psMetadataItem *satItem = psMetadataLookup(inCell->concepts, "CELL.SATURATION");
+	    float inSat = satItem->data.F32;
+	    satItem->data.F32 *= pow(10.0, -0.4*options->norm->data.F32[i]);
+            psVectorAppend (satValues, satItem->data.F32);
+
+	    psLogMsg("ppStack", PS_LOG_INFO, "image %d mods : zp %f -> %f, exptime %f -> %f, sat %f -> %f", 
+		     i, inZP, zpItem->data.F32, inExptime, expItem->data.F32, inSat, satItem->data.F32);
         }
 
@@ -237,4 +266,6 @@
             psFree(rng);
             psFree(target);
+	    psFree(renorms);
+	    psFree(satValues);
             return false;
         }
@@ -257,4 +288,6 @@
         psFree(fpaList);
         psFree(cellList);
+	psFree(renorms);
+	psFree(satValues);
         return true;
     }
@@ -281,14 +314,30 @@
         psFree(fpaList);
         psFree(cellList);
-
-        UPDATE_CONCEPT(outFPA,  "FPA.EXPOSURE",  options->sumExposure);
-        UPDATE_CONCEPT(outCell, "CELL.EXPOSURE", options->sumExposure);
-        UPDATE_CONCEPT(outCell, "CELL.DARKTIME", NAN);
-        UPDATE_CONCEPT(outFPA,  "FPA.ZP",        options->zp);
-
-        UPDATE_CONCEPT(unconvFPA,  "FPA.EXPOSURE",  options->sumExposure);
-        UPDATE_CONCEPT(unconvCell, "CELL.EXPOSURE", options->sumExposure);
-        UPDATE_CONCEPT(unconvCell, "CELL.DARKTIME", NAN);
-        UPDATE_CONCEPT(unconvFPA,  "FPA.ZP",        options->zp);
+	
+	// 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 mdok = false;
+	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);
+        UPDATE_CONCEPT(outFPA,  "FPA.AIRMASS",     options->airmass);
+
+        UPDATE_CONCEPT(unconvFPA,  "FPA.EXPOSURE",    options->sumExposure);
+        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);
+        UPDATE_CONCEPT(unconvFPA,  "FPA.AIRMASS",     options->airmass);
+
+	psLogMsg("ppStack", PS_LOG_INFO, "stack adjust metadata values : zp %f, exptime %f, sat %f", options->zp, options->sumExposure, satBest);
 
         if (options->stats) {
@@ -297,9 +346,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
@@ -316,4 +370,6 @@
             psError(PPSTACK_ERR_PROG, false, "Unable to sort vector.");
             psFree(values);
+	    psFree(renorms);
+	    psFree(satValues);
             return false;
         }
@@ -356,4 +412,6 @@
         psErrorClear();
         psWarning("No good images survived convolution stage.");
+	psFree(renorms);
+	psFree(satValues);
         return true;
     }
@@ -366,4 +424,5 @@
     }
     psFree(renorms);
+    psFree(satValues);
 
     if (options->stats) {