Index: trunk/ppSub/src/ppSubReadoutSubtract.c
===================================================================
--- trunk/ppSub/src/ppSubReadoutSubtract.c	(revision 21424)
+++ trunk/ppSub/src/ppSubReadoutSubtract.c	(revision 21524)
@@ -6,18 +6,30 @@
  *
  *  @author IfA
- *  @version $Revision: 1.4 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2009-02-09 21:26:05 $
+ *  @version $Revision: 1.5 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2009-02-18 00:31:20 $
  *  Copyright 2009 Institute for Astronomy, University of Hawaii
  */
 
+#ifdef HAVE_CONFIG_H
+#include <config.h>
+#endif
+
+#include <stdio.h>
+#include <pslib.h>
+#include <psmodules.h>
+
 #include "ppSub.h"
+
 #define WCS_TOLERANCE 0.001             // Tolerance for WCS
 
-bool ppSubReadoutSubtract (pmConfig *config, const pmFPAview *view) {
+bool ppSubReadoutSubtract(pmConfig *config, const pmFPAview *view)
+{
+    psAssert(config, "Require configuration");
+    psAssert(view, "Require view");
 
-    bool mdok = false;
 
     // Look up recipe values
-    psMetadata *recipe = psMetadataLookupMetadata(&mdok, config->recipes, PPSUB_RECIPE); // Recipe for ppSim
+    bool mdok = false;                  // Status of MD lookup
+    psMetadata *recipe = psMetadataLookupMetadata(&mdok, config->recipes, PPSUB_RECIPE); // Recipe for ppSub
     psAssert(recipe, "We checked this earlier, so it should be here.");
 
@@ -25,6 +37,6 @@
 
     // Subtraction is: minuend - subtrahend
-    pmReadout *minuend = NULL;
-    pmReadout *subtrahend = NULL;
+    pmReadout *minuend = NULL;          // Positive image
+    pmReadout *subtrahend = NULL;       // Negative image
     if (reverse) {
         minuend = pmFPAfileThisReadout(config->files, view, "PPSUB.REF.CONV");
@@ -59,4 +71,10 @@
     outRO->mask = (psImage*)psBinaryOp(outRO->mask, minuend->mask, "|", subtrahend->mask);
 
+    psArray *covars = psArrayAlloc(2);  // Covariance pseudo-matrices
+    covars->data[0] = psMemIncrRefCounter(minuend->covariance);
+    covars->data[1] = psMemIncrRefCounter(subtrahend->covariance);
+    outRO->covariance = psImageCovarianceSum(covars);
+    psFree(covars);
+
     outRO->data_exists = true;
     outRO->parent->data_exists = true;
@@ -65,11 +83,9 @@
     pmSubtractionVisualShowSubtraction(minuend->image, subtrahend->image, outRO->image);
 
-    // copy concepts from the input to the output (XXX should this always use minuend?)
-    pmFPAfile *inFile = psMetadataLookupPtr (&mdok, config->files, "PPSUB.INPUT");
-    pmFPA *inFPA = inFile->fpa;
-
-    pmFPAfile *outFile = psMetadataLookupPtr (&mdok, config->files, "PPSUB.OUTPUT");
-    pmFPA *outFPA = outFile->fpa;
-
+    // Copy concepts from the input to the output
+    pmFPAfile *inFile = psMetadataLookupPtr(&mdok, config->files, "PPSUB.INPUT"); // Input file
+    pmFPA *inFPA = inFile->fpa;         // Input FPA
+    pmFPAfile *outFile = psMetadataLookupPtr(&mdok, config->files, "PPSUB.OUTPUT"); // Output file
+    pmFPA *outFPA = outFile->fpa;       // Output FPA
     if (!pmFPACopyConcepts(outFPA, inFPA)) {
         psError(PS_ERR_UNKNOWN, false, "Unable to copy concepts from input to output.");
@@ -78,21 +94,17 @@
     }
 
-    // get the HDUs and output Chip to copy the astrometry
-    pmHDU *inHDU = inFPA->hdu;        // input HDU
-    pmHDU *outHDU = outFPA->hdu;
-    pmChip *outChip = pmFPAfileThisChip(config->files, view, "PPSUB.OUTPUT");
-
+    // Copy astrometry over
+    // It should find its way into the output images and photometry
+    pmHDU *inHDU = inFPA->hdu;          // Input HDU
+    pmHDU *outHDU = outFPA->hdu;        // Output HDU
+    pmChip *outChip = pmFPAfileThisChip(config->files, view, "PPSUB.OUTPUT"); // Output chip
     if (!outHDU || !inHDU) {
         psError(PS_ERR_UNKNOWN, false, "Unable to find HDU at FPA level to copy astrometry.");
         return false;
     }
-
-    // Copy astrometry over
-    // It should get into the output images and photometry
     if (!pmAstromReadWCS(outFPA, outChip, inHDU->header, 1.0)) {
         psError(PS_ERR_UNKNOWN, false, "Unable to read WCS astrometry from input FPA.");
         return false;
     }
-
     if (!pmAstromWriteWCS(outHDU->header, outFPA, outChip, WCS_TOLERANCE)) {
         psError(PS_ERR_UNKNOWN, false, "Unable to write WCS astrometry to output FPA.");
@@ -102,17 +114,2 @@
     return true;
 }
-
-// XXX this test code was in place to check for and squash unexpected NANs
-
-#if 0
-    pmReadoutMaskApply(outRO, maskBad);
-
-    for (int y = 0; y < outRO->image->numRows; y++) {
-        for (int x = 0; x < outRO->image->numCols; x++) {
-            if (isnan(outRO->image->data.F32[y][x]) && !(outRO->mask->data.PS_TYPE_IMAGE_MASK_DATA[y][x] & maskVal)) {
-                printf("Unmasked NAN at %d %d --> %d\n", x, y, outRO->mask->data.PS_TYPE_IMAGE_MASK_DATA[y][x]);
-            }
-        }
-    }
-#endif
-