Index: trunk/ppSub/src/ppSubReadout.c
===================================================================
--- trunk/ppSub/src/ppSubReadout.c	(revision 17835)
+++ trunk/ppSub/src/ppSubReadout.c	(revision 18289)
@@ -12,5 +12,5 @@
 
 #define WCS_TOLERANCE 0.001             // Tolerance for WCS
-//#define TESTING                         // For test output
+#define TESTING                         // For test output
 
 
@@ -73,4 +73,5 @@
     int ringsOrder = psMetadataLookupS32(NULL, config->arguments, "RINGS.ORDER"); // RINGS polynomial order
     int binning = psMetadataLookupS32(NULL, config->arguments, "SPAM.BINNING"); // Binning for SPAM kernel
+    float penalty = psMetadataLookupF32(NULL, config->arguments, "PENALTY"); // Penalty for wideness
     psMaskType maskBad = pmConfigMask(psMetadataLookupStr(NULL, config->arguments, "MASK.BAD"),
                                     config); // Value to mask
@@ -131,6 +132,6 @@
 
     if (!pmSubtractionMatch(inConv, refConv, inRO, refRO, footprint, regionSize, spacing, threshold, sources,
-                            stampsName, type, size, order, widths, orders, inner, ringsOrder,
-                            binning, optimum, optWidths, optOrder, optThresh, iter, rej, maskBad,
+                            stampsName, type, size, order, widths, orders, inner, ringsOrder, binning,
+                            penalty, optimum, optWidths, optOrder, optThresh, iter, rej, maskBad,
                             maskBlank, badFrac, mode)) {
         psError(PS_ERR_UNKNOWN, false, "Unable to match images.");
@@ -184,10 +185,31 @@
 #endif
 
-    // "Subtract" the mask and weight map
+    // Subtraction is: minuend - subtrahend
+    pmReadout *minuend = inConv;
+    pmReadout *subtrahend = refConv;
+
+    if (reverse) {
+        pmReadout *temp = subtrahend;
+        subtrahend = minuend;
+        minuend = temp;
+    }
+
+#ifdef TESTING
+    {
+        pmReadoutMaskApply(minuend, maskBlank);
+        psFits *fits = psFitsOpen("minuend.fits", "w");
+        psFitsWriteImage(fits, NULL, minuend->image, 0, NULL);
+        psFitsClose(fits);
+    }
+    {
+        pmReadoutMaskApply(subtrahend, maskBlank);
+        psFits *fits = psFitsOpen("subtrahend.fits", "w");
+        psFitsWriteImage(fits, NULL, subtrahend->image, 0, NULL);
+            psFitsClose(fits);
+    }
+#endif
+
     outRO->mask = (psImage*)psBinaryOp(outRO->mask, inConv->mask, "|", refConv->mask);
-    if (inConv->weight && refConv->weight) {
-        outRO->weight = (psImage*)psBinaryOp(outRO->weight, inConv->weight, "+", refConv->weight);
-    }
-    outRO->data_exists = outCell->data_exists = outCell->parent->data_exists = true;
+    outRO->data_exists = outCell->data_exists = outCell->parent->data_exists = true; // It'll be there soon
 
     // Photometry is to be performed in two stages:
@@ -198,6 +220,8 @@
     pmPSF *psf = NULL;                  // PSF for photometry
     if (psMetadataLookupBool(NULL, config->arguments, "PHOTOMETRY")) {
-        // We use a summed image as the basis for the PSF: this will have the maximum S/N.
-        outRO->image = (psImage*)psBinaryOp(outRO->image, inConv->image, "+", refConv->image);
+        outRO->image = psImageCopy(outRO->image, minuend->image, PS_TYPE_F32);
+        if (minuend->weight) {
+            outRO->weight = psImageCopy(outRO->weight, minuend->weight, PS_TYPE_F32);
+        }
 
         pmFPAfile *photFile = psMetadataLookupPtr(NULL, config->files, "PSPHOT.INPUT");
@@ -244,33 +268,12 @@
     }
 
-    // Do the subtraction
-    {
-        // Subtraction is: minuend - subtrahend
-        pmReadout *minuend = inConv;
-        pmReadout *subtrahend = refConv;
-
-        if (reverse) {
-            pmReadout *temp = subtrahend;
-            subtrahend = minuend;
-            minuend = temp;
-        }
-
-        outRO->image = (psImage*)psBinaryOp(outRO->image, minuend->image, "-", subtrahend->image);
-
-#ifdef TESTING
-        {
-            psFits *fits = psFitsOpen("minuend.fits", "w");
-            psFitsWriteImage(fits, NULL, minuend->image, 0, NULL);
-            psFitsClose(fits);
-        }
-        {
-            psFits *fits = psFitsOpen("subtrahend.fits", "w");
-            psFitsWriteImage(fits, NULL, subtrahend->image, 0, NULL);
-            psFitsClose(fits);
-        }
-#endif
-
-        pmReadoutMaskApply(outRO, maskBlank);
-    }
+    // Do the actual subtraction
+    outRO->image = (psImage*)psBinaryOp(outRO->image, minuend->image, "-", subtrahend->image);
+    if (inConv->weight && refConv->weight) {
+        outRO->weight = (psImage*)psBinaryOp(outRO->weight, inConv->weight, "+", refConv->weight);
+    }
+    outRO->data_exists = outCell->data_exists = outCell->parent->data_exists = true;
+
+    pmReadoutMaskApply(outRO, maskBlank);
 
     psFree(inConv);
@@ -304,6 +307,4 @@
     }
 
-    psTraceSetLevel("psphot", 6);
-
     // Photometry stage 2: find and measure sources on the subtracted image
     if (psMetadataLookupBool(NULL, config->arguments, "PHOTOMETRY")) {
@@ -314,5 +315,5 @@
         pmReadout *psfRO = pmFPAfileThisReadout(config->files, view, "PSPHOT.PSF.LOAD");
         if (!psfRO) {
-            psError(PS_ERR_UNEXPECTED_NULL, true, "Unable to find file PSPHOT.PSF.LOAD");
+            psError(PS_ERR_UNEXPECTED_NULL, true, "Unable to find readout for file PSPHOT.PSF.LOAD");
             return false;
         }