Index: trunk/psModules/src/imcombine/pmSubtraction.c
===================================================================
--- trunk/psModules/src/imcombine/pmSubtraction.c	(revision 29598)
+++ trunk/psModules/src/imcombine/pmSubtraction.c	(revision 29599)
@@ -27,4 +27,5 @@
 #include "pmSubtraction.h"
 
+# define FFT_WINDOW 0
 //#define TESTING
 
@@ -96,5 +97,4 @@
     return;
 }
-
 
 // Generate an image of the solved kernel
@@ -187,9 +187,23 @@
     }
 
+    /* we have three possible implementations for handling flux conservation and photometric scaling.
+       
+       1) the original implementation of Alard-Lupton subtracted the 0th kernel gaussian from
+          the rest of the kernels.  the non-zero integral of the 0th kernel allows for some
+          flux re-scaling if needed.
+
+       2) the original ppSub implementation subtracted a delta function from the kernels.  this
+          seems to make the solution sensitive to noise terms.
+
+       3) accept the measured flux normalization and make all kernels have zero integral, but
+          use a gaussian to zero the flux.  
+
+	  *** the operation below is only needed if we use option (2) ***
+    */
+
     if (normalise) {
         // Put in the normalisation component
         kernel->kernel[0][0] += (wantDual ? 1.0 : p_pmSubtractionSolutionNorm(kernels));
     }
-
     return kernel;
 }
@@ -243,5 +257,9 @@
     psImage *subMask = mask ? psImageSubset(mask, border) : NULL; // Subimage mask
 
+# if (FFT_WINDOW)
+    psImage *convolved = psImageConvolveFFTwithWindow(NULL, subImage, subMask, maskVal, kernel); // Convolution
+# else
     psImage *convolved = psImageConvolveFFT(NULL, subImage, subMask, maskVal, kernel); // Convolution
+# endif
 
     psFree(subImage);
@@ -287,6 +305,11 @@
 
     // XXX Can trim this a little by combining the convolution: only have to take the FFT of the kernel once
+# if (FFT_WINDOW)
+    psImage *convVariance = psImageConvolveFFTwithWindow(NULL, subVariance, subMask, maskVal, kernel); // Convolved variance
+    psImage *convKE = subKE ? psImageConvolveFFTwithWindow(NULL, subKE, subMask, maskVal, kernel) : NULL; // Conv KE
+# else
     psImage *convVariance = psImageConvolveFFT(NULL, subVariance, subMask, maskVal, kernel); // Convolved variance
     psImage *convKE = subKE ? psImageConvolveFFT(NULL, subKE, subMask, maskVal, kernel) : NULL; // Conv KE
+# endif
 
     psFree(subVariance);
@@ -539,5 +562,9 @@
     PS_ASSERT_KERNEL_NON_NULL(kernel, NULL);
 
+# if (FFT_WINDOW)
+    psImage *conv = psImageConvolveFFTwithWindow(NULL, image->image, NULL, 0, kernel); // Convolved image
+# else
     psImage *conv = psImageConvolveFFT(NULL, image->image, NULL, 0, kernel); // Convolved image
+# endif
 
     // note: do not attempt to renormalize kernels here: cannot have different stars with
@@ -978,4 +1005,6 @@
 }
 
+// generate an image of the convolution kernel realized at the given coordinate
+// if 'wantDual' is set, solution2 is supplied
 psImage *pmSubtractionKernelImage(const pmSubtractionKernels *kernels, float x, float y, bool wantDual)
 {
@@ -1285,5 +1314,7 @@
 
             psRegion *subRegion = psRegionAlloc(i, xSubMax, j, ySubMax); // Bounds of subtraction
-            if (threaded) {
+	    // for a TEST, do not run threaded for testing
+            // if (false && threaded) {
+	    if (threaded) {
                 psThreadJob *job = psThreadJobAlloc("PSMODULES_SUBTRACTION_CONVOLVE");
                 psArray *args = job->args;