Index: /branches/eam_branches/ipp-20100823/psModules/src/imcombine/pmSubtractionEquation.c
===================================================================
--- /branches/eam_branches/ipp-20100823/psModules/src/imcombine/pmSubtractionEquation.c	(revision 29505)
+++ /branches/eam_branches/ipp-20100823/psModules/src/imcombine/pmSubtractionEquation.c	(revision 29506)
@@ -21,5 +21,5 @@
 
 // XXX TEST:
-# define USE_WINDOW                      // window to avoid neighbor contamination
+//# define USE_WINDOW                      // window to avoid neighbor contamination
 
 # define PENALTY false
@@ -1091,5 +1091,5 @@
         }
 
-#if 1
+#if 0
 	psImage *save = psImageCopy(NULL, sumMatrix, PS_TYPE_F32);
         psFitsWriteImageSimple ("sumMatrix.fits", save, NULL);
@@ -1115,7 +1115,9 @@
 	}
 
+# if (0)
         for (int i = 0; i < solution->n; i++) {
 	    psLogMsg("psModules.imcombine", PS_LOG_DETAIL, "Single solution %d: %lf\n", i, solution->data.F64[i]);
         }
+# endif
 
         if (!kernels->solution1) {
@@ -1169,5 +1171,5 @@
         }
 
-#if 1
+#if 0
 	psImage *save = psImageCopy(NULL, sumMatrix, PS_TYPE_F32);
         psFitsWriteImageSimple ("sumMatrix.fits", save, NULL);
@@ -1187,5 +1189,5 @@
 	solution = psMatrixSolveSVD(solution, sumMatrix, sumVector, NAN);
 
-#if (1)
+#if (0)
         for (int i = 0; i < solution->n; i++) {
             fprintf(stderr, "Dual solution %d: %lf\n", i, solution->data.F64[i]);
Index: /branches/eam_branches/ipp-20100823/psModules/src/imcombine/pmSubtractionKernels.c
===================================================================
--- /branches/eam_branches/ipp-20100823/psModules/src/imcombine/pmSubtractionKernels.c	(revision 29505)
+++ /branches/eam_branches/ipp-20100823/psModules/src/imcombine/pmSubtractionKernels.c	(revision 29506)
@@ -183,7 +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,
@@ -290,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,