Index: trunk/psModules/src/imcombine/pmSubtractionKernels.c
===================================================================
--- trunk/psModules/src/imcombine/pmSubtractionKernels.c	(revision 29004)
+++ trunk/psModules/src/imcombine/pmSubtractionKernels.c	(revision 29543)
@@ -183,4 +183,10 @@
 }
 
+# define CENTRAL_DELTA 0
+
+# if (CENTRAL_DELTA)
+
+// XXX *** this code used the central pixel to force zero net flux,
+// Alard actually uses kernel(0) for this purpose (for even order, kernel[i] = kernel'[i] - kernel[0])
 static bool pmSubtractionKernelPreCalcNormalize(pmSubtractionKernels *kernels, pmSubtractionKernelPreCalc *preCalc,
 						int index, int uOrder, int vOrder, float fwhm,
@@ -220,10 +226,10 @@
 	    // Re-normalize 
             // scale2D  = 1.0 / fabs(sum);
-            scale2D  = 1.0 / sqrt(sum2);
+            scale2D  = 1.0 / sqrt(sum2) / PS_SQR(fwhm);
             zeroNull = true;
         } else {
             // Odd functions: choose normalisation so that parameters have about the same strength as for even
             // functions, no subtraction of null pixel because the sum is already (near) zero
-            scale2D = 1.0 / sqrt(sum2);
+            scale2D = 1.0 / sqrt(sum2) / PS_SQR(fwhm);
             zeroNull = false;
         }
@@ -235,10 +241,10 @@
     if (forceZeroNull) {
         // Force rescaling and subtraction of null pixel even though the order doesn't indicate it's even
-        scale2D = 1.0 / fabs(sum);
+        scale2D = 1.0 / fabs(sum) / PS_SQR(fwhm);
         zeroNull = true;
     }
     if (!forceZeroNull && ((uOrder % 2) || (vOrder % 2))) {
         // Odd function
-        scale2D = 1.0 / sqrt(sum2);
+        scale2D = 1.0 / sqrt(sum2) / PS_SQR(fwhm);
     }
 
@@ -255,8 +261,8 @@
     if (zeroNull) {
         // preCalc->kernel->kernel[0][0] -= 1.0;
-        preCalc->kernel->kernel[0][0] -= sum / sqrt (sum2);
-    }
-
-#if 1
+        preCalc->kernel->kernel[0][0] -= sum * scale2D;
+    }
+
+#if 0
     {
         double Sum = 0.0;   // Sum of kernel component
@@ -287,4 +293,102 @@
     return true;
 }
+
+# else /* CENTRAL_DELTA */
+
+static bool zeroIsNormal = false;
+
+// this code uses kernel(0) to force zero flux, and is invalid for other kinds of normalizations
+static bool pmSubtractionKernelPreCalcNormalize(pmSubtractionKernels *kernels, pmSubtractionKernelPreCalc *preCalc,
+						int index, int uOrder, int vOrder, float fwhm,
+						bool AlardLuptonStyle, bool forceZeroNull)
+{
+    // 1) for odd functions: no renormalization 
+    // 2) for even functions, normalize the kernel to unity
+    // 3) for even functions & index > 0, subtract kernel(0)
+
+    // Calculate moments
+    double sum = 0.0, sum2 = 0.0;           // Sum of kernel component
+    float min = INFINITY, max = -INFINITY;  // Minimum and maximum kernel value
+
+    for (int v = preCalc->kernel->yMin; v <= preCalc->kernel->yMax; v++) {
+        for (int u = preCalc->kernel->xMin; u <= preCalc->kernel->xMax; u++) {
+            double value = preCalc->kernel->kernel[v][u];
+            double value2 = PS_SQR(value);
+            sum += value;
+            sum2 += value2;
+            min = PS_MIN(value, min);
+            max = PS_MAX(value, max);
+        }
+    }
+
+#if 0
+    fprintf(stderr, "%d raw: %lf, null: %f, min: %lf, max: %lf\n", index, sum, preCalc->kernel->kernel[0][0], min, max);
+#endif
+
+    float scale2D = 1.0;		// Scaling for 2-D kernels
+    float scale1D = 1.0;		// Scaling for 1-D kernels
+
+    if (uOrder % 2 == 0 && vOrder % 2 == 0) {
+
+	scale2D = 1.0 / sum;		// Scaling for 2-D kernels
+	scale1D = sqrtf(scale2D);		// Scaling for 1-D kernels
+
+	for (int v = preCalc->kernel->yMin; v <= preCalc->kernel->yMax; v++) {
+	    for (int u = preCalc->kernel->xMin; u <= preCalc->kernel->xMax; u++) {
+		preCalc->kernel->kernel[v][u] *= scale2D;
+	    }
+	}
+	if (index == 0) {
+	    zeroIsNormal = true;
+	} else {
+	    psAssert(zeroIsNormal, "failed to normalize zero kernel first");
+	    pmSubtractionKernelPreCalc *zeroKernel = kernels->preCalc->data[0];
+	    psAssert(zeroKernel, "failed to supply zero kernel");
+	    for (int v = preCalc->kernel->yMin; v <= preCalc->kernel->yMax; v++) {
+		for (int u = preCalc->kernel->xMin; u <= preCalc->kernel->xMax; u++) {
+		    preCalc->kernel->kernel[v][u] -= zeroKernel->kernel->kernel[v][u];
+		}
+	    }
+	}
+
+	// XXX why do we bother renormlizing the 1D kernels?  I don't think we use them again...
+	if (preCalc->xKernel) {
+	    psBinaryOp(preCalc->xKernel, preCalc->xKernel, "*", psScalarAlloc(scale1D, PS_TYPE_F32));
+	}
+	if (preCalc->yKernel) {
+	    psBinaryOp(preCalc->yKernel, preCalc->yKernel, "*", psScalarAlloc(scale1D, PS_TYPE_F32));
+	}
+    }
+
+#if 0
+    {
+        double Sum = 0.0;   // Sum of kernel component
+        double Sum2 = 0.0;   // Sum of kernel component
+        float min = INFINITY, max = -INFINITY;  // Minimum and maximum kernel value
+	for (int v = preCalc->kernel->yMin; v <= preCalc->kernel->yMax; v++) {
+	    for (int u = preCalc->kernel->xMin; u <= preCalc->kernel->xMax; u++) {
+		double value = preCalc->kernel->kernel[v][u];
+                Sum += value;
+		Sum2 += PS_SQR(value);
+                min = PS_MIN(preCalc->kernel->kernel[v][u], min);
+                max = PS_MAX(preCalc->kernel->kernel[v][u], max);
+            }
+        }
+        fprintf(stderr, "%d sum: %lf, sum2: %lf, null: %f, min: %lf, max: %lf, scale: %f\n", index, Sum, Sum2, preCalc->kernel->kernel[0][0], min, max, scale2D);
+    }
+#endif
+
+    kernels->widths->data.F32[index] = fwhm;
+    kernels->u->data.S32[index] = uOrder;
+    kernels->v->data.S32[index] = vOrder;
+    if (kernels->preCalc->data[index]) {
+        psFree(kernels->preCalc->data[index]);
+    }
+    kernels->preCalc->data[index] = preCalc;
+    psTrace("psModules.imcombine", 7, "Kernel %d: %f %d %d\n", index, fwhm, uOrder, vOrder);
+
+    return true;
+}
+# endif /* Central Delta */
 
 pmSubtractionKernels *p_pmSubtractionKernelsRawISIS(int size, int spatialOrder,
@@ -603,10 +707,10 @@
     kernels->sampleStamps = NULL;
 
-    kernels->fSigResMean  = NAN;
-    kernels->fSigResStdev = NAN;
-    kernels->fMaxResMean  = NAN;
-    kernels->fMaxResStdev = NAN;
-    kernels->fMinResMean  = NAN;
-    kernels->fMinResStdev = NAN;
+    kernels->fResSigmaMean  = NAN;
+    kernels->fResSigmaStdev = NAN;
+    kernels->fResOuterMean  = NAN;
+    kernels->fResOuterStdev = NAN;
+    kernels->fResTotalMean  = NAN;
+    kernels->fResTotalStdev = NAN;
 
     return kernels;