Index: trunk/ppStack/src/ppStackArguments.c
===================================================================
--- trunk/ppStack/src/ppStackArguments.c	(revision 31054)
+++ trunk/ppStack/src/ppStackArguments.c	(revision 31158)
@@ -110,4 +110,10 @@
 {
     assert(config);
+
+    // generic arguments (version, dumpconfig)
+    PS_ARGUMENTS_GENERIC( ppStack, config, argc, argv );
+
+    // thread arguments
+    PS_ARGUMENTS_THREADS( ppStack, config, argc, argv )
 
     {
@@ -181,5 +187,4 @@
     psMetadataAddStr(arguments, PS_LIST_TAIL, "-temp-variance", 0, "Suffix for temporary variance maps", NULL);
     psMetadataAddBool(arguments, PS_LIST_TAIL, "-temp-delete", 0, "Delete temporary files on completion?", false);
-    psMetadataAddS32(arguments, PS_LIST_TAIL, "-threads", 0, "Number of threads to use", 0);
     psMetadataAddBool(arguments, PS_LIST_TAIL, "-visual", 0, "visualisation", false);
 
@@ -238,10 +243,4 @@
     valueArgStr(arguments, "-stats", "STATS", arguments);
 
-    int numThreads = psMetadataLookupS32(NULL, arguments, "-threads"); // Number of threads
-    if (numThreads > 0 && !psThreadPoolInit(numThreads)) {
-        psError(PPSTACK_ERR_ARGUMENTS, false, "Unable to setup %d threads", numThreads);
-        return false;
-    }
-
     psMetadataAddBool(arguments, PS_LIST_TAIL, "PPSTACK.DEBUG.STACK", 0,
                       "Read old convolved images to debug stack?", debugStack);
Index: trunk/ppStack/src/ppStackCombineFinal.c
===================================================================
--- trunk/ppStack/src/ppStackCombineFinal.c	(revision 31054)
+++ trunk/ppStack/src/ppStackCombineFinal.c	(revision 31158)
@@ -61,5 +61,5 @@
         }
 
-        // call: ppStackReadoutFinal(config, outRO, readouts, rejected)
+        // calls ppStackReadoutFinal(config, outRO, readouts, rejected) in ppStackReadout.c
         psThreadJob *job = psThreadJobAlloc("PPSTACK_FINAL_COMBINE"); // Job to start
         psArrayAdd(job->args, 1, thread);
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) {
Index: trunk/ppStack/src/ppStackFiles.c
===================================================================
--- trunk/ppStack/src/ppStackFiles.c	(revision 31054)
+++ trunk/ppStack/src/ppStackFiles.c	(revision 31158)
@@ -65,6 +65,5 @@
     for (int i = 0; i < numLeaks; i++) {
         psMemBlock *mb = leaks[i];
-        fprintf(memFile, "%12lu\t%12zd\t%s:%d\n", mb->id, mb->userMemorySize,
-                mb->file, mb->lineno);
+        fprintf(memFile, "%12lu\t%12zd\t%s:%d\n", mb->id, mb->userMemorySize, mb->file, mb->lineno);
         total += mb->userMemorySize;
     }
@@ -75,6 +74,4 @@
     num++;
 }
-
-
 
 // Activate/deactivate a list of files
Index: trunk/ppStack/src/ppStackLoop.c
===================================================================
--- trunk/ppStack/src/ppStackLoop.c	(revision 31054)
+++ trunk/ppStack/src/ppStackLoop.c	(revision 31158)
@@ -55,6 +55,5 @@
     }
 
-    // Start threading
-    ppStackThreadInit();
+    // Define threading elements
     ppStackThreadData *stack = ppStackThreadDataSetup(options, config, true);
     if (!stack) {
@@ -114,5 +113,5 @@
     }
 
-    // Final combination
+    // Final combination.  This one does NOT need to be normalized since the convolution takes care of that
     psTrace("ppStack", 2, "Final stack of convolved images....\n");
     if (!ppStackCombineFinal(stack, options->convCovars, options, config, false, false, true)) {
@@ -181,8 +180,7 @@
         }
 
-	// generate the unconvolved stack
+	// generate the unconvolved stack. NOTE: this one must be normalized since the inputs have not been
         psTrace("ppStack", 2, "Stack of unconvolved images....\n");
-        if (!ppStackCombineFinal(stack, options->origCovars, options, config,
-                                 false, true, false)) {
+        if (!ppStackCombineFinal(stack, options->origCovars, options, config, false, true, false)) {
             psError(psErrorCodeLast(), false, "Unable to perform unconvolved combination.");
             psFree(stack);
Index: trunk/ppStack/src/ppStackMatch.c
===================================================================
--- trunk/ppStack/src/ppStackMatch.c	(revision 31054)
+++ trunk/ppStack/src/ppStackMatch.c	(revision 31158)
@@ -142,5 +142,5 @@
     bool mdok;                          // Status of MD lookup
     float penalty = psMetadataLookupF32(NULL, ppsub, "PENALTY"); // Penalty for wideness
-    int threads = psMetadataLookupS32(NULL, config->arguments, "-threads"); // Number of threads
+    int threads = psMetadataLookupS32(NULL, config->arguments, "NTHREADS"); // Number of threads
 
     // Replaced pmReadoutMaskNonfinite with pmReadoutMaskInvalid (tests for already masked pixels)
Index: trunk/ppStack/src/ppStackOptions.h
===================================================================
--- trunk/ppStack/src/ppStackOptions.h	(revision 31054)
+++ trunk/ppStack/src/ppStackOptions.h	(revision 31158)
@@ -22,7 +22,13 @@
     float sumExposure;                  // Sum of exposure times
     float zp;                           // Zero point for output
+    float zpErr;			// Zero point error for output
+    float airmass;			// airmass for output image
+    float airmassSlope;			// airmass slope used for zp analysis
     psVector *inputMask;                // Mask for inputs
     psArray *sourceLists;               // Individual lists of sources for matching
     psVector *norm;                     // Normalisation for each image
+    psVector *zpInput;			// reported zero point of input exposure
+    psVector *expTimeInput;		// exposure time of input exposure
+    psVector *airmassInput;		// airmass of input exposure
     psArray *sources;                   // Matched sources
     // Convolve
Index: trunk/ppStack/src/ppStackPrepare.c
===================================================================
--- trunk/ppStack/src/ppStackPrepare.c	(revision 31054)
+++ trunk/ppStack/src/ppStackPrepare.c	(revision 31158)
@@ -1,7 +1,9 @@
-#include "ppStack.h"
-
-#define PHOT_SOURCE_MASK (PM_SOURCE_MODE_FAIL | PM_SOURCE_MODE_SATSTAR | PM_SOURCE_MODE_BLEND | \
-                          PM_SOURCE_MODE_BADPSF | PM_SOURCE_MODE_DEFECT | PM_SOURCE_MODE_SATURATED | \
-                          PM_SOURCE_MODE_CR_LIMIT | PM_SOURCE_MODE_EXT_LIMIT) // Mask to apply to sources
+# include "ppStack.h"
+
+# define RE_PHOTOMETER 0
+
+# define PHOT_SOURCE_MASK (PM_SOURCE_MODE_FAIL | PM_SOURCE_MODE_SATSTAR | PM_SOURCE_MODE_BLEND | \
+			   PM_SOURCE_MODE_BADPSF | PM_SOURCE_MODE_DEFECT | PM_SOURCE_MODE_SATURATED | \
+			   PM_SOURCE_MODE_CR_LIMIT | PM_SOURCE_MODE_EXT_LIMIT) // Mask to apply to sources
 
 bool ppStackInputPhotometer(const pmReadout *ro, const psArray *sources, const pmConfig *config)
@@ -13,4 +15,5 @@
     psAssert(recipe, "We've thrown an error on this before.");
 
+# if (RE_PHOTOMETER)
     bool mdok;                          // Status of MD lookup
 
@@ -38,7 +41,6 @@
     psImageMaskType maskVal = pmConfigMaskGet(maskValStr, config); // Bits to mask
 
+    // Measure background
     psImage *image = ro->image, *mask = ro->mask; // Image and mask from readout
-
-    // Measure background
     psRandom *rng = psRandomAlloc(PS_RANDOM_TAUS); // Random number generator
     psStats *stats = psStatsAlloc(PS_STAT_ROBUST_MEDIAN | PS_STAT_ROBUST_STDEV); // Statistics
@@ -52,7 +54,7 @@
     psFree(stats);
     psFree(rng);
-
-    // Photometer sources
-    int numCols = image->numCols, numRows = image->numRows; // Size of image
+# endif
+
+    // Examine the sources (photometer if RE_PHOTOMETER is set)
     int numSources = sources->n; // Number of sources
     int numGood = 0;            // Number of good sources
@@ -64,5 +66,13 @@
             continue;
         }
-        float x = source->peak->xf, y = source->peak->yf; // Coordinates of source
+
+// This code re-measures the magnitude in the assumption that the warp analysis failed It would
+// be better to catch and correct such cases.  In any case, we should check on the validity of
+// this assumption (eg, compare the mean psf-ap offset) and raise an error on that?
+
+# if (RE_PHOTOMETER)
+	
+	int numCols = image->numCols, numRows = image->numRows; // Size of image
+	float x = source->peak->xf, y = source->peak->yf; // Coordinates of source
         int xCentral = x + 0.5, yCentral = y + 0.5; // Central pixel
         // Range of integration
@@ -96,4 +106,8 @@
             }
         }
+# else
+	numGood++;
+	maxMag = PS_MAX(maxMag, source->psfMag);
+# endif
     }
     psLogMsg("ppStack", PS_LOG_INFO, "Photometered %d good sources in image; max mag %f", numGood, maxMag);
@@ -217,4 +231,8 @@
         }
     }
+
+    bool mdok = false;
+    float maxFWHM = psMetadataLookupF32(&mdok, recipe, "PSF.INPUT.MAX"); // max allowed input fwhm
+    float clipFWHMnSig = psMetadataLookupF32(&mdok, recipe, "PSF.INPUT.CLIP.NSIGMA"); // sigma clipping of inputs
 
     psString log = psStringCopy("Input seeing FWHMs:\n"); // Log message
@@ -251,4 +269,10 @@
         }
 
+	// reject any input images which exceed the specified max FWHM
+	if (isfinite(maxFWHM) && (options->inputSeeing->data.F32[i] > maxFWHM)) {
+            options->inputMask->data.PS_TYPE_VECTOR_MASK_DATA[i] = PPSTACK_MASK_PSF;
+            psLogMsg("ppStack", PS_LOG_INFO, "PSF FWHM for image %d is too large (%f vs %f maxFWHM) --- rejected.", i, options->inputSeeing->data.F32[i], maxFWHM);
+	}
+
         psStringAppend(&log, "Input %d: %f\n", i, options->inputSeeing->data.F32[i]);
     }
@@ -257,4 +281,25 @@
     }
     psFree(log);
+
+    // We should have the ability to filter the input list based on the seeing:
+    // * reject above some max value and/or min value
+    // * reject images out-of-line with the rest of the images
+
+    // measure stats
+    psStats *fwhmStats = psStatsAlloc(PS_STAT_CLIPPED_MEAN | PS_STAT_CLIPPED_STDEV);
+    psVectorStats (fwhmStats, options->inputSeeing, NULL, options->inputMask, 0xff);
+    psLogMsg("ppStack", PS_LOG_INFO, "Input FWHMs : %f +/- %f", fwhmStats->clippedMean, fwhmStats->clippedStdev);
+
+    if (isfinite(clipFWHMnSig)) {
+	float fwhmLimit = fwhmStats->clippedMean + clipFWHMnSig * fwhmStats->clippedStdev;
+	fwhmLimit = isfinite(maxFWHM) ? PS_MIN (maxFWHM, fwhmLimit) : fwhmLimit;
+	for (int i = 0; i < num; i++) {
+            if (options->inputMask->data.PS_TYPE_VECTOR_MASK_DATA[i]) continue;
+            if (options->inputSeeing->data.F32[i] > fwhmLimit) {
+		options->inputMask->data.PS_TYPE_VECTOR_MASK_DATA[i] = PPSTACK_MASK_PSF;
+		psLogMsg("ppStack", PS_LOG_INFO, "PSF FWHM for image %d is too large (%f vs %f fwhmLimit) --- rejected.", i, options->inputSeeing->data.F32[i], fwhmLimit);
+	    }
+        }
+    }
 
     // Generate target PSF
Index: trunk/ppStack/src/ppStackReject.c
===================================================================
--- trunk/ppStack/src/ppStackReject.c	(revision 31054)
+++ trunk/ppStack/src/ppStackReject.c	(revision 31158)
@@ -55,5 +55,5 @@
 
     // Reject bad pixels
-    if (psMetadataLookupS32(NULL, config->arguments, "-threads") > 0) {
+    if (psMetadataLookupS32(NULL, config->arguments, "NTHREADS") > 0) {
         pmStackRejectThreadsInit();
     }
Index: trunk/ppStack/src/ppStackSetup.c
===================================================================
--- trunk/ppStack/src/ppStackSetup.c	(revision 31054)
+++ trunk/ppStack/src/ppStackSetup.c	(revision 31158)
@@ -51,5 +51,5 @@
     psString outputName = psStringCopy(psMetadataLookupStr(NULL, config->arguments,
                                                            "OUTPUT")); // Name for temporary files
-    const char *tempName = basename(outputName);
+    const char *tempName = psStringFileBasename(outputName);
     if (!tempName) {
         psError(PPSTACK_ERR_ARGUMENTS, false, "Unable to construct basename for temporary files.");
@@ -65,4 +65,5 @@
                 "Unable to find TEMP.IMAGE, TEMP.MASK and TEMP.VARIANCE in recipe");
         psFree(outputName);
+	psFree(tempName);
         return false;
     }
@@ -82,4 +83,5 @@
     }
     psFree(outputName);
+    psFree(tempName);
 
     // Original images
Index: trunk/ppStack/src/ppStackSources.c
===================================================================
--- trunk/ppStack/src/ppStackSources.c	(revision 31054)
+++ trunk/ppStack/src/ppStackSources.c	(revision 31158)
@@ -115,5 +115,11 @@
     float fracMatch = psMetadataLookupF32(NULL, recipe, "ZP.MATCH"); // Fraction of images to match for star
 
-    psMetadata *airmassZP = psMetadataLookupMetadata(NULL, recipe, "ZP.AIRMASS"); // Airmass terms
+    bool mdok = false;
+    float airmassTarget = psMetadataLookupF32(&mdok, recipe, "ZP.AIRMASS.TARGET"); // output airmass value 
+    if (!mdok) {
+	airmassTarget = 1.0;
+    }
+
+    psMetadata *airmassZP = psMetadataLookupMetadata(NULL, recipe, "ZP.AIRMASS"); // Airmass terms (slopes) by filter
     if (!airmassZP) {
         psError(PPSTACK_ERR_CONFIG, false, "Unable to find ZP.AIRMASS in recipe.");
@@ -128,4 +134,9 @@
     int num = options->num;             // Number of inputs
     psAssert(num == sourceLists->n, "Wrong number of source lists: %ld\n", sourceLists->n);
+
+    // vectors to store these inputs so they may be recorded in the output headers
+    options->zpInput      = psVectorAlloc(num, PS_TYPE_F32);
+    options->expTimeInput = psVectorAlloc(num, PS_TYPE_F32);
+    options->airmassInput = psVectorAlloc(num, PS_TYPE_F32);
 
     psVector *zp = psVectorAlloc(num, PS_TYPE_F32); // Relative zero points for each image
@@ -167,4 +178,10 @@
         const char *expFilter = psMetadataLookupStr(NULL, file->fpa->concepts, "FPA.FILTER"); // Filter name
         zpExp->data.F32[i] = psMetadataLookupF32(NULL, file->fpa->concepts, "FPA.ZP"); // Exposure zero point
+	// XXX need to get the zero point error values and propagate to get the FPA.ZP.ERR value
+
+	options->zpInput->data.F32[i] = zpExp->data.F32[i]; // NOTE zpExp may be re-assigned below using relative photometry
+	options->expTimeInput->data.F32[i] = exptime;
+	options->airmassInput->data.F32[i] = airmass;
+
         psLogMsg("ppStack", PS_LOG_INFO,
                  "Image %d: %.2f sec exposure in %s at airmass %.2f with zero point %.2f",
@@ -185,5 +202,4 @@
         if (!filter) {
             filter = expFilter;
-            bool mdok;
             airmassTerm = psMetadataLookupF32(&mdok, airmassZP, filter);
             if (!mdok || !isfinite(airmassTerm)) {
@@ -206,4 +222,6 @@
         }
 
+	// XXX this is wrong, or at least inconsistent with the above: this needs to include 
+	// a value for the nominal system zero point to be consistent with zpExp
         zp->data.F32[i] = airmassTerm * airmass + 2.5 * log10(exptime);
     }
@@ -236,6 +254,12 @@
 
     if (zpExpNum == numGoodImages) {
+	psLogMsg ("ppStack", PS_LOG_INFO, "all zero points are finite; using reported zero points listed above");
         for (int i = 0; i < num; i++) {
             zp->data.F32[i] = zpExp->data.F32[i];
+        }
+    } else {
+	psLogMsg ("ppStack", PS_LOG_INFO, "missing some zero points; using guess values:");
+        for (int i = 0; i < num; i++) {
+	    psLogMsg("ppStack", PS_LOG_INFO, "Image %d: %.2f sec exposure with zero point %.2f", i, options->exposures->data.F32[i], zp->data.F32[i]);
         }
     }
@@ -265,9 +289,40 @@
             }
         }
-        // M = m + c0 + c1 * airmass - 2.5log(t) + transparency
-        // Want sources to have m corresponding to airmass = 1 and t = sumExpTime and transparency = 0
-        // m_0 + c1 * airmass_0 + 2.5log(t_0) - trans_0 = m_1 + c1 * airmass_1 + 2.5log(t_1) - trans_1
-        // m_0 = m_1 + zp_1 + trans_1 - c1 * airmass_0 - 2.5log(t_0)
-        // We don't need to know the magnitude zero point for the filter, since it cancels out
+
+	// EAM : the discussion here was not quite right (or at least sloppy).  Here is a replacement explanation:
+
+	// For any star, the observed instrumental magnitude on an image and the apparent magnitude are related by:
+	// M_app = m_inst + zp + c1 * airmass + 2.5log(t) - transparency
+	// NOTE the sign of 'transparency'  this must agree with the definition in pmSourceMatch.c. see, eg, line 457 where 
+	// transparency = m_inst + zp + c1 * airmass + 2.5log(t) - M_app 
+
+	// we want to adjust the input images to be in a consistent flux system so that the
+	// final stack can be generated with a specific target zero point.  Any adjustment to
+	// the flux scale of the image must be made in coordination with the resulting
+	// zeropoint, exposure time, and airmass such that the above relationship yields the
+	// same apparent magnitude for a given star:
+
+	// m_inst_i : instrumental mags on input image (in)
+	// m_inst_o : instrumental mags on re-normalized image (out)
+
+	// m_inst_o + zp_o + c1 * airmass_o + 2.5log(t_o) - trans_o = m_inst_i + zp_i + c1 * airmass_i + 2.5log(t_i) - trans_i
+
+	// m_inst_o = m_inst_i + (zp_i - zp_o) + c1 * (airmass_i - airmass_o) + 2.5log(t_i) - 2.5log(t_o) - trans_i + trans_o
+
+	// zp_i, airmass_i, t_i, trans_i : reported or measured for input image
+
+	// zp_o      = zpTarget      (from recipe)
+	// airmass_o = airmassTarget (from recipe)
+	// t_o       = sumExpTime    [sum of input exposure times: once images are scale to this time, they can be avereaged]
+	// trans_o   = 0.0           [obviously!]
+
+	// we have 2 cases: (a) all reported ZPs are good or (b) some are bad:
+	// (a) FPA.ZP = zp_i + c1 * airmass_i
+	//  --> zp[i] = zp_i + c1 * airmass_i + 2.5log(exptime_i)
+	// (b)  zp[i] = c1 * airmass_i + 2.5log(exptime_i)
+	// NOTE: in case (b), the current code is equating the TARGET zp with the NOMINAL zp, which is wrong.
+
+	// m_inst_o - m_inst_i = zp[i] - zpTarget - c1 * airmassTarget - 2.5log(sumExpTime) - trans_i
+
         if (options->matchZPs) {
             options->norm = psVectorAlloc(num, PS_TYPE_F32);
@@ -277,7 +332,9 @@
                 }
                 psArray *sources = sourceLists->data[i]; // Sources of interest
-                float magCorr = zp->data.F32[i] - trans->data.F32[i] - 2.5*log10(options->sumExposure);
-                if (zpExpNum == numGoodImages) {
+                float magCorr = zp->data.F32[i] - trans->data.F32[i] - 2.5*log10(options->sumExposure) - airmassTerm * airmassTarget;
+                if (zpExpNum == numGoodImages) { // case (a)
                     // Using measured zero points, so attempt to set target zero point
+		    // XXX see NOTE above regarding case (b) : this is wrong.  the code should load a nominal zero point and supply it above
+		    // 
                     magCorr -= zpTarget;
                 }
@@ -302,4 +359,6 @@
             // Producing image with target zero point
             options->zp = zpTarget;
+            options->airmass = airmassTarget;
+            options->airmassSlope = airmassTerm;
         } else {
             options->zp = NAN;
Index: trunk/ppStack/src/ppStackThread.c
===================================================================
--- trunk/ppStack/src/ppStackThread.c	(revision 31054)
+++ trunk/ppStack/src/ppStackThread.c	(revision 31158)
@@ -101,5 +101,5 @@
     }
 
-    int numThreads = psMetadataLookupS32(NULL, config->arguments, "-threads"); // Number of threads
+    int numThreads = psMetadataLookupS32(NULL, config->arguments, "NTHREADS"); // Number of threads
 
     // Generate readouts for each input file in each file group
@@ -248,5 +248,5 @@
 
 
-void ppStackThreadInit(void)
+void ppStackSetThreads(void)
 {
     static bool threaded = false;       // Are we running threaded?
Index: trunk/ppStack/src/ppStackThread.h
===================================================================
--- trunk/ppStack/src/ppStackThread.h	(revision 31054)
+++ trunk/ppStack/src/ppStackThread.h	(revision 31158)
@@ -43,5 +43,5 @@
 
 // Initialise the threads
-void ppStackThreadInit(void);
+void ppStackSetThreads(void);
 
 
Index: trunk/ppStack/src/ppStackUpdateHeader.c
===================================================================
--- trunk/ppStack/src/ppStackUpdateHeader.c	(revision 31054)
+++ trunk/ppStack/src/ppStackUpdateHeader.c	(revision 31158)
@@ -70,4 +70,58 @@
     psMetadataAddStr(hdu->header, PS_LIST_TAIL, "TESS_ID", PS_META_REPLACE, "type of stack", tessID);
 
+    psMetadataAddF32(hdu->header, PS_LIST_TAIL, "AIRM_SLP", PS_META_REPLACE, "airmass slope", options->airmassSlope);
+    
+
+    // write the following fields to the output headers:
+    // INP_NNNN : input file names
+    // SCL_NNNN : scale factor for each input file
+    // ZPT_NNNN : zero point of input file 
+    // EXT_NNNN : exptime of input file
+    // AIR_NNNN : airmass of input file
+
+    // use separate loops so they are blocked together in the headers (more readable)
+
+    // save the input filenames
+    for (int i = 0; i < options->num; i++) {
+	char field[64];
+	snprintf (field, 64, "INP_%04d", i);
+
+	psString basename = psStringFileBasename(options->origImages->data[i]);
+	psMetadataAddStr(hdu->header, PS_LIST_TAIL, field, PS_META_REPLACE, "input image name", basename);
+	psFree(basename);
+    }	
+
+    // save the input normalizations
+    for (int i = 0; i < options->num; i++) {
+	char field[64];
+
+	float value = options->norm ? pow(10.0, -0.4*options->norm->data.F32[i]) : NAN;
+	snprintf (field, 64, "SCL_%04d", i);
+	psMetadataAddF32(hdu->header, PS_LIST_TAIL, field, PS_META_REPLACE, "input image scale factor", value);
+    }	
+    // save the input exptimes
+    for (int i = 0; i < options->num; i++) {
+	char field[64];
+
+	float value = options->zpInput ? options->zpInput->data.F32[i] : NAN;
+	snprintf (field, 64, "ZPT_%04d", i);
+	psMetadataAddF32(hdu->header, PS_LIST_TAIL, field, PS_META_REPLACE, "input image zero point", value);
+    }	
+    // save the input normalizations
+    for (int i = 0; i < options->num; i++) {
+	char field[64];
+
+	float value = options->expTimeInput ? options->expTimeInput->data.F32[i] : NAN;
+	snprintf (field, 64, "EXP_%04d", i);
+	psMetadataAddF32(hdu->header, PS_LIST_TAIL, field, PS_META_REPLACE, "input image exptime", value);
+    }	
+    // save the input normalizations
+    for (int i = 0; i < options->num; i++) {
+	char field[64];
+
+	float value = options->airmassInput ? options->airmassInput->data.F32[i] : NAN;
+	snprintf (field, 64, "AIR_%04d", i);
+	psMetadataAddF32(hdu->header, PS_LIST_TAIL, field, PS_META_REPLACE, "input image airmass", value);
+    }	
     return true;
 }