Index: branches/ipp-magic-v0/ppSub/src/ppSubDefineOutput.c
===================================================================
--- branches/ipp-magic-v0/ppSub/src/ppSubDefineOutput.c	(revision 21361)
+++ branches/ipp-magic-v0/ppSub/src/ppSubDefineOutput.c	(revision 21361)
@@ -0,0 +1,70 @@
+/** @file ppSubDefineOutput.c
+ *
+ *  @brief
+ *
+ *  @ingroup ppSub
+ *
+ *  @author IfA
+ *  @version $Revision: 1.2 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2009-02-06 01:37:17 $
+ *  Copyright 2009 Institute for Astronomy, University of Hawaii
+ */
+
+#include "ppSub.h"
+
+bool ppSubDefineOutput (pmConfig *config, const pmFPAview *view) {
+
+    bool mdok;
+
+    // generate an output readout 
+    pmCell *outCell = pmFPAfileThisCell(config->files, view, "PPSUB.OUTPUT"); // Output cell
+    pmReadout *outRO = pmReadoutAlloc(outCell); // Output readout: subtraction
+    pmFPA *outFPA = outCell->parent->parent; // Output FPA
+    pmHDU *outHDU = outFPA->hdu; // Output HDU
+    if (!outHDU->header) {
+        outHDU->header = psMetadataAlloc();
+    }
+
+    // convolved input images
+    pmReadout *inConv = pmFPAfileThisReadout(config->files, view, "PPSUB.INPUT.CONV"); // Input readout
+    pmReadout *refConv = pmFPAfileThisReadout(config->files, view, "PPSUB.REF.CONV"); // Reference readout
+
+    // Add kernel descrption to header.  We don't know which readout gets the kernels because
+    // it depends on which is larger (the choice is set in ppSubMatchPSFs)
+    bool inputHasKernel = true;
+    pmSubtractionKernels *kernels = psMetadataLookupPtr(&mdok, inConv->analysis, PM_SUBTRACTION_ANALYSIS_KERNEL); // Subtraction kernel
+    if (!kernels) {
+        kernels = psMetadataLookupPtr(&mdok, refConv->analysis, PM_SUBTRACTION_ANALYSIS_KERNEL);
+	inputHasKernel = false;
+    }
+    if (!kernels) {
+        psError(PS_ERR_UNEXPECTED_NULL, true, "Unable to find SUBTRACTION.KERNEL");
+        psFree(inConv);
+        psFree(refConv);
+        psFree(outRO);
+        return false;
+    }
+    psMetadataAddStr(outHDU->header, PS_LIST_TAIL, "PPSUB.KERNEL", 0, "Subtraction kernel", kernels->description);
+
+    if (inputHasKernel) {
+	outRO->analysis = psMetadataCopy(outRO->analysis, inConv->analysis);
+    } else {
+	outRO->analysis = psMetadataCopy(outRO->analysis, refConv->analysis);
+    } 
+
+#ifdef TESTING
+    psImage *kernelImage = NULL;
+
+    if (inputHasKernel) {
+	kernelImage = psMetadataLookupPtr(&mdok, inConv->analysis, "SUBTRACTION.KERNEL.IMAGE"); // Image of the kernels
+    } else {
+        kernelImage = psMetadataLookupPtr(&mdok, refConv->analysis, "SUBTRACTION.KERNEL.IMAGE");
+    }
+
+    psFits *fits = psFitsOpen("kernel.fits", "w");
+    psFitsWriteImage(fits, NULL, kernelImage, 0, NULL);
+    psFitsClose(fits);
+#endif
+
+    return true;
+}
Index: branches/ipp-magic-v0/ppSub/src/ppSubExtras.c
===================================================================
--- branches/ipp-magic-v0/ppSub/src/ppSubExtras.c	(revision 21361)
+++ branches/ipp-magic-v0/ppSub/src/ppSubExtras.c	(revision 21361)
@@ -0,0 +1,38 @@
+/** @file ppSubExtras.c
+ *
+ *  @brief
+ *
+ *  @ingroup ppSub
+ *
+ *  @author IfA
+ *  @version $Revision: 1.2 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2009-02-06 01:37:17 $
+ *  Copyright 2009 Institute for Astronomy, University of Hawaii
+ */
+
+#include "ppSub.h"
+
+// This file contains minor functions used in ppSub.  some of these should be pushed into
+// psModules or psLib.
+
+/**
+ * Copy every instance of a single keyword from one metadata to another
+ */
+bool psMetadataCopySingle(psMetadata *target, psMetadata *source, const char *name)
+{
+    PS_ASSERT_METADATA_NON_NULL(target, false);
+    PS_ASSERT_METADATA_NON_NULL(source, false);
+    PS_ASSERT_STRING_NON_EMPTY(name, false);
+
+    psString regex = NULL;      // Regular expression
+    psStringAppend(&regex, "^%s$", name);
+    psMetadataIterator *iter = psMetadataIteratorAlloc(source, PS_LIST_HEAD, regex); // Iterator
+    psFree(regex);
+    psMetadataItem *item;       // Item from iteration
+    while ((item = psMetadataGetAndIncrement(iter))) {
+        psMetadataAddItem(target, item, PS_LIST_TAIL, PS_META_DUPLICATE_OK);
+    }
+    psFree(iter);
+
+    return true;
+}
Index: branches/ipp-magic-v0/ppSub/src/ppSubMakePSF.c
===================================================================
--- branches/ipp-magic-v0/ppSub/src/ppSubMakePSF.c	(revision 21361)
+++ branches/ipp-magic-v0/ppSub/src/ppSubMakePSF.c	(revision 21361)
@@ -0,0 +1,93 @@
+/** @file ppSubMakePSF.c
+ *
+ *  @brief
+ *
+ *  @ingroup ppSub
+ *
+ *  @author IfA
+ *  @version $Revision: 1.2 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2009-02-06 01:37:17 $
+ *  Copyright 2009 Institute for Astronomy, University of Hawaii
+ */
+
+#include "ppSub.h"
+
+/**
+ * Photometry stage 1: measure the PSF from the minuend image 
+ */
+bool ppSubMakePSF (pmConfig *config, const pmFPAview *view) {
+
+    bool mdok = false;
+
+    // Photometry is to be performed in two stages:
+    // 1. Measure the PSF using the PSF-matched images
+    // 2. Find and measure sources on the subtracted image
+    psTimerStart("PPSUB_PHOT");
+
+    // Look up recipe values
+    psMetadata *recipe = psMetadataLookupMetadata(NULL, config->recipes, PPSUB_RECIPE); // Recipe for ppSub
+    psAssert(recipe, "We checked this earlier, so it should be here.");
+
+    psMetadata *psphotRecipe = psMetadataLookupMetadata(NULL, config->recipes, PSPHOT_RECIPE); // Recipe for psphot
+    psAssert(recipe, "We checked this earlier, so it should be here.");
+
+    if (!psMetadataLookupBool(NULL, recipe, "PHOTOMETRY")) return true;
+
+    bool reverse = psMetadataLookupBool(NULL, config->arguments, "REVERSE"); // Reverse sense of subtraction?
+
+    pmReadout *minuend = NULL;
+    pmFPAfile *minuendFile = NULL;
+    if (reverse) {
+	minuend = pmFPAfileThisReadout(config->files, view, "PPSUB.REF.CONV");
+	minuendFile = psMetadataLookupPtr(&mdok, config->files, "PPSUB.REF.CONV");
+    } else {
+	minuend = pmFPAfileThisReadout(config->files, view, "PPSUB.INPUT.CONV");
+	minuendFile = psMetadataLookupPtr(&mdok, config->files, "PPSUB.INPUT.CONV");
+    }
+
+    // supply the minuend pmFPAfile to psphot as PSPHOT.INPUT:
+    psMetadataAddPtr (config->files, PS_LIST_TAIL, "PSPHOT.INPUT", PS_DATA_UNKNOWN | PS_META_REPLACE, "psphot input : view on another pmFPAfile", minuendFile);
+
+    // old-style variance renormalization
+    if (!ppSubReadoutRenormPhot (config, recipe, minuend)) {
+	psError(PS_ERR_UNKNOWN, false, "failure in renormalization");
+	return false;
+    }
+
+    // extract the loaded sources from the associated readout
+    pmReadout *sourcesRO = pmFPAfileThisReadout(config->files, view, "PPSUB.SOURCES");
+    psArray *sources = psMetadataLookupPtr(&mdok, sourcesRO->analysis, "PSPHOT.SOURCES");
+
+    // generate PSF from the supplied sources (assumes image is background-subtracted)
+    if (!psphotReadoutFindPSF(config, view, sources)) {
+	psError(PS_ERR_UNKNOWN, false, "Unable to perform photometry on subtracted image.");
+	return false;
+    }
+
+    // Record the FWHM in the output header
+    pmReadout *outRO = pmFPAfileThisReadout(config->files, view, "PPSUB.OUTPUT"); // Output readout
+    pmHDU *hdu = pmHDUFromCell(outRO->parent); // HDU with header
+    psMetadataItemSupplement(hdu->header, psphotRecipe, "FWHM_MAJ");
+    psMetadataItemSupplement(hdu->header, psphotRecipe, "FWHM_MIN");
+
+    return true;
+}
+
+// XXX we used to need this, is it still needed?
+
+// pmCell *photCell = pmFPAfileThisCell(config->files, view, "PSPHOT.INPUT");
+// pmCellFreeReadouts(photCell);
+
+// if (!pmFPAfileDropInternal(config->files, "PSPHOT.BACKMDL") ||
+//     !pmFPAfileDropInternal (config->files, "PSPHOT.BACKMDL.STDEV") ||
+//     !pmFPAfileDropInternal (config->files, "PSPHOT.BACKGND")) {
+//     psError(PS_ERR_UNKNOWN, false, "Unable to drop PSPHOT internal files.");
+//     return false;
+// }
+
+// Blow away the sources psphot found --- they're irrelevant for the subtraction
+// XXX is this still needed?  These are now probably not being set 
+// pmReadout *photRO = pmFPAviewThisReadout(view, photFile->fpa); // Readout with sources
+// psMetadataRemoveKey(photRO->analysis, "PSPHOT.SOURCES");
+// psMetadataRemoveKey(photRO->analysis, "PSPHOT.HEADER");
+
Index: branches/ipp-magic-v0/ppSub/src/ppSubMatchPSFs.c
===================================================================
--- branches/ipp-magic-v0/ppSub/src/ppSubMatchPSFs.c	(revision 21361)
+++ branches/ipp-magic-v0/ppSub/src/ppSubMatchPSFs.c	(revision 21361)
@@ -0,0 +1,134 @@
+/** @file ppSubMatchPSFs.c
+ *
+ *  @brief
+ *
+ *  @ingroup ppSub
+ *
+ *  @author IfA
+ *  @version $Revision: 1.2 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2009-02-06 01:37:17 $
+ *  Copyright 2009 Institute for Astronomy, University of Hawaii
+ */
+
+#include "ppSub.h"
+
+/** Generate the PSF-matching kernel and convolve the images as needed.  Most of this function
+ * involves looking up the parameters in the recipe and supplying them to the function
+ * pmSubtractionMatch();
+ */
+
+bool ppSubMatchPSFs (pmConfig *config, const pmFPAview *view) {
+
+    bool mdok;                          // Status of MD lookup
+
+    // Look up recipe values
+    psMetadata *recipe = psMetadataLookupMetadata(NULL, config->recipes, PPSUB_RECIPE); // Recipe for ppSim
+    psAssert(recipe, "We checked this earlier, so it should be here.");
+
+    // input images
+    pmReadout *inRO = pmFPAfileThisReadout(config->files, view, "PPSUB.INPUT"); // Input readout
+    pmReadout *refRO = pmFPAfileThisReadout(config->files, view, "PPSUB.REF"); // Reference readout
+
+    // get or generate output image holders
+    pmReadout *inConv = pmFPAfileThisReadout(config->files, view, "PPSUB.INPUT.CONV"); // Input convolved readout
+    if (!inConv) {
+	pmCell *cell = pmFPAfileThisCell(config->files, view, "PPSUB.INPUT.CONV"); // Input convolved readout
+	inConv = pmReadoutAlloc(cell); // Convolved version of input
+    }
+    pmReadout *refConv = pmFPAfileThisReadout(config->files, view, "PPSUB.REF.CONV"); // Reference convolved readout
+    if (!refConv) {
+	pmCell *cell = pmFPAfileThisCell(config->files, view, "PPSUB.REF.CONV"); // Input convolved readout
+	refConv = pmReadoutAlloc(cell); // Convolved version of input
+    }
+
+    // options which control the generation of optimal parameters
+    bool optimum = psMetadataLookupBool(&mdok, recipe, "OPTIMUM"); // Derive optimum parameters?
+    float optMin = psMetadataLookupF32(&mdok, recipe, "OPTIMUM.MIN"); // Minimum width for search
+    float optMax = psMetadataLookupF32(&mdok, recipe, "OPTIMUM.MAX"); // Maximum width for search
+    float optStep = psMetadataLookupF32(&mdok, recipe, "OPTIMUM.STEP"); // Step for search
+
+    // Vector with FWHMs for optimum search
+    psVector *optWidths = NULL;
+    if (optimum) {
+        optWidths = psVectorCreate(optWidths, optMin, optMax, optStep, PS_TYPE_F32);
+    }
+
+    // Sources in image : these must be loaded from previous analysis stages
+    pmReadout *sourcesRO = pmFPAfileThisReadout(config->files, view, "PPSUB.SOURCES"); // Readout with sources 
+    psArray *sources = NULL;            
+    if (sourcesRO) {
+        sources = psMetadataLookupPtr(&mdok, sourcesRO->analysis, "PSPHOT.SOURCES");
+    }
+
+    int footprint = psMetadataLookupS32(NULL, recipe, "STAMP.FOOTPRINT"); // Stamp half-size
+    int stride = psMetadataLookupS32(NULL, recipe, "STRIDE"); // Size of convolution patches
+    float regionSize = psMetadataLookupF32(NULL, recipe, "REGION.SIZE"); // Size of iso-kernel regs
+    float spacing = psMetadataLookupF32(NULL, recipe, "STAMP.SPACING"); // Typical stamp spacing
+    float threshold = psMetadataLookupF32(NULL, recipe, "STAMP.THRESHOLD"); // Threshold for stmps
+
+    const char *stampsName = psMetadataLookupStr(&mdok, config->arguments, "STAMPS"); // Filename for stamps
+
+    const char *typeStr = psMetadataLookupStr(NULL, recipe, "KERNEL.TYPE"); // Kernel type
+    psAssert(typeStr, "We put it here in ppSubArguments.c");
+
+    pmSubtractionKernelsType type = pmSubtractionKernelsTypeFromString(typeStr); // Type of kernel
+    if (type == PM_SUBTRACTION_KERNEL_NONE) {
+        psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Unrecognised kernel type: %s", typeStr);
+        return false;
+    }
+
+    int size = psMetadataLookupS32(NULL, recipe, "KERNEL.SIZE"); // Kernel half-size
+    int order = psMetadataLookupS32(NULL, recipe, "SPATIAL.ORDER"); // Spatial polynomial order
+    psVector *widths = psMetadataLookupPtr(NULL, recipe, "ISIS.WIDTHS"); // ISIS Gaussian widths
+    psVector *orders = psMetadataLookupPtr(NULL, recipe, "ISIS.ORDERS"); // ISIS Polynomial orders
+
+    int inner = psMetadataLookupS32(NULL, recipe, "INNER"); // Inner radius
+    int ringsOrder = psMetadataLookupS32(NULL, recipe, "RINGS.ORDER"); // RINGS polynomial order
+    int binning = psMetadataLookupS32(NULL, recipe, "SPAM.BINNING"); // Binning for SPAM kernel
+    float penalty = psMetadataLookupF32(NULL, recipe, "PENALTY"); // Penalty for wideness
+
+    int optOrder = psMetadataLookupS32(&mdok, recipe, "OPTIMUM.ORDER"); // Order for search
+    float optThresh = psMetadataLookupF32(&mdok, recipe, "OPTIMUM.TOL"); // Tolerance for search
+
+    int iter = psMetadataLookupS32(NULL, recipe, "ITER"); // Rejection iterations
+    float rej = psMetadataLookupF32(NULL, recipe, "REJ"); // Rejection threshold
+    float sys = psMetadataLookupF32(NULL, recipe, "SYS"); // Relative systematic error in kernel
+
+    float badFrac = psMetadataLookupF32(NULL, recipe, "BADFRAC"); // Maximum bad fraction
+    float poorFrac = psMetadataLookupF32(&mdok, recipe, "POOR.FRACTION"); // Fraction for "poor"
+
+    // XXX EAM : do we need to / want to define different values for BAD and POOR subtraction vs BAD and POOR warp?
+    psImageMaskType maskVal = pmConfigMaskGet("MASK.VALUE", config); // Bits to mask going in to pmSubtractionMatch
+    psImageMaskType maskPoor = pmConfigMaskGet("POOR.WARP", config); // Bits to mask for poor pixels
+    psImageMaskType maskBad = pmConfigMaskGet("BLANK", config); // Bits to mask for bad pixels
+
+    bool dual = psMetadataLookupBool(&mdok, recipe, "DUAL"); // Dual convolution?
+    pmSubtractionMode subMode = dual ? PM_SUBTRACTION_MODE_DUAL : PM_SUBTRACTION_MODE_UNSURE; // Subtracn mode
+
+    int threads = psMetadataLookupS32(NULL, config->arguments, "-threads"); // Number of threads
+    if (threads > 0) {
+        pmSubtractionThreadsInit(inRO, refRO);
+    }
+
+    // Match the PSFs
+    if (!pmSubtractionMatch(inConv, refConv, inRO, refRO, footprint, stride, regionSize, spacing, threshold,
+                            sources, stampsName, type, size, order, widths, orders, inner, ringsOrder,
+                            binning, penalty, optimum, optWidths, optOrder, optThresh, iter, rej, sys,
+                            maskVal, maskBad, maskPoor, poorFrac, badFrac, subMode)) {
+        psError(PS_ERR_UNKNOWN, false, "Unable to match images.");
+        return false;
+    }
+
+    psFree(optWidths);
+
+    pmSubtractionThreadsFinalize(inRO, refRO);
+
+    // XXX drop the pixels associated with inRO and refRO (now that we have inConv and refConf)
+
+    psphotSaveImage (NULL, inRO->image, "inRO.fits");
+    psphotSaveImage (NULL, refRO->image, "refRO.fits");
+    psphotSaveImage (NULL, inConv->image, "inConv.fits");
+    psphotSaveImage (NULL, refConv->image, "refConv.fits");
+
+    return true;
+}
Index: branches/ipp-magic-v0/ppSub/src/ppSubSetMasks.c
===================================================================
--- branches/ipp-magic-v0/ppSub/src/ppSubSetMasks.c	(revision 21361)
+++ branches/ipp-magic-v0/ppSub/src/ppSubSetMasks.c	(revision 21361)
@@ -0,0 +1,115 @@
+/** @file ppSubSetMasks.c
+ *
+ *  @brief
+ *
+ *  @ingroup ppSub
+ *
+ *  @author IfA
+ *  @version $Revision: 1.2 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2009-02-06 01:37:17 $
+ *  Copyright 2009 Institute for Astronomy, University of Hawaii
+ */
+
+#include "ppSub.h"
+
+/** this function generates (if needed) and sets or updates the masks for both the input and
+ * reference images.  this function also has the code for interpolation over bad pixels, but it
+ * is currently if-def-ed out
+ */
+
+bool ppSubSetMasks (pmConfig *config, const pmFPAview *view) {
+
+    psImageMaskType maskValue;
+    psImageMaskType markValue;
+    bool mdok = false;
+
+    if (!pmConfigMaskSetBits(&maskValue, &markValue, config)) {
+        psError(PS_ERR_UNKNOWN, false, "Unable to determine mask value.");
+        return false;
+    }
+
+    // set the mask bits needed by psphot (in psphot recipe) 
+    psphotSetMaskRecipe (config, maskValue, markValue);
+
+    // 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.");
+
+    psImageMaskType satValue = pmConfigMaskGet("SAT", config);
+    psAssert (satValue, "SAT must be non-zero");
+
+    psImageMaskType badValue = pmConfigMaskGet("BAD", config);
+    psAssert (badValue, "BAD must be non-zero");
+
+    // input images
+    pmReadout *inRO = pmFPAfileThisReadout(config->files, view, "PPSUB.INPUT"); // Input readout
+    pmReadout *refRO = pmFPAfileThisReadout(config->files, view, "PPSUB.REF"); // Reference readout
+
+    psImage *input = inRO->image;       // Input image
+    psImage *reference = refRO->image;  // Reference image
+    PS_ASSERT_IMAGES_SIZE_EQUAL(input, reference, false);
+
+    // XXX use psImageCopy below to avoid caring about image dimensions
+    int numCols = input->numCols;
+    int numRows = input->numRows;
+
+    // Generate masks if they don't exist
+    if (!inRO->mask) {
+        if (psMetadataLookupBool(NULL, recipe, "MASK.GENERATE")) {
+            pmReadoutSetMask(inRO, satValue, badValue);
+        } else {
+            inRO->mask = psImageAlloc(numCols, numRows, PS_TYPE_IMAGE_MASK);
+            psImageInit(inRO->mask, 0);
+        }
+    }
+    if (!refRO->mask) {
+        if (psMetadataLookupBool(NULL, recipe, "MASK.GENERATE")) {
+            pmReadoutSetMask(refRO, satValue, badValue);
+        } else {
+            refRO->mask = psImageAlloc(numCols, numRows, PS_TYPE_IMAGE_MASK);
+            psImageInit(refRO->mask, 0);
+        }
+    }
+
+    // Mask the NAN values
+    if (!pmReadoutMaskNonfinite(inRO, satValue)) {
+        psError(PS_ERR_UNKNOWN, false, "Unable to mask non-finite pixels in input.");
+        return false;
+    }
+    if (!pmReadoutMaskNonfinite(refRO, satValue)) {
+	psError(PS_ERR_UNKNOWN, false, "Unable to mask non-finite pixels in reference.");
+	return false;
+    }
+
+#if (0)
+    // Look up recipe values
+    psMetadata *recipe = psMetadataLookupMetadata(NULL, config->recipes, PPSUB_RECIPE); // Recipe for ppSim
+    psAssert(recipe, "We checked this earlier, so it should be here.");
+
+    bool mdok = false;
+
+    psString interpModeStr = psMetadataLookupStr(&mdok, recipe, "INTERPOLATION"); // Interpolation mode
+    psImageInterpolateMode interpMode = psImageInterpolateModeFromString(interpModeStr); // Interp
+    if (interpMode == PS_INTERPOLATE_NONE) {
+        psError(PS_ERR_BAD_PARAMETER_VALUE, false, "Unknown interpolation mode: %s", interpModeStr);
+        return false;
+    }
+    float poorFrac = psMetadataLookupF32(&mdok, recipe, "POOR.FRACTION"); // Fraction for "poor"
+
+    psImageMaskType maskPoor = pmConfigMaskGet("POOR.WARP", config); // Bits to mask for poor pixels
+    psImageMaskType maskBad = pmConfigMaskGet("BLANK", config); // Bits to mask for bad pixels
+
+    // Interpolate over bad pixels, so the bad pixels don't explode
+    if (!pmReadoutInterpolateBadPixels(inRO, maskVal, interpMode, poorFrac, maskPoor, maskBad)) {
+	psError(PS_ERR_UNKNOWN, false, "Unable to interpolate bad pixels for input image.");
+	return false;
+    }
+    if (!pmReadoutInterpolateBadPixels(refRO, maskVal, interpMode, poorFrac, maskPoor, maskBad)) {
+	psError(PS_ERR_UNKNOWN, false, "Unable to interpolate bad pixels for reference image.");
+	return false;
+    }
+    maskVal |= maskBad;
+#endif
+
+    return true;
+}