Index: trunk/ppSub/src/ppSubMatchPSFs.c
===================================================================
--- trunk/ppSub/src/ppSubMatchPSFs.c	(revision 21360)
+++ trunk/ppSub/src/ppSubMatchPSFs.c	(revision 21524)
@@ -6,19 +6,24 @@
  *
  *  @author IfA
- *  @version $Revision: 1.2 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2009-02-06 01:37:17 $
+ *  @version $Revision: 1.3 $ $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 <psphot.h>
+
 #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
+bool ppSubMatchPSFs(pmConfig *config, const pmFPAview *view)
+{
+    psAssert(config, "Require configuration");
+    psAssert(view, "Require view");
 
     // Look up recipe values
@@ -26,35 +31,25 @@
     psAssert(recipe, "We checked this earlier, so it should be here.");
 
-    // input images
+    // 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
+    // Output image holders
+    pmReadout *inConv = pmFPAfileThisReadout(config->files, view, "PPSUB.INPUT.CONV"); // Input convolved
     if (!inConv) {
-	pmCell *cell = pmFPAfileThisCell(config->files, view, "PPSUB.INPUT.CONV"); // Input convolved readout
-	inConv = pmReadoutAlloc(cell); // Convolved version of input
+        pmCell *cell = pmFPAfileThisCell(config->files, view, "PPSUB.INPUT.CONV"); // Cell for convolved input
+        inConv = pmReadoutAlloc(cell);
     }
-    pmReadout *refConv = pmFPAfileThisReadout(config->files, view, "PPSUB.REF.CONV"); // Reference convolved readout
+    pmReadout *refConv = pmFPAfileThisReadout(config->files, view, "PPSUB.REF.CONV"); // Reference convolved
     if (!refConv) {
-	pmCell *cell = pmFPAfileThisCell(config->files, view, "PPSUB.REF.CONV"); // Input convolved readout
-	refConv = pmReadoutAlloc(cell); // Convolved version of input
+        pmCell *cell = pmFPAfileThisCell(config->files, view, "PPSUB.REF.CONV"); // Cell for convolved ref.
+        refConv = pmReadoutAlloc(cell);
     }
 
-    // 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
+    bool mdok;                          // Status of MD lookup
 
-    // 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;            
+    // Sources in image, used for stamps: these must be loaded from previous analysis stages
+    pmReadout *sourcesRO = pmFPAfileThisReadout(config->files, view, "PPSUB.SOURCES"); // Readout with sources
+    psArray *sources = NULL;            // Sources in image; used for stamps
     if (sourcesRO) {
         sources = psMetadataLookupPtr(&mdok, sourcesRO->analysis, "PSPHOT.SOURCES");
@@ -66,10 +61,8 @@
     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) {
@@ -87,8 +80,4 @@
     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
@@ -98,6 +87,18 @@
     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
+    // 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
+    psVector *optWidths = NULL;         // Vector with FWHMs for optimum search
+    if (optimum) {
+        optWidths = psVectorCreate(optWidths, optMin, optMax, optStep, PS_TYPE_F32);
+    }
+    int optOrder = psMetadataLookupS32(&mdok, recipe, "OPTIMUM.ORDER"); // Order for search
+    float optThresh = psMetadataLookupF32(&mdok, recipe, "OPTIMUM.TOL"); // Tolerance for search
+
+    // XXX Do we need/want to define different values for BAD and POOR subtraction vs BAD and POOR warp?
+    psImageMaskType maskVal = pmConfigMaskGet("MASK.VALUE", config); // Bits to mask in inputs
     psImageMaskType maskPoor = pmConfigMaskGet("POOR.WARP", config); // Bits to mask for poor pixels
     psImageMaskType maskBad = pmConfigMaskGet("BLANK", config); // Bits to mask for bad pixels
@@ -124,4 +125,7 @@
     pmSubtractionThreadsFinalize(inRO, refRO);
 
+    psImageCovarianceTransfer(inConv->variance, inConv->covariance);
+    psImageCovarianceTransfer(refConv->variance, refConv->covariance);
+
     // XXX drop the pixels associated with inRO and refRO (now that we have inConv and refConf)