Index: branches/ipp-magic-v0/ppSub/src/ppSubReadoutSubtract.c
===================================================================
--- branches/ipp-magic-v0/ppSub/src/ppSubReadoutSubtract.c	(revision 21425)
+++ branches/ipp-magic-v0/ppSub/src/ppSubReadoutSubtract.c	(revision 21425)
@@ -0,0 +1,118 @@
+/** @file ppSubReadoutSubtract.c
+ *
+ *  @brief
+ *
+ *  @ingroup ppSub
+ *
+ *  @author IfA
+ *  @version $Revision: 1.4 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2009-02-09 21:26:05 $
+ *  Copyright 2009 Institute for Astronomy, University of Hawaii
+ */
+
+#include "ppSub.h"
+#define WCS_TOLERANCE 0.001             // Tolerance for WCS
+
+bool ppSubReadoutSubtract (pmConfig *config, const pmFPAview *view) {
+
+    bool mdok = false;
+
+    // Look up recipe values
+    psMetadata *recipe = psMetadataLookupMetadata(&mdok, config->recipes, PPSUB_RECIPE); // Recipe for ppSim
+    psAssert(recipe, "We checked this earlier, so it should be here.");
+
+    bool reverse = psMetadataLookupBool(&mdok, config->arguments, "REVERSE"); // Reverse sense of subtraction?
+
+    // Subtraction is: minuend - subtrahend
+    pmReadout *minuend = NULL;
+    pmReadout *subtrahend = NULL;
+    if (reverse) {
+        minuend = pmFPAfileThisReadout(config->files, view, "PPSUB.REF.CONV");
+        subtrahend = pmFPAfileThisReadout(config->files, view, "PPSUB.INPUT.CONV");
+    } else {
+        minuend = pmFPAfileThisReadout(config->files, view, "PPSUB.INPUT.CONV");
+        subtrahend = pmFPAfileThisReadout(config->files, view, "PPSUB.REF.CONV");
+    }
+
+// XXX this should be replaced by SAVE options to the inConv and refConf pmFPAfiles
+#ifdef TESTING
+    {
+        pmReadoutMaskApply(minuend, maskVal);
+        psFits *fits = psFitsOpen("minuend.fits", "w");
+        psFitsWriteImage(fits, NULL, minuend->image, 0, NULL);
+        psFitsClose(fits);
+    }
+    {
+        pmReadoutMaskApply(subtrahend, maskVal);
+        psFits *fits = psFitsOpen("subtrahend.fits", "w");
+        psFitsWriteImage(fits, NULL, subtrahend->image, 0, NULL);
+            psFitsClose(fits);
+    }
+#endif
+
+    // Do the actual subtraction
+    pmReadout *outRO = pmFPAfileThisReadout(config->files, view, "PPSUB.OUTPUT");
+    outRO->image = (psImage*)psBinaryOp(outRO->image, minuend->image, "-", subtrahend->image);
+    if (minuend->variance && subtrahend->variance) {
+        outRO->variance = (psImage*)psBinaryOp(outRO->variance, minuend->variance, "+", subtrahend->variance);
+    }
+    outRO->mask = (psImage*)psBinaryOp(outRO->mask, minuend->mask, "|", subtrahend->mask);
+
+    outRO->data_exists = true;
+    outRO->parent->data_exists = true;
+    outRO->parent->parent->data_exists = true;
+
+    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;
+
+    if (!pmFPACopyConcepts(outFPA, inFPA)) {
+        psError(PS_ERR_UNKNOWN, false, "Unable to copy concepts from input to output.");
+        psFree(outRO);
+        return false;
+    }
+
+    // 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");
+
+    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.");
+        return false;
+    }
+
+    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
+