Index: branches/meh_branches/ppstack_test/psModules/src/camera/Makefile.am
===================================================================
--- branches/meh_branches/ppstack_test/psModules/src/camera/Makefile.am	(revision 33596)
+++ branches/meh_branches/ppstack_test/psModules/src/camera/Makefile.am	(revision 34041)
@@ -30,5 +30,6 @@
 	pmReadoutStack.c \
 	pmReadoutFake.c \
-	pmFPABin.c
+	pmFPABin.c \
+	pmFPAExpNumIO.c
 
 pkginclude_HEADERS = \
@@ -59,5 +60,6 @@
 	pmReadoutStack.h \
 	pmReadoutFake.h \
-	pmFPABin.h
+	pmFPABin.h \
+	pmFPAExpNumIO.h
 
 CLEANFILES = *~
Index: branches/meh_branches/ppstack_test/psModules/src/camera/pmFPAExpNumIO.c
===================================================================
--- branches/meh_branches/ppstack_test/psModules/src/camera/pmFPAExpNumIO.c	(revision 34041)
+++ branches/meh_branches/ppstack_test/psModules/src/camera/pmFPAExpNumIO.c	(revision 34041)
@@ -0,0 +1,450 @@
+#include <stdio.h>
+#include <pslib.h>
+#include <string.h>
+
+#include "pmHDU.h"
+#include "pmHDUUtils.h"
+#include "pmFPA.h"
+#include "pmFPAview.h"
+#include "pmFPAfile.h"
+#include "pmFPAfileFitsIO.h"
+#include "pmFPAHeader.h"
+#include "pmConceptsRead.h"
+#include "pmConceptsWrite.h"
+
+#include "pmFPAExpNumIO.h"
+
+//////////////////////////////////////////////////////////////////////////////////////////////////////////////
+
+
+static bool pmReadoutReadExpNum(pmReadout *ro, psFits *fits)
+{
+    PM_ASSERT_READOUT_NON_NULL(ro, false);
+    PS_ASSERT_FITS_NON_NULL(fits, false);
+
+    pmCell *cell = ro->parent;          // Cell of interest
+    if (!cell) {
+        psError(PS_ERR_UNEXPECTED_NULL, true, "No cell associated with readout.");
+        return false;
+    }
+    pmChip *chip = cell->parent;    // Chip of interest
+    pmFPA *fpa = chip->parent;      // FPA of interest
+    pmHDU *hdu = pmHDUGetLowest(fpa, chip, cell); // HDU for readout
+    if (!hdu) {
+        psError(PS_ERR_UNEXPECTED_NULL, true, "No HDU associated with readout.");
+        return false;
+    }
+
+#ifdef notdef
+    if (!psFitsMoveExtName(fits, hdu->extname)) {
+        psError(PS_ERR_IO, false, "Unable to move to EXPNUM image.");
+        return false;
+    }
+#endif
+    psMetadata *header = psFitsReadHeader(NULL, fits); // Header
+    if (!header) {
+        psError(PS_ERR_UNEXPECTED_NULL, false, "Unable to read header for EXPNUM image.");
+        return false;
+    }
+
+    bool data = false;                  // Did we find any data?
+
+    int naxis = psMetadataLookupS32(NULL, header, "NAXIS"); // Number of axes
+    if (naxis > 0) {
+        psImage *expnum = psFitsReadImage(fits, psRegionSet(0, 0, 0, 0), 0);
+        psMetadataAddImage(ro->analysis, PS_LIST_TAIL, "EXPNUM", PS_META_REPLACE,
+                           "EXPNUM", expnum);
+        psFree(expnum);
+        data = true;
+    }
+
+    if (data) {
+        ro->data_exists = true;
+        ro->parent->data_exists = true;
+        ro->parent->parent->data_exists = true;
+    }
+
+    psFree(header);
+
+    return true;
+}
+
+static bool pmCellReadExpNum(pmCell *cell, const pmFPAview *view,
+                              pmFPAfile *file, pmConfig *config)
+{
+    PS_ASSERT_PTR_NON_NULL(view, false);
+    PS_ASSERT_PTR_NON_NULL(cell, false);
+    PS_ASSERT_PTR_NON_NULL(cell->readouts, false);
+
+    pmFPAview *thisView = pmFPAviewAlloc(view->nRows); // Copy of input view
+    *thisView = *view;
+
+    // Create a readout if none exists
+    if (!cell->readouts || cell->readouts->n == 0) {
+        pmReadout *readout = pmReadoutAlloc(cell); // New readout
+        psFree(readout);                // Drop reference
+    }
+
+    cell->data_exists = false;
+    for (int i = 0; i < cell->readouts->n; i++) {
+        pmReadout *readout = cell->readouts->data[i];
+        thisView->readout = i;
+        if (!pmReadoutReadExpNum(readout, file->fits)) {
+            psError(PS_ERR_IO, false, "Unable to read pattern correction.");
+            return false;
+        }
+    }
+    psFree(thisView);
+
+    if (!pmCellReadHeader(cell, file->fits, config)) {
+        psError(PS_ERR_IO, false, "Unable to read header for cell.");
+        return false;
+    }
+    // load in the concept information for this cell
+    if (!pmConceptsReadCell(cell, PM_CONCEPT_SOURCE_HEADER, true, NULL)) {
+        psErrorClear();
+        psWarning("Difficulty reading concepts for cell; attempting to proceed.");
+    }
+
+    return true;
+}
+
+static bool pmChipReadExpNum(pmChip *chip, const pmFPAview *view,
+                              pmFPAfile *file, pmConfig *config)
+{
+    PS_ASSERT_PTR_NON_NULL(view, false);
+    PS_ASSERT_PTR_NON_NULL(chip, false);
+    PS_ASSERT_PTR_NON_NULL(chip->cells, false);
+
+    pmFPAview *thisView = pmFPAviewAlloc(view->nRows); // Copy of input view
+    *thisView = *view;
+
+    chip->data_exists = false;
+    for (int i = 0; i < chip->cells->n; i++) {
+        pmCell *cell = chip->cells->data[i];
+        thisView->cell = i;
+        pmCellReadExpNum(cell, thisView, file, config);
+        if (!cell->data_exists) {
+            continue;
+        }
+        chip->data_exists = true;
+    }
+    psFree(thisView);
+
+    if (!pmChipReadHeader(chip, file->fits, config)) {
+        psError(PS_ERR_IO, false, "Unable to read header for cell.");
+        return false;
+    }
+    if (!pmConceptsReadChip(chip, PM_CONCEPT_SOURCE_HEADER, true, true, NULL)) {
+        psError(PS_ERR_IO, false, "Failed to read concepts for chip.\n");
+        return false;
+    }
+
+    return true;
+}
+
+static bool pmFPAReadExpNum(pmFPA *fpa, const pmFPAview *view,
+                             pmFPAfile *file, pmConfig *config)
+{
+    PS_ASSERT_PTR_NON_NULL(view, false);
+    PS_ASSERT_PTR_NON_NULL(file, false);
+    PS_ASSERT_PTR_NON_NULL(file->fpa, false);
+    PS_ASSERT_PTR_NON_NULL(file->fpa->chips, false);
+
+    pmFPAview *thisView = pmFPAviewAlloc(view->nRows); // Copy of input view
+    *thisView = *view;
+
+    for (int i = 0; i < fpa->chips->n; i++) {
+        pmChip *chip = fpa->chips->data[i];
+        thisView->chip = i;
+        pmChipReadExpNum(chip, thisView, file, config);
+    }
+    psFree(thisView);
+
+    if (!pmFPAReadHeader(fpa, file->fits, config)) {
+        psError(PS_ERR_IO, false, "Unable to read header for cell.");
+        return false;
+    }
+    if (!pmConceptsReadFPA(fpa, PM_CONCEPT_SOURCE_HEADER, true, NULL)) {
+        psError(PS_ERR_IO, false, "Failed to read concepts for fpa.\n");
+        return false;
+    }
+
+    return true;
+}
+
+//////////////////////////////////////////////////////////////////////////////////////////////////////////////
+
+bool pmExpNumRead(const pmFPAview *view, pmFPAfile *file, pmConfig *config)
+{
+    PS_ASSERT_PTR_NON_NULL(view, false);
+    PS_ASSERT_PTR_NON_NULL(file, false);
+    PS_ASSERT_PTR_NON_NULL(file->fpa, false);
+
+    pmFPA *fpa = file->fpa;
+
+    if (view->chip == -1) {
+        pmFPAReadExpNum(fpa, view, file, config);
+        return true;
+    }
+
+    if (view->chip >= fpa->chips->n) {
+        return false;
+    }
+    pmChip *chip = fpa->chips->data[view->chip];
+
+    if (view->cell == -1) {
+        pmChipReadExpNum(chip, view, file, config);
+        return true;
+    }
+
+    if (view->cell >= chip->cells->n) {
+        return false;
+    }
+    pmCell *cell = chip->cells->data[view->cell];
+
+    if (view->readout == -1) {
+        pmCellReadExpNum(cell, view, file, config);
+        return true;
+    }
+
+    if (view->readout >= cell->readouts->n) {
+        return false;
+    }
+    pmReadout *readout = cell->readouts->data[view->readout];
+
+    return pmReadoutReadExpNum(readout, file->fits);
+}
+#ifdef notyet
+// Currently only ppStack writes an EXPNUM image and it uses file type MASK
+bool pmReadoutWriteExpNum(pmReadout *ro, psFits *fits)
+{
+    PM_ASSERT_READOUT_NON_NULL(ro, false);
+    PS_ASSERT_FITS_NON_NULL(fits, false);
+
+    bool gotRow, gotCell;
+    psImage *rowCorr = psMetadataLookupPtr(&gotRow, ro->analysis, PM_PATTERN_ROW_CORRECTION); // Row correction
+    float cellCorr = psMetadataLookupF32(&gotCell, ro->analysis, PM_PATTERN_CELL_CORRECTION); // Cell corr.
+
+    pmCell *cell = ro->parent;          // Cell of interest
+    if (!cell) {
+        psError(PS_ERR_UNEXPECTED_NULL, true, "No cell associated with readout.");
+        return false;
+    }
+    pmChip *chip = cell->parent;    // Chip of interest
+    pmFPA *fpa = chip->parent;      // FPA of interest
+    pmHDU *hdu = pmHDUGetLowest(fpa, chip, cell); // HDU for readout
+    if (!hdu) {
+        psError(PS_ERR_UNEXPECTED_NULL, true, "No HDU associated with readout.");
+        return false;
+    }
+
+    psMetadata *header = psMetadataCopy(NULL, hdu->header); // Header for output
+
+    if (gotCell) {
+        psMetadataAddF32(header, PS_LIST_TAIL, PM_PATTERN_CELL_CORRECTION, PS_META_REPLACE,
+                         "ExpNum cell correction value", cellCorr);
+    }
+
+    if (gotRow) {
+        if (!psFitsWriteImage(fits, header, rowCorr, 0, hdu->extname)) {
+            psError(PS_ERR_IO, false, "Unable to write pattern row correction.");
+            psFree(header);
+            return false;
+        }
+    } else {
+        if (!psFitsWriteBlank(fits, header, hdu->extname)) {
+            psError(PS_ERR_IO, false, "Unable to write pattern cell correction.");
+            psFree(header);
+            return false;
+        }
+    }
+
+    psFree(header);
+    return true;
+}
+
+static bool pmCellWriteExpNum(pmCell *cell, const pmFPAview *view,
+                               pmFPAfile *file, pmConfig *config)
+{
+    PS_ASSERT_PTR_NON_NULL(cell, false);
+    PS_ASSERT_PTR_NON_NULL(cell->readouts, false);
+    PS_ASSERT_PTR_NON_NULL(view, false);
+
+    if (!pmConceptsWriteCell(cell, true, config)) {
+        psError(PS_ERR_IO, false, "Unable to write concepts for cell.");
+        return false;
+    }
+
+    // Only do the FIRST readout --- don't want to write lots of headers
+    pmReadout *readout = cell->readouts->data[0];
+    if (!pmReadoutWriteExpNum(readout, file->fits)) {
+        psError(PS_ERR_IO, false, "Failed to write readout");
+        return false;
+    }
+    return true;
+}
+
+static bool pmChipWriteExpNum(pmChip *chip, const pmFPAview *view,
+                               pmFPAfile *file, pmConfig *config)
+{
+    PS_ASSERT_PTR_NON_NULL(chip, false);
+    PS_ASSERT_PTR_NON_NULL(chip->cells, false);
+    PS_ASSERT_PTR_NON_NULL(view, false);
+
+    if (!pmConceptsWriteChip(chip, true, true, config)) {
+        psError(PS_ERR_IO, false, "Unable to write concepts for chip.\n");
+        return false;
+    }
+
+    pmFPAview *thisView = pmFPAviewAlloc(view->nRows); // Copy of input view
+    *thisView = *view;
+
+    for (int i = 0; i < chip->cells->n; i++) {
+        pmCell *cell = chip->cells->data[i];
+        thisView->cell = i;
+        if (!pmCellWriteExpNum(cell, thisView, file, config)) {
+            psError(PS_ERR_IO, false, "Failed to write %dth cell", i);
+            psFree(thisView);
+            return false;
+        }
+    }
+    psFree(thisView);
+    return true;
+}
+
+static bool pmFPAWriteExpNum(pmFPA *fpa, const pmFPAview *view,
+                              pmFPAfile *file, pmConfig *config)
+{
+    PS_ASSERT_PTR_NON_NULL(view, false);
+    PS_ASSERT_PTR_NON_NULL(fpa, false);
+    PS_ASSERT_PTR_NON_NULL(fpa->chips, false);
+
+    if (!pmConceptsWriteFPA(fpa, true, config)) {
+        psError(PS_ERR_IO, false, "Unable to write concepts for FPA.\n");
+        return false;
+    }
+
+    pmFPAview *thisView = pmFPAviewAlloc(view->nRows); // Copy of input view
+    *thisView = *view;
+
+    for (int i = 0; i < fpa->chips->n; i++) {
+        pmChip *chip = fpa->chips->data[i];
+        thisView->chip = i;
+        if (!pmChipWriteExpNum(chip, thisView, file, config)) {
+            psError(PS_ERR_IO, false, "Failed to write %dth chip", i);
+            psFree(thisView);
+            return false;
+        }
+    }
+    psFree(thisView);
+    return true;
+}
+//////////////////////////////////////////////////////////////////////////////////////////////////////////////
+
+// Currently only ppStack writes an EXPNUM image and it uses file type MASK
+bool pmExpNumWrite(const pmFPAview *view, pmFPAfile *file, pmConfig *config)
+{
+    PS_ASSERT_PTR_NON_NULL(view, false);
+    PS_ASSERT_PTR_NON_NULL(file, false);
+    PS_ASSERT_PTR_NON_NULL(file->fpa, false);
+
+    pmFPA *fpa = pmFPAfileSuitableFPA(file, view, config, false); // Suitable FPA for writing
+
+    if (view->chip == -1) {
+        if (!pmFPAWriteExpNum(fpa, view, file, config)) {
+            psError(PS_ERR_IO, false, "Failed to write pattern correction from fpa");
+            psFree(fpa);
+            return false;
+        }
+        psFree(fpa);
+        return true;
+    }
+
+    if (view->chip >= fpa->chips->n) {
+        psError(PS_ERR_UNKNOWN, false, "Writing chip == %d (>= chips->n == %ld)", view->chip, fpa->chips->n);
+        psFree(fpa);
+        return false;
+    }
+    pmChip *chip = fpa->chips->data[view->chip];
+
+    if (view->cell == -1) {
+        if (!pmChipWriteExpNum(chip, view, file, config)) {
+            psError(PS_ERR_IO, false, "Failed to write pattern correction from chip");
+            psFree(fpa);
+            return false;
+        }
+        psFree(fpa);
+        return true;
+    }
+
+    if (view->cell >= chip->cells->n) {
+        psError(PS_ERR_UNKNOWN, false, "Writing cell == %d (>= cells->n == %ld)",
+                view->cell, chip->cells->n);
+        psFree(fpa);
+        return false;
+    }
+    pmCell *cell = chip->cells->data[view->cell];
+
+    if (view->readout == -1) {
+        if (!pmCellWriteExpNum(cell, view, file, config)) {
+            psError(PS_ERR_IO, false, "Failed to write pattern correction from cell");
+            psFree(fpa);
+            return false;
+        }
+        psFree(fpa);
+        return true;
+    }
+
+    if (view->readout >= cell->readouts->n) {
+        psError(PS_ERR_UNKNOWN, false, "Writing readout == %d (>= readouts->n == %ld)",
+                view->readout, cell->readouts->n);
+        psFree(fpa);
+        return false;
+    }
+    pmReadout *readout = cell->readouts->data[view->readout];
+
+    if (!pmReadoutWriteExpNum(readout, file->fits)) {
+        psError(PS_ERR_IO, false, "Failed to write pattern correction from readout %d", view->readout);
+        psFree(fpa);
+        return false;
+    }
+
+    psFree(fpa);
+    return true;
+}
+
+bool pmExpNumWritePHU(const pmFPAview *view, pmFPAfile *file, pmConfig *config)
+{
+    PS_ASSERT_PTR_NON_NULL(view, false);
+    PS_ASSERT_PTR_NON_NULL(file, false);
+
+    if (file->wrote_phu) {
+        return true;
+    }
+
+    // find the FPA phu
+    pmFPA *fpa = pmFPAfileSuitableFPA(file, view, config, false); // Suitable FPA for writing
+    pmHDU *phu = psMemIncrRefCounter(pmFPAviewThisPHU(view, fpa));
+    psFree(fpa);
+
+    // if there is no PHU, this is a single header+image (extension-less) file. This could be the case for an
+    // input SPLIT set of files being written out as a MEF.  if there is a PHU, write it out as a 'blank'
+    psMetadata *outhead = psMetadataAlloc();
+    if (phu) {
+        psMetadataCopy (outhead, phu->header);
+    }
+    psFree(phu);
+
+    pmConfigConformHeader(outhead, file->format);
+
+    psFitsWriteBlank(file->fits, outhead, "");
+    file->wrote_phu = true;
+
+    psTrace("pmFPAfile", 5, "wrote phu %s (type: %d)\n", file->filename, file->type);
+    psFree(outhead);
+
+    return true;
+}
+#endif // notyet
+
Index: branches/meh_branches/ppstack_test/psModules/src/camera/pmFPAExpNumIO.h
===================================================================
--- branches/meh_branches/ppstack_test/psModules/src/camera/pmFPAExpNumIO.h	(revision 34041)
+++ branches/meh_branches/ppstack_test/psModules/src/camera/pmFPAExpNumIO.h	(revision 34041)
@@ -0,0 +1,36 @@
+#ifndef PM_EXPNUM_IO_H
+#define PM_EXPNUM_IO_H
+
+#include <pslib.h>
+
+#include <pmHDU.h>
+#include <pmFPA.h>
+
+
+bool pmExpNumRead(const pmFPAview *view, ///< View into which to read
+                   pmFPAfile *file, ///< File from which to read
+                   pmConfig *config ///< Configuration
+    );
+
+#ifdef notyet
+// Currently we don't need write functions for the EXPNUM images. 
+// Only ppStack creates them and it uses a MASK file type
+
+/// Write pattern correction within a readout to a FITS file
+bool pmReadoutWritePattern(
+    pmReadout *readout,                 ///< Readout for which to write pattern correction (in analysis MD)
+    psFits *fits                        ///< FITS file to which to write
+    );
+bool pmPatternWrite(const pmFPAview *view, ///< View from which to write
+                    pmFPAfile *file, ///< File to which to write
+                    pmConfig *config ///< Configuration
+    );
+
+bool pmPatternWritePHU(const pmFPAview *view, // View to PHU
+                       pmFPAfile *file, ///< File to which to write
+                       pmConfig *config ///< Configuration
+    );
+
+#endif // notyet
+
+#endif
Index: branches/meh_branches/ppstack_test/psModules/src/camera/pmFPAfile.c
===================================================================
--- branches/meh_branches/ppstack_test/psModules/src/camera/pmFPAfile.c	(revision 33596)
+++ branches/meh_branches/ppstack_test/psModules/src/camera/pmFPAfile.c	(revision 34041)
@@ -382,5 +382,5 @@
         // Number of the file in list
         psString num = NULL;            // Number to use
-        psStringAppend(&num, "%03d", file->fileID);
+        psStringAppend(&num, "%" PRId64, file->fileID);
         psStringSubstitute(&newRule, num, "{FILE.ID}");
         psFree(num);
@@ -543,4 +543,7 @@
     if (!strcasecmp(type, "PATTERN")) {
         return PM_FPA_FILE_PATTERN;
+    }
+    if (!strcasecmp(type, "EXPNUM")) {
+        return PM_FPA_FILE_EXPNUM;
     }
 
@@ -589,4 +592,6 @@
       case PM_FPA_FILE_PATTERN:
         return "PATTERN";
+      case PM_FPA_FILE_EXPNUM:
+        return "EXPNUM";
       default:
         return ("NONE");
Index: branches/meh_branches/ppstack_test/psModules/src/camera/pmFPAfile.h
===================================================================
--- branches/meh_branches/ppstack_test/psModules/src/camera/pmFPAfile.h	(revision 33596)
+++ branches/meh_branches/ppstack_test/psModules/src/camera/pmFPAfile.h	(revision 34041)
@@ -51,4 +51,5 @@
     PM_FPA_FILE_PATTERN,
     PM_FPA_FILE_LINEARITY,
+    PM_FPA_FILE_EXPNUM,
 } pmFPAfileType;
 
@@ -113,5 +114,5 @@
 
     int fileIndex;			// Index of file
-    int fileID;				// internal sequence number
+    psS64 fileID;		        // internal sequence number
 
     psS64 imageId, sourceId;            // Image and source identifiers
Index: branches/meh_branches/ppstack_test/psModules/src/camera/pmFPAfileIO.c
===================================================================
--- branches/meh_branches/ppstack_test/psModules/src/camera/pmFPAfileIO.c	(revision 33596)
+++ branches/meh_branches/ppstack_test/psModules/src/camera/pmFPAfileIO.c	(revision 34041)
@@ -53,4 +53,5 @@
 #include "pmSubtractionIO.h"
 #include "pmPatternIO.h"
+#include "pmFPAExpNumIO.h"
 #include "pmConcepts.h"
 #include "pmConfigRun.h"
@@ -231,4 +232,7 @@
         status = pmAstromModelReadForView (view, file, config);
         break;
+      case PM_FPA_FILE_EXPNUM:
+        status = pmExpNumRead(view, file, config);
+        break;
       case PM_FPA_FILE_ASTROM_REFSTARS:
       case PM_FPA_FILE_JPEG:
@@ -291,4 +295,5 @@
       case PM_FPA_FILE_DARK:
       case PM_FPA_FILE_PATTERN:
+      case PM_FPA_FILE_EXPNUM:
         {
             // create FPA structure component based on view
@@ -500,4 +505,9 @@
         status = pmFPAviewWriteSourcePlot (view, file, config);
         break;
+
+      case PM_FPA_FILE_EXPNUM:
+        // when ppStack output's EXPNUM file it uses a file rule where the file type is MASK
+        psError(PS_ERR_IO, true, "cannot write type EXPNUM (%s)", file->name);
+        return false;
 
       case PM_FPA_FILE_WCS:
@@ -562,4 +572,5 @@
       case PM_FPA_FILE_ASTROM_REFSTARS:
       case PM_FPA_FILE_LINEARITY:
+      case PM_FPA_FILE_EXPNUM:
         psTrace ("psModules.camera", 5, "closing %s (%s) (%d:%d:%d)\n", file->filename, file->name, view->chip, view->cell, view->readout);
         status = psFitsClose (file->fits);
@@ -636,4 +647,5 @@
       case PM_FPA_FILE_ASTROM_MODEL:
       case PM_FPA_FILE_ASTROM_REFSTARS:
+      case PM_FPA_FILE_EXPNUM:
         psTrace ("psModules.camera", 6, "NOT freeing %s (%s) : save for further analysis\n", file->filename, file->name);
         return true;
@@ -797,4 +809,5 @@
       case PM_FPA_FILE_ASTROM_REFSTARS:
       case PM_FPA_FILE_LINEARITY:
+      case PM_FPA_FILE_EXPNUM:
         psTrace ("psModules.camera", 5, "opening %s (%s) (%d:%d:%d)\n",
                  file->filename, file->name, view->chip, view->cell, view->readout);
@@ -996,4 +1009,5 @@
         status = pmAstromRefstarsWritePHU (view, file, config);
         break;
+      case PM_FPA_FILE_EXPNUM:
       case PM_FPA_FILE_ASTROM_MODEL:
       case PM_FPA_FILE_SX:
Index: branches/meh_branches/ppstack_test/psModules/src/camera/pmReadoutFake.c
===================================================================
--- branches/meh_branches/ppstack_test/psModules/src/camera/pmReadoutFake.c	(revision 33596)
+++ branches/meh_branches/ppstack_test/psModules/src/camera/pmReadoutFake.c	(revision 34041)
@@ -162,4 +162,9 @@
             fakeRadius = PS_MAX(fakeRadius, fakeModel->modelRadius(fakeModel->params, minFlux));
         }
+	//MEH - add some outputs to test --
+        psTrace("psModules.camera", 10, "TESTING:: source at %f,%f: radius %f fakeradius %f minflux %f psf params later \n",
+                fakeModel->params->data.F32[PM_PAR_XPOS],fakeModel->params->data.F32[PM_PAR_YPOS],radius,fakeRadius,minFlux);
+	//MEH -- seems radius is passed as float 0.0.. hardcode test - 
+
         if (radius > 0) {
             fakeRadius = PS_MAX(fakeRadius, radius);
@@ -210,9 +215,12 @@
     const pmPSF *psf = args->data[7];         // PSF
     float minFlux = PS_SCALAR_VALUE(args->data[8], F32); // Minimum flux
-    float radius = PS_SCALAR_VALUE(args->data[9], F32);  // Minimum radius
+    float radius = PS_SCALAR_VALUE(args->data[9], S32);  // Minimum radius - typecast to float from S32 outside of PS_SCALAR_VALUE otherwise sets 0.0 
     bool circularise = PS_SCALAR_VALUE(args->data[10], U8); // Circularise PSF?
     bool normalisePeak = PS_SCALAR_VALUE(args->data[11], U8); // Normalise for peak?
     int groupIndex = PS_SCALAR_VALUE(args->data[12], S32); // Group index
     int cellIndex = PS_SCALAR_VALUE(args->data[13], S32);  // Cell index
+
+    //MEH
+    psTrace("psModules.camera", 10, "TESTING:thread: radius %f minflux %f  \n",radius,minFlux);
 
     return readoutFake(readout, groups, x, y, mag, xOffset, yOffset, psf, minFlux, radius, circularise,
@@ -242,5 +250,4 @@
     return old;
 }
-
 
 bool pmReadoutFakeFromVectors(pmReadout *readout, int numCols, int numRows,
@@ -287,4 +294,5 @@
     }
 
+    psTrace("psModules.camera", 10, "TESTING:prethread: radius %d minflux %f threaded %d \n",radius,minFlux,threaded);
     if (threaded) {
         for (int i = 0; i < groups->groups->n; i++) {
@@ -303,4 +311,5 @@
                 psArrayAdd(args, 1, (pmPSF*)psf);
                 PS_ARRAY_ADD_SCALAR(args, minFlux, PS_TYPE_F32);
+		//MEH - S32 is proper type here
                 PS_ARRAY_ADD_SCALAR(args, radius, PS_TYPE_S32);
                 PS_ARRAY_ADD_SCALAR(args, circularise, PS_TYPE_U8);
@@ -349,5 +358,4 @@
 }
 
-
 bool pmReadoutFakeFromSources(pmReadout *readout, int numCols, int numRows, const psArray *sources,
                               pmSourceMode sourceMask, const psVector *xOffset, const psVector *yOffset,
@@ -392,4 +400,7 @@
     mag->n = numGood;
 
+    //MEH -- convert int->float for calling downstream since nothing calls Vectors directly? psphot does.. argh!
+    psTrace("psModules.camera", 10, "TESTING:: radius %d minflux %f  \n",radius,minFlux);
+
     bool status = pmReadoutFakeFromVectors(readout, numCols, numRows, x, y, mag, xOffset, yOffset, psf,
                                            minFlux, radius, circularise, normalisePeak);
Index: branches/meh_branches/ppstack_test/psModules/src/detrend/pmDetrendDB.c
===================================================================
--- branches/meh_branches/ppstack_test/psModules/src/detrend/pmDetrendDB.c	(revision 33596)
+++ branches/meh_branches/ppstack_test/psModules/src/detrend/pmDetrendDB.c	(revision 34041)
@@ -107,4 +107,6 @@
         DETREND_STRING_CASE(ASTROM);
         DETREND_STRING_CASE(NOISEMAP);
+	DETREND_STRING_CASE(VIDEOMASK);
+	DETREND_STRING_CASE(VIDEODARK);
 	DETREND_STRING_CASE(LINEARITY);
     default:
Index: branches/meh_branches/ppstack_test/psModules/src/detrend/pmDetrendDB.h
===================================================================
--- branches/meh_branches/ppstack_test/psModules/src/detrend/pmDetrendDB.h	(revision 33596)
+++ branches/meh_branches/ppstack_test/psModules/src/detrend/pmDetrendDB.h	(revision 34041)
@@ -36,4 +36,6 @@
     PM_DETREND_TYPE_ASTROM,
     PM_DETREND_TYPE_NOISEMAP,
+    PM_DETREND_TYPE_VIDEOMASK,
+    PM_DETREND_TYPE_VIDEODARK,
     PM_DETREND_TYPE_LINEARITY,
 } pmDetrendType;
Index: branches/meh_branches/ppstack_test/psModules/src/extras/psVectorBracket.c
===================================================================
--- branches/meh_branches/ppstack_test/psModules/src/extras/psVectorBracket.c	(revision 33596)
+++ branches/meh_branches/ppstack_test/psModules/src/extras/psVectorBracket.c	(revision 34041)
@@ -51,6 +51,6 @@
         }
     }
-    // at this point, index[Nhi] >= key > index[Nlo]
-    N = Nhi;
+    N = (Nhi >= index->n) ? Nhi - 1 : Nhi;
+    // at this point, index[N] >= key > index[Nlo]
     while ((index->data.F32[N] >= key) && (N > Nlo)) {
         N--;
Index: branches/meh_branches/ppstack_test/psModules/src/imcombine/pmSubtraction.c
===================================================================
--- branches/meh_branches/ppstack_test/psModules/src/imcombine/pmSubtraction.c	(revision 33596)
+++ branches/meh_branches/ppstack_test/psModules/src/imcombine/pmSubtraction.c	(revision 34041)
@@ -835,4 +835,5 @@
 //#ifdef TESTING
     //MEH fix index issue - index is supposed to be an image index but from where?
+    //MEH - index conflict or changed in past?
     for (int j = 0; j < kernels->num; j++) {
         if (stamp->convolutions1) {
Index: branches/meh_branches/ppstack_test/psModules/src/imcombine/pmSubtractionStamps.c
===================================================================
--- branches/meh_branches/ppstack_test/psModules/src/imcombine/pmSubtractionStamps.c	(revision 33596)
+++ branches/meh_branches/ppstack_test/psModules/src/imcombine/pmSubtractionStamps.c	(revision 34041)
@@ -1232,5 +1232,13 @@
 	// XXX this is somewhat arbitrary...
 	if (source->psfMagErr > 0.05) continue;
-	if (fabs(source->psfMag - source->apMag) > 0.5) continue;
+        if (isfinite(source->apMag)) {
+            if (fabs(source->psfMag - source->apMag) > 0.5) continue;
+        } else if (isfinite(source->apMagRaw)) {
+            if (fabs(source->psfMag - source->apMagRaw) > 0.5) continue;
+        } else {
+            // XXX: Should we carry on or drop this source?
+            // drop it for now
+            continue;
+        }
 
         if (source->modelPSF) {
Index: branches/meh_branches/ppstack_test/psModules/src/objects/models/pmModel_QGAUSS.c
===================================================================
--- branches/meh_branches/ppstack_test/psModules/src/objects/models/pmModel_QGAUSS.c	(revision 33596)
+++ branches/meh_branches/ppstack_test/psModules/src/objects/models/pmModel_QGAUSS.c	(revision 34041)
@@ -402,4 +402,9 @@
     assert (psf->params->n > PM_PAR_YPOS);
     assert (psf->params->n > PM_PAR_XPOS);
+
+    if (! isfinite(Io)) {
+        fprintf(stderr, "non-finite Io passed to PM_MODEL_PARAMS_FROM_PSF\n");
+        return false;
+    }
 
     PAR[PM_PAR_SKY]  = 0.0;
Index: branches/meh_branches/ppstack_test/psModules/src/objects/models/pmModel_SERSIC.c
===================================================================
--- branches/meh_branches/ppstack_test/psModules/src/objects/models/pmModel_SERSIC.c	(revision 33596)
+++ branches/meh_branches/ppstack_test/psModules/src/objects/models/pmModel_SERSIC.c	(revision 34041)
@@ -192,4 +192,10 @@
     psF32 z0 = PAR[PM_PAR_I0]*f1;
     psF32 f0 = PAR[PM_PAR_SKY] + z0;
+
+    if (!isfinite(z0)) {
+        fprintf(stderr, "z0 is not finite for %f %f %f %f %f.  Parameters: \n", X, Y, radius, z, f1);
+        fprintf(stderr, "%f %f %f %f %f %f %f %f\n", PAR[0], PAR[1], PAR[2], PAR[3], PAR[4],
+            PAR[5], PAR[6], PAR[7]);
+    }
 
     assert (isfinite(f2));
Index: branches/meh_branches/ppstack_test/psModules/src/objects/pmFootprintCullPeaks.c
===================================================================
--- branches/meh_branches/ppstack_test/psModules/src/objects/pmFootprintCullPeaks.c	(revision 33596)
+++ branches/meh_branches/ppstack_test/psModules/src/objects/pmFootprintCullPeaks.c	(revision 34041)
@@ -178,9 +178,9 @@
 	    psArray *myFP = pmFootprintsFind(subImg, threshold, 5);
 	    if (!myFP) {
-		psWarning ("missing footprint?");
+		psWarning ("missing footprint? threshold: %.f", threshold);
 		continue;
 	    }
 	    if (!myFP->n) {
-		psWarning ("empty footprint?");
+		psWarning ("empty footprint? threshold: %.f", threshold);
 		psFree (myFP);
 		continue;
Index: branches/meh_branches/ppstack_test/psModules/src/objects/pmModelFuncs.h
===================================================================
--- branches/meh_branches/ppstack_test/psModules/src/objects/pmModelFuncs.h	(revision 33596)
+++ branches/meh_branches/ppstack_test/psModules/src/objects/pmModelFuncs.h	(revision 34041)
@@ -29,4 +29,6 @@
 # define PM_MODEL_FUNCS_H
 
+# define HAVE_MODEL_VAR 1
+
 /// @addtogroup Objects Object Detection / Analysis Functions
 /// @{
@@ -46,13 +48,16 @@
 
 typedef enum {
-    PM_MODEL_OP_NONE    = 0x00,
-    PM_MODEL_OP_FUNC    = 0x01,
-    PM_MODEL_OP_RES0    = 0x02,
-    PM_MODEL_OP_RES1    = 0x04,
-    PM_MODEL_OP_FULL    = 0x07,
-    PM_MODEL_OP_SKY     = 0x08,
-    PM_MODEL_OP_CENTER  = 0x10,
-    PM_MODEL_OP_NORM    = 0x20,
-    PM_MODEL_OP_NOISE   = 0x40,
+    PM_MODEL_OP_NONE     = 0x00,
+    PM_MODEL_OP_FUNC     = 0x01,
+    PM_MODEL_OP_RES0     = 0x02,
+    PM_MODEL_OP_RES1     = 0x04,
+    PM_MODEL_OP_FULL     = 0x07,
+    PM_MODEL_OP_SKY      = 0x08,
+    PM_MODEL_OP_CENTER   = 0x10,
+    PM_MODEL_OP_NORM     = 0x20,
+    PM_MODEL_OP_NOISE    = 0x40,
+# if (HAVE_MODEL_VAR)
+    PM_MODEL_OP_MODELVAR = 0x80,
+# endif
 } pmModelOpMode;
 
Index: branches/meh_branches/ppstack_test/psModules/src/objects/pmModelUtils.c
===================================================================
--- branches/meh_branches/ppstack_test/psModules/src/objects/pmModelUtils.c	(revision 33596)
+++ branches/meh_branches/ppstack_test/psModules/src/objects/pmModelUtils.c	(revision 34041)
@@ -144,8 +144,16 @@
 
     *Io = source->peak->rawFlux;
+
+#ifndef ALLOW_NONFINITE_PEAK
+    // Gene says fail of peak !finite
+    if (!isfinite(*Io)) return false;
+#else 
+    // This is the way it used to be. Somtimes an infinite value Io made it's way down the pipeline
+    // causing assertion failures
     if (!isfinite(*Io) && !source->moments) return false;
 
     *Io = source->moments->Peak;
     if (!isfinite(*Io)) return false;
+#endif
 
     return true;
Index: branches/meh_branches/ppstack_test/psModules/src/objects/pmPCMdata.c
===================================================================
--- branches/meh_branches/ppstack_test/psModules/src/objects/pmPCMdata.c	(revision 33596)
+++ branches/meh_branches/ppstack_test/psModules/src/objects/pmPCMdata.c	(revision 34041)
@@ -136,4 +136,18 @@
 	    sum += value;
 	}
+    }
+
+    if (!(sum > 0.0)) {
+        // Crazy PSF image print out some debugging information ...
+        fprintf(stderr, "invalid kernel sum %f found by pmPCMkernelFromPSF\n", sum);    for (int j = psf->yMin; j <= psf->yMax; j++) {
+            fprintf(stderr, "Row %d\n", j);
+            for (int i = psf->xMin; i <= psf->xMax; i++) {
+                double value = source->psfImage->data.F32[y0 + j][x0 + i];
+                fprintf(stderr, "  %d %f\n", i, value);
+            }
+        }
+        fflush(stderr);
+        // ... but avoid the asssertion two lines down by escaping
+        goto escape;
     }
     assert (sum > 0.0);
Index: branches/meh_branches/ppstack_test/psModules/src/objects/pmPSF.c
===================================================================
--- branches/meh_branches/ppstack_test/psModules/src/objects/pmPSF.c	(revision 33596)
+++ branches/meh_branches/ppstack_test/psModules/src/objects/pmPSF.c	(revision 34041)
@@ -455,5 +455,6 @@
         return NAN;
     }
-
+    //MEH -- multiple ways of doing this in past? add param details to log useful for ppStack?
+   
     // get the model full-width at half-max
     float fwhmMajor = 2*model->modelRadius (model->params, 0.5);
Index: branches/meh_branches/ppstack_test/psModules/src/objects/pmPSFtryFitEXT.c
===================================================================
--- branches/meh_branches/ppstack_test/psModules/src/objects/pmPSFtryFitEXT.c	(revision 33596)
+++ branches/meh_branches/ppstack_test/psModules/src/objects/pmPSFtryFitEXT.c	(revision 34041)
@@ -73,4 +73,12 @@
             continue;
         }
+        // If mask object does not exist, mark the source as bad.
+        // We cannot proceed with it because psImageMaskPixels leaves an uncleared error code last which causes
+        // psphot to exit with a fault. 
+        if (source->maskObj == NULL) {
+            psTrace ("psModules.objects", 4, "source %d (%d,%d) : null maskObj\n", i, source->peak->x, source->peak->y);
+            psfTry->mask->data.PS_TYPE_VECTOR_MASK_DATA[i] = PSFTRY_MASK_EXT_FAIL;
+            continue;
+        }
 
         source->modelEXT = pmSourceModelGuess (source, options->type);
@@ -89,5 +97,5 @@
 
         // clear object mask to define valid pixels
-	psImageMaskPixels (source->maskObj, "AND", PS_NOT_IMAGE_MASK(markVal)); // clear the circular mask
+        psImageMaskPixels (source->maskObj, "AND", PS_NOT_IMAGE_MASK(markVal)); // clear the circular mask
 
         // exclude the poor fits
Index: branches/meh_branches/ppstack_test/psModules/src/objects/pmSource.c
===================================================================
--- branches/meh_branches/ppstack_test/psModules/src/objects/pmSource.c	(revision 33596)
+++ branches/meh_branches/ppstack_test/psModules/src/objects/pmSource.c	(revision 34041)
@@ -39,4 +39,5 @@
 #include "pmSourceExtendedPars.h"
 #include "pmSourceDiffStats.h"
+#include "pmSourcePhotometry.h"
 #include "pmSource.h"
 
@@ -50,4 +51,7 @@
     psFree(tmp->pixels);
     psFree(tmp->variance);
+# if (HAVE_MODEL_VAR)
+    psFree(tmp->modelVar);
+# endif
     psFree(tmp->maskObj);
     psFree(tmp->maskView);
@@ -76,4 +80,7 @@
     psFree (source->pixels);
     psFree (source->variance);
+# if (HAVE_MODEL_VAR)
+    psFree (source->modelVar);
+# endif
     psFree (source->maskObj);
     psFree (source->maskView);
@@ -83,4 +90,7 @@
     source->pixels = NULL;
     source->variance = NULL;
+# if (HAVE_MODEL_VAR)
+    source->modelVar = NULL;
+# endif
     source->maskObj = NULL;
     source->maskView = NULL;
@@ -112,4 +122,7 @@
     source->pixels = NULL;
     source->variance = NULL;
+# if (HAVE_MODEL_VAR)
+    source->modelVar = NULL;
+# endif
     source->maskObj = NULL;
     source->maskView = NULL;
@@ -159,5 +172,7 @@
     source->radialAper = NULL;
     source->parent = NULL;
+    source->tmpPtr = NULL;
     source->imageID = -1;
+    source->nFrames = 0;
 
     psMemSetDeallocator(source, (psFreeFunc) sourceFree);
@@ -197,4 +212,7 @@
     source->pixels   = in->pixels   ? psImageCopyView(NULL, in->pixels)   : NULL;
     source->variance = in->variance ? psImageCopyView(NULL, in->variance) : NULL;
+# if (HAVE_MODEL_VAR)
+    source->modelVar = NULL;
+# endif
     source->maskView = in->maskView ? psImageCopyView(NULL, in->maskView) : NULL;
 
@@ -1026,4 +1044,9 @@
     bool addNoise = mode & PM_MODEL_OP_NOISE;
 
+# if (HAVE_MODEL_VAR)
+    bool addModelVar = mode & PM_MODEL_OP_MODELVAR;
+    if (addModelVar) psAssert (source->modelVar, "programming error");
+# endif
+
     // require the use of pmModelAddWithOffset if we are adding noise (because the model size and norm are rescaled)
     if (!addNoise && source->modelFlux) {
@@ -1047,5 +1070,9 @@
         }
 
+# if (HAVE_MODEL_VAR)
+        psF32 **target = addModelVar ? source->modelVar->data.F32 : source->pixels->data.F32;
+# else
         psF32 **target = source->pixels->data.F32;
+# endif
 
         for (int iy = 0; iy < source->modelFlux->numRows; iy++) {
@@ -1062,4 +1089,13 @@
             }
         }
+# if (HAVE_MODEL_VAR)
+	if (!addModelVar) {
+	    if (add) {
+		source->tmpFlags &= ~PM_SOURCE_TMPF_SUBTRACTED;
+	    } else {
+		source->tmpFlags |= PM_SOURCE_TMPF_SUBTRACTED;
+	    }
+	}
+# else
 	if (add) {
 	    source->tmpFlags &= ~PM_SOURCE_TMPF_SUBTRACTED;
@@ -1067,4 +1103,5 @@
 	    source->tmpFlags |= PM_SOURCE_TMPF_SUBTRACTED;
 	}
+# endif
         return true;
     }
@@ -1074,11 +1111,24 @@
         target = source->variance;
     }
+# if (HAVE_MODEL_VAR)
+    if (addModelVar) {
+        target = source->modelVar;
+    }
+# endif
 
     if (add) {
 	status = pmModelAddWithOffset (target, source->maskObj, model, PM_MODEL_OP_FULL, maskVal, dx, dy);
+# if (HAVE_MODEL_VAR)
+	if (!addNoise && !addModelVar) source->tmpFlags &= ~PM_SOURCE_TMPF_SUBTRACTED;
+# else
 	source->tmpFlags &= ~PM_SOURCE_TMPF_SUBTRACTED;
+# endif
     } else {
 	status = pmModelSubWithOffset (target, source->maskObj, model, PM_MODEL_OP_FULL, maskVal, dx, dy);
+# if (HAVE_MODEL_VAR)
+	if (!addNoise && !addModelVar) source->tmpFlags |= PM_SOURCE_TMPF_SUBTRACTED;
+# else
 	source->tmpFlags |= PM_SOURCE_TMPF_SUBTRACTED;
+# endif
     }
 
@@ -1156,4 +1206,92 @@
     PAR[PM_PAR_SYY] = oldshape.sy;
     PAR[PM_PAR_SXY] = oldshape.sxy;
+
+    return true;
+}
+
+bool pmSourceSmoothOp (pmSource *source, pmModelOpMode mode, psImage *target, float sigma, bool add, psImageMaskType maskVal, int dx, int dy)
+{
+//    assert (mode & PM_MODEL_OP_NOISE);
+    PS_ASSERT_PTR_NON_NULL(source, false);
+    PS_ASSERT_PTR_NON_NULL(source->peak, false);
+    PS_ASSERT_PTR_NON_NULL(target, false);
+
+    if (add) {
+        psTrace ("psphot", 3, "adding smoothed object at %f,%f\n", source->peak->xf, source->peak->yf);
+    } else {
+        psTrace ("psphot", 3, "removing smooted object at %f,%f\n", source->peak->xf, source->peak->yf);
+    }
+
+    pmModel *model = pmSourceGetModel(NULL, source);
+    if (!model) {
+        return false;
+    }
+    pmSourceSmoothOpModel (model, source, mode, target, sigma, add, maskVal, dx, dy);
+
+    return true;
+}
+
+bool pmSourceSmoothOpModel (pmModel *model, pmSource *source, pmModelOpMode mode, psImage *target, float sigma, bool add, psImageMaskType maskVal, int dx, int dy) 
+{
+    bool status;
+    psEllipseShape oldshape;
+    psEllipseShape newshape;
+    psEllipseAxes axes;
+
+    if (add) {
+	psTrace ("psphot", 4, "adding smoothed object at %f,%f\n", model->params->data.F32[PM_PAR_XPOS], model->params->data.F32[PM_PAR_YPOS]);
+    } else {
+	psTrace ("psphot", 4, "removing smoothed object at %f,%f\n", model->params->data.F32[PM_PAR_XPOS], model->params->data.F32[PM_PAR_YPOS]);
+    }
+
+    psF32 *PAR = model->params->data.F32;
+
+    // Isn't this hanging around somewhere?
+    float oldFlux = NAN;
+    if (!pmSourcePhotometryModel (NULL, &oldFlux, model)) return false;
+
+    // save original values
+    float oldI0  = PAR[PM_PAR_I0];
+    float oldsxx = PAR[PM_PAR_SXX];
+    float oldsyy = PAR[PM_PAR_SYY];
+    float oldsxy = PAR[PM_PAR_SXY];
+
+    // Since we are going to scale the flux correctly we need to get our
+    // factors of sqrt(2) right
+    oldshape.sx  = oldsxx / M_SQRT2;
+    oldshape.sy  = oldsyy / M_SQRT2;
+    oldshape.sxy = oldsxy;
+
+    // XXX can this be done more intelligently?
+    if (oldI0 == 0.0) return false;
+    if (!isfinite(oldI0)) return false;
+
+    // increase size and height of source
+    axes = psEllipseShapeToAxes (oldshape, 20.0);
+    axes.major = sqrt(PS_SQR(axes.major) + PS_SQR(sigma));
+    axes.minor = sqrt(PS_SQR(axes.minor) + PS_SQR(sigma));
+    newshape = psEllipseAxesToShape (axes);
+    PAR[PM_PAR_SXX] = newshape.sx * M_SQRT2;
+    PAR[PM_PAR_SYY] = newshape.sy * M_SQRT2;
+    PAR[PM_PAR_SXY] = newshape.sxy;
+
+    PAR[PM_PAR_I0]  = 1.0;
+
+    float newFlux;
+    if (!pmSourcePhotometryModel (NULL, &newFlux, model)) return false;
+
+    PAR[PM_PAR_I0]  = oldFlux / newFlux;
+
+    if (add) {
+	status = pmModelAddWithOffset (target, source->maskObj, model, mode, maskVal, dx, dy);
+    } else {
+	status = pmModelSubWithOffset (target, source->maskObj, model, mode, maskVal, dx, dy);
+    }
+
+    // restore original values
+    PAR[PM_PAR_I0]  = oldI0;
+    PAR[PM_PAR_SXX] = oldsxx;
+    PAR[PM_PAR_SYY] = oldsyy;
+    PAR[PM_PAR_SXY] = oldsxy;
 
     return true;
Index: branches/meh_branches/ppstack_test/psModules/src/objects/pmSource.h
===================================================================
--- branches/meh_branches/ppstack_test/psModules/src/objects/pmSource.h	(revision 33596)
+++ branches/meh_branches/ppstack_test/psModules/src/objects/pmSource.h	(revision 34041)
@@ -72,4 +72,7 @@
     psImage *pixels;                    ///< Rectangular region including object pixels.
     psImage *variance;			///< Image variance.
+# if (HAVE_MODEL_VAR)
+    psImage *modelVar;			///< variance based on current models
+# endif
     psImage *maskObj;                   ///< unique mask for this object which marks included pixels associated with objects.
     psImage *maskView;                  ///< view into global image mask for this object region
@@ -117,5 +120,7 @@
     psArray *radialAper;		///< radial flux in circular apertures
     pmSource *parent;			///< reference to the master source from which this is derived
+    psPtr *tmpPtr;                      ///< pointer that may be used to store data in a particular module. e.g. psphotKronIterate.
     int imageID;
+    psU16 nFrames;
 };
 
@@ -301,4 +306,7 @@
 bool pmSourceNoiseOp (pmSource *source, pmModelOpMode mode, float FACTOR, float SIZE, bool add, psImageMaskType maskVal, int dx, int dy);
 
+bool pmSourceSmoothOp (pmSource *source, pmModelOpMode mode, psImage *target, float sigma, bool add, psImageMaskType maskVal, int dx, int dy);
+bool pmSourceSmoothOpModel (pmModel *model, pmSource *source, pmModelOpMode mode, psImage *target, float sigma, bool add, psImageMaskType maskVal, int dx, int dy);
+
 bool pmSourceOp (pmSource *source, pmModelOpMode mode, bool add, psImageMaskType maskVal, int dx, int dy);
 bool pmSourceCacheModel (pmSource *source, psImageMaskType maskVal);
Index: branches/meh_branches/ppstack_test/psModules/src/objects/pmSourceIO.c
===================================================================
--- branches/meh_branches/ppstack_test/psModules/src/objects/pmSourceIO.c	(revision 33596)
+++ branches/meh_branches/ppstack_test/psModules/src/objects/pmSourceIO.c	(revision 34041)
@@ -1069,4 +1069,9 @@
                 for (long i = sources->n -1; i >= 0; i--) {
                     pmSource *source = sources->data[i];
+                    if (source->seq < 0) {
+                        // This can happen cmf files that have been corrupted
+                        psError(PS_ERR_IO, true, "seq < 0 for source %ld: Suspect %s is corrupt", i, file->origname);
+                        return false;
+                    }
                     if (source->seq > seq_max) {
                         seq_max = source->seq;
Index: branches/meh_branches/ppstack_test/psModules/src/objects/pmSourceIO_CMF.c.in
===================================================================
--- branches/meh_branches/ppstack_test/psModules/src/objects/pmSourceIO_CMF.c.in	(revision 33596)
+++ branches/meh_branches/ppstack_test/psModules/src/objects/pmSourceIO_CMF.c.in	(revision 34041)
@@ -80,10 +80,9 @@
     table = psArrayAllocEmpty (sources->n);
 
-    short nImageOverlap; 
     float magOffset; 
     float zeroptErr; 
     float fwhmMajor; 
     float fwhmMinor;
-    pmSourceOutputsCommonValues (&nImageOverlap, &magOffset, &zeroptErr, &fwhmMajor, &fwhmMinor, readout, imageHeader);
+    pmSourceOutputsCommonValues (&magOffset, &zeroptErr, &fwhmMajor, &fwhmMinor, readout, imageHeader);
 
     // we write out PSF-fits for all sources, regardless of quality.  the source flags tell us the state
@@ -179,5 +178,5 @@
 
         // XXX not sure how to get this : need to load Nimages with weight?
-        @ALL@     psMetadataAdd (row, PS_LIST_TAIL, "N_FRAMES",         PS_DATA_U16, "Number of frames overlapping source center", nImageOverlap);
+        @ALL@     psMetadataAdd (row, PS_LIST_TAIL, "N_FRAMES",         PS_DATA_U16, "Number of frames overlapping source center", source->nFrames);
         @ALL@     psMetadataAdd (row, PS_LIST_TAIL, "PADDING",          PS_DATA_S16, "padding", 0);
 
Index: branches/meh_branches/ppstack_test/psModules/src/objects/pmSourceIO_CMF_PS1_DV1.c
===================================================================
--- branches/meh_branches/ppstack_test/psModules/src/objects/pmSourceIO_CMF_PS1_DV1.c	(revision 33596)
+++ branches/meh_branches/ppstack_test/psModules/src/objects/pmSourceIO_CMF_PS1_DV1.c	(revision 34041)
@@ -81,10 +81,9 @@
     }
 
-    short nImageOverlap; 
     float magOffset; 
     float zeroptErr; 
     float fwhmMajor; 
     float fwhmMinor;
-    pmSourceOutputsCommonValues (&nImageOverlap, &magOffset, &zeroptErr, &fwhmMajor, &fwhmMinor, readout, imageHeader);
+    pmSourceOutputsCommonValues (&magOffset, &zeroptErr, &fwhmMajor, &fwhmMinor, readout, imageHeader);
 
     table = psArrayAllocEmpty (sources->n);
@@ -161,6 +160,5 @@
         psMetadataAdd (row, PS_LIST_TAIL, "FLAGS",            PS_DATA_U32, "psphot analysis flags",                      source->mode);
 
-        // XXX not sure how to get this : need to load Nimages with weight?
-        psMetadataAdd (row, PS_LIST_TAIL, "N_FRAMES",         PS_DATA_U16, "Number of frames overlapping source center", nImageOverlap);
+        psMetadataAdd (row, PS_LIST_TAIL, "N_FRAMES",         PS_DATA_U16, "Number of frames overlapping source center", source->nFrames);
         psMetadataAdd (row, PS_LIST_TAIL, "PADDING",          PS_DATA_S16, "padding", 0);
 
Index: branches/meh_branches/ppstack_test/psModules/src/objects/pmSourceIO_CMF_PS1_DV2.c
===================================================================
--- branches/meh_branches/ppstack_test/psModules/src/objects/pmSourceIO_CMF_PS1_DV2.c	(revision 33596)
+++ branches/meh_branches/ppstack_test/psModules/src/objects/pmSourceIO_CMF_PS1_DV2.c	(revision 34041)
@@ -82,10 +82,9 @@
     table = psArrayAllocEmpty (sources->n);
 
-    short nImageOverlap; 
     float magOffset; 
     float zeroptErr; 
     float fwhmMajor; 
     float fwhmMinor;
-    pmSourceOutputsCommonValues (&nImageOverlap, &magOffset, &zeroptErr, &fwhmMajor, &fwhmMinor, readout, imageHeader);
+    pmSourceOutputsCommonValues (&magOffset, &zeroptErr, &fwhmMajor, &fwhmMinor, readout, imageHeader);
 
     // we write out PSF-fits for all sources, regardless of quality.  the source flags tell us the state
@@ -180,5 +179,5 @@
 
         // XXX not sure how to get this : need to load Nimages with weight?
-        psMetadataAdd (row, PS_LIST_TAIL, "N_FRAMES",         PS_DATA_U16, "Number of frames overlapping source center", nImageOverlap);
+        psMetadataAdd (row, PS_LIST_TAIL, "N_FRAMES",         PS_DATA_U16, "Number of frames overlapping source center", source->nFrames);
         psMetadataAdd (row, PS_LIST_TAIL, "PADDING",          PS_DATA_S16, "padding", 0);
 
Index: branches/meh_branches/ppstack_test/psModules/src/objects/pmSourceIO_CMF_PS1_SV1.c
===================================================================
--- branches/meh_branches/ppstack_test/psModules/src/objects/pmSourceIO_CMF_PS1_SV1.c	(revision 33596)
+++ branches/meh_branches/ppstack_test/psModules/src/objects/pmSourceIO_CMF_PS1_SV1.c	(revision 34041)
@@ -84,10 +84,9 @@
     table = psArrayAllocEmpty (sources->n);
 
-    short nImageOverlap; 
     float magOffset; 
     float zeroptErr; 
     float fwhmMajor; 
     float fwhmMinor;
-    pmSourceOutputsCommonValues (&nImageOverlap, &magOffset, &zeroptErr, &fwhmMajor, &fwhmMinor, readout, imageHeader);
+    pmSourceOutputsCommonValues (&magOffset, &zeroptErr, &fwhmMajor, &fwhmMinor, readout, imageHeader);
 
     // we write out PSF-fits for all sources, regardless of quality.  the source flags tell us the state
@@ -170,6 +169,5 @@
         psMetadataAdd (row, PS_LIST_TAIL, "FLAGS2",           PS_DATA_U32, "psphot analysis flags",                     source->mode2);
 
-        // XXX not sure how to get this : need to load Nimages with weight?
-        psMetadataAdd (row, PS_LIST_TAIL, "N_FRAMES",         PS_DATA_U16, "Number of frames overlapping source center", nImageOverlap);
+        psMetadataAdd (row, PS_LIST_TAIL, "N_FRAMES",         PS_DATA_U16, "Number of frames overlapping source center", source->nFrames);
         psMetadataAdd (row, PS_LIST_TAIL, "PADDING",          PS_DATA_S16, "padding", 0);
 
Index: branches/meh_branches/ppstack_test/psModules/src/objects/pmSourceMoments.c
===================================================================
--- branches/meh_branches/ppstack_test/psModules/src/objects/pmSourceMoments.c	(revision 33596)
+++ branches/meh_branches/ppstack_test/psModules/src/objects/pmSourceMoments.c	(revision 34041)
@@ -319,14 +319,14 @@
     }
 
-    source->moments->Mrf = RF/RS;
     source->moments->Mrh = RH/RS;
 
-    // if Mrf (first radial moment) is very small, we are getting into low-significance
+    // if Mrf = RF/RS (first radial moment) is very small, we are getting into low-significance
     // territory.  saturate at minKronRadius.  conversely, if Mrf is >> radius for faint
     // sources, we are clearly making an error.  saturate at radius.
-    float kronRefRadius = MAX(minKronRadius, source->moments->Mrf);
+    float kronRefRadius = MAX(minKronRadius, RF/RS);
     if (source->moments->SN < 10) {
 	kronRefRadius = MIN(radius, kronRefRadius);
     }
+    source->moments->Mrf = kronRefRadius;
 
     // *** now calculate the kron flux values using the 1st radial moment
Index: branches/meh_branches/ppstack_test/psModules/src/objects/pmSourceOutputs.c
===================================================================
--- branches/meh_branches/ppstack_test/psModules/src/objects/pmSourceOutputs.c	(revision 33596)
+++ branches/meh_branches/ppstack_test/psModules/src/objects/pmSourceOutputs.c	(revision 34041)
@@ -40,5 +40,5 @@
 #include "pmSourceOutputs.h"
 
-bool pmSourceOutputsCommonValues (short *nImageOverlap, float *magOffset, float *zeroptErr, float *fwhmMajor, float *fwhmMinor, pmReadout *readout, psMetadata *header) {
+bool pmSourceOutputsCommonValues (float *magOffset, float *zeroptErr, float *fwhmMajor, float *fwhmMinor, pmReadout *readout, psMetadata *header) {
 
     pmFPA  *fpa  = readout->parent->parent->parent;
@@ -72,5 +72,4 @@
     }
 
-    *nImageOverlap = psMetadataLookupS32 (&status2, header, "NINPUTS");
     return true;
 
Index: branches/meh_branches/ppstack_test/psModules/src/objects/pmSourceOutputs.h
===================================================================
--- branches/meh_branches/ppstack_test/psModules/src/objects/pmSourceOutputs.h	(revision 33596)
+++ branches/meh_branches/ppstack_test/psModules/src/objects/pmSourceOutputs.h	(revision 34041)
@@ -56,5 +56,5 @@
 } pmSourceOutputsMoments;
 
-bool pmSourceOutputsCommonValues (short *nImageOverlap, float *magOffset, float *zeroptErr, float *fwhmMajor, float *fwhmMinor, pmReadout *readout, psMetadata *header);
+bool pmSourceOutputsCommonValues (float *magOffset, float *zeroptErr, float *fwhmMajor, float *fwhmMinor, pmReadout *readout, psMetadata *header);
 
 bool pmSourceOutputsSetValues (pmSourceOutputs *outputs, pmSource *source, pmChip *chip, float fwhmMajor, float fwhmMinor, float magOffset);
Index: branches/meh_branches/ppstack_test/psModules/src/objects/pmSourcePhotometry.c
===================================================================
--- branches/meh_branches/ppstack_test/psModules/src/objects/pmSourcePhotometry.c	(revision 33596)
+++ branches/meh_branches/ppstack_test/psModules/src/objects/pmSourcePhotometry.c	(revision 34041)
@@ -899,15 +899,42 @@
 }
 
+# if (HAVE_MODEL_VAR) 
+double pmSourceModelWeight(const pmSource *Mi, int term, const pmSourceFitVarMode fitVarMode, const float covarFactor, psImageMaskType maskVal)
+# else
 double pmSourceModelWeight(const pmSource *Mi, int term, const bool unweighted_sum, const float covarFactor, psImageMaskType maskVal)
+# endif
 {
     PS_ASSERT_PTR_NON_NULL(Mi, NAN);
+# if (HAVE_MODEL_VAR) 
+    double flux = 0, wt = 1.0, factor = 0;
+# else
     double flux = 0, wt = 0, factor = 0;
+# endif
 
     const psImage *Pi = Mi->modelFlux;
     assert (Pi != NULL);
+
+# if (HAVE_MODEL_VAR)
+    const psImage *Wi = NULL;
+    switch (fitVarMode) {
+      case PM_SOURCE_PHOTFIT_CONST:
+	break;
+      case PM_SOURCE_PHOTFIT_IMAGE_VAR:
+	Wi = Mi->variance;
+        psAssert (Wi, "programming error");
+	break;
+      case PM_SOURCE_PHOTFIT_MODEL_VAR:
+	Wi = Mi->modelVar;
+        psAssert (Wi, "programming error");
+	break;
+      case PM_SOURCE_PHOTFIT_NONE:
+	psAbort("programming error");
+    }	
+# else
     const psImage *Wi = Mi->variance;
     if (!unweighted_sum) {
         assert (Wi != NULL);
     }
+# endif
     const psImage *Ti = Mi->maskObj;
     assert (Ti != NULL);
@@ -917,4 +944,10 @@
             if (Ti->data.PS_TYPE_IMAGE_MASK_DATA[yi][xi] & maskVal)
                 continue;
+# if (HAVE_MODEL_VAR)
+            if (fitVarMode != PM_SOURCE_PHOTFIT_CONST) {
+                wt = covarFactor * Wi->data.F32[yi][xi];
+                if (wt == 0) continue;
+            }
+# else
             if (!unweighted_sum) {
                 wt = covarFactor * Wi->data.F32[yi][xi];
@@ -922,4 +955,5 @@
                     continue;
             }
+# endif
 
             switch (term) {
@@ -937,4 +971,8 @@
             }
 
+# if (HAVE_MODEL_VAR)
+	    // wt is 1.0 for CONST
+	    flux += (factor * Pi->data.F32[yi][xi]) / wt;
+# else
             if (unweighted_sum) {
                 flux += (factor * Pi->data.F32[yi][xi]);
@@ -942,4 +980,5 @@
                 flux += (factor * Pi->data.F32[yi][xi]) / wt;
             }
+# endif
         }
     }
@@ -947,5 +986,9 @@
 }
 
+# if (HAVE_MODEL_VAR)
+double pmSourceModelDotModel (const pmSource *Mi, const pmSource *Mj, const pmSourceFitVarMode fitVarMode, const float covarFactor, psImageMaskType maskVal)
+# else
 double pmSourceModelDotModel (const pmSource *Mi, const pmSource *Mj, const bool unweighted_sum, const float covarFactor, psImageMaskType maskVal)
+# endif
 {
     PS_ASSERT_PTR_NON_NULL(Mi, NAN);
@@ -955,5 +998,10 @@
     int xIs, xJs, yIs, yJs;
     int xIe, yIe;
+# if (HAVE_MODEL_VAR)
+    double flux;
+    double wt = 1.0;
+# else
     double flux, wt;
+# endif
 
     const psImage *Pi = Mi->modelFlux;
@@ -962,8 +1010,26 @@
     assert (Pj != NULL);
 
+# if (HAVE_MODEL_VAR)
+    const psImage *Wi = NULL;
+    switch (fitVarMode) {
+      case PM_SOURCE_PHOTFIT_CONST:
+	break;
+      case PM_SOURCE_PHOTFIT_IMAGE_VAR:
+	Wi = Mi->variance;
+        psAssert (Wi, "programming error");
+	break;
+      case PM_SOURCE_PHOTFIT_MODEL_VAR:
+	Wi = Mi->modelVar;
+        psAssert (Wi, "programming error");
+	break;
+      case PM_SOURCE_PHOTFIT_NONE:
+	psAbort("programming error");
+    }	
+# else
     const psImage *Wi = Mi->variance;
     if (!unweighted_sum) {
         assert (Wi != NULL);
     }
+# endif
 
     const psImage *Ti = Mi->maskObj;
@@ -995,4 +1061,22 @@
                 continue;
 
+# if (HAVE_MODEL_VAR)
+	    float value = (Pi->data.F32[yi][xi] * Pj->data.F32[yj][xj]);
+	    switch (fitVarMode) {
+	      case PM_SOURCE_PHOTFIT_CONST:
+		wt = 1.0;
+		break;
+	      case PM_SOURCE_PHOTFIT_IMAGE_VAR:
+	      case PM_SOURCE_PHOTFIT_MODEL_VAR:
+                wt = covarFactor * Wi->data.F32[yi][xi];
+		break;
+	      case PM_SOURCE_PHOTFIT_NONE:
+		psAbort("programming error");
+	    }
+            // skip pixels with nonsense weight values
+	    if (wt <= 0) continue;
+
+	    flux += value / wt;
+# else
             // XXX skip the nonsense weight pixels?
             if (unweighted_sum) {
@@ -1004,4 +1088,5 @@
                 }
             }
+# endif
         }
     }
@@ -1009,5 +1094,9 @@
 }
 
+# if (HAVE_MODEL_VAR)
+double pmSourceDataDotModel (const pmSource *Mi, const pmSource *Mj, const pmSourceFitVarMode fitVarMode, const float covarFactor, psImageMaskType maskVal)
+# else
 double pmSourceDataDotModel (const pmSource *Mi, const pmSource *Mj, const bool unweighted_sum, const float covarFactor, psImageMaskType maskVal)
+# endif
 {
     PS_ASSERT_PTR_NON_NULL(Mi, NAN);
@@ -1017,5 +1106,10 @@
     int xIs, xJs, yIs, yJs;
     int xIe, yIe;
+# if (HAVE_MODEL_VAR)
+    double flux;
+    double wt = 1.0;
+# else
     double flux, wt;
+# endif
 
     const psImage *Pi = Mi->pixels;
@@ -1024,8 +1118,26 @@
     assert (Pj != NULL);
 
+# if (HAVE_MODEL_VAR)
+    const psImage *Wi = NULL;
+    switch (fitVarMode) {
+      case PM_SOURCE_PHOTFIT_CONST:
+	break;
+      case PM_SOURCE_PHOTFIT_IMAGE_VAR:
+	Wi = Mi->variance;
+        psAssert (Wi, "programming error");
+	break;
+      case PM_SOURCE_PHOTFIT_MODEL_VAR:
+	Wi = Mi->modelVar;
+        psAssert (Wi, "programming error");
+	break;
+      case PM_SOURCE_PHOTFIT_NONE:
+	psAbort("programming error");
+    }	
+# else
     const psImage *Wi = Mi->variance;
     if (!unweighted_sum) {
         assert (Wi != NULL);
     }
+# endif
 
     const psImage *Ti = Mi->maskObj;
@@ -1057,4 +1169,23 @@
                 continue;
 
+# if (HAVE_MODEL_VAR)
+	    float value = (Pi->data.F32[yi][xi] * Pj->data.F32[yj][xj]);
+	    switch (fitVarMode) {
+	      case PM_SOURCE_PHOTFIT_CONST:
+		wt = 1.0;
+		break;
+	      case PM_SOURCE_PHOTFIT_IMAGE_VAR:
+	      case PM_SOURCE_PHOTFIT_MODEL_VAR:
+                wt = covarFactor * Wi->data.F32[yi][xi];
+		break;
+	      case PM_SOURCE_PHOTFIT_NONE:
+		psAbort("programming error");
+	    }
+            // skip pixels with nonsense weight values
+	    if (wt <= 0) continue;
+
+	    flux += value / wt;
+
+# else
             // XXX skip the nonsense weight pixels?
             if (unweighted_sum) {
@@ -1066,4 +1197,5 @@
                 }
             }
+# endif
         }
     }
Index: branches/meh_branches/ppstack_test/psModules/src/objects/pmSourcePhotometry.h
===================================================================
--- branches/meh_branches/ppstack_test/psModules/src/objects/pmSourcePhotometry.h	(revision 33596)
+++ branches/meh_branches/ppstack_test/psModules/src/objects/pmSourcePhotometry.h	(revision 34041)
@@ -38,4 +38,13 @@
 } pmSourcePhotometryMode;
 
+# if (HAVE_MODEL_VAR)
+typedef enum {
+    PM_SOURCE_PHOTFIT_NONE       = 0,
+    PM_SOURCE_PHOTFIT_CONST      = 1,
+    PM_SOURCE_PHOTFIT_IMAGE_VAR  = 2,
+    PM_SOURCE_PHOTFIT_MODEL_VAR  = 3,
+} pmSourceFitVarMode;
+# endif
+
 bool pmSourcePhotometryModel(
     float *fitMag,                      ///< integrated fit magnitude
@@ -75,7 +84,13 @@
 bool pmSourceMeasureDiffStats (pmSource *source, psImageMaskType maskVal, psImageMaskType markVal);
 
+# if (HAVE_MODEL_VAR)
+double pmSourceDataDotModel (const pmSource *Mi, const pmSource *Mj, const pmSourceFitVarMode fitVarMode, const float covarFactor, psImageMaskType maskVal);
+double pmSourceModelDotModel (const pmSource *Mi, const pmSource *Mj, const pmSourceFitVarMode fitVarMode, const float covarFactor, psImageMaskType maskVal);
+double pmSourceModelWeight(const pmSource *Mi, int term, const pmSourceFitVarMode fitVarMode, const float covarFactor, psImageMaskType maskVal);
+# else
 double pmSourceDataDotModel (const pmSource *Mi, const pmSource *Mj, const bool unweighted_sum, const float covarFactor, psImageMaskType maskVal);
 double pmSourceModelDotModel (const pmSource *Mi, const pmSource *Mj, const bool unweighted_sum, const float covarFactor, psImageMaskType maskVal);
 double pmSourceModelWeight(const pmSource *Mi, int term, const bool unweighted_sum, const float covarFactor, psImageMaskType maskVal);
+# endif
 
 bool pmSourceNeighborFlags (pmSource *source);