Index: /branches/eam_branches/ipp-20100621/psModules/src/imcombine/pmSubtractionEquation.c
===================================================================
--- /branches/eam_branches/ipp-20100621/psModules/src/imcombine/pmSubtractionEquation.c	(revision 28893)
+++ /branches/eam_branches/ipp-20100621/psModules/src/imcombine/pmSubtractionEquation.c	(revision 28894)
@@ -603,5 +603,6 @@
                 matrix->data.F64[index][index] += norm * penalties1->data.F32[i];
                 if (kernels->mode == PM_SUBTRACTION_MODE_DUAL) {
-		    fprintf (stderr, "penalty: %f + %f (%f * %f)\n", matrix->data.F64[index + numParams + 2][index + numParams + 2], norm * penalties2->data.F32[i], norm, penalties2->data.F32[i]);
+		    fprintf (stderr, "penalty: (x^%d y^%d fwhm %f) : %f + %f (%f * %f)\n", kernels->u->data.S32[index], kernels->v->data.S32[index], kernels->widths->data.F32[index], 
+			     matrix->data.F64[index + numParams + 2][index + numParams + 2], norm * penalties2->data.F32[i], norm, penalties2->data.F32[i]);
 		    matrix->data.F64[index + numParams + 2][index + numParams + 2] += norm * penalties2->data.F32[i];			     
                     // matrix[i][i] is ~ (k_i * I_1)(k_i * I_1)
@@ -1116,5 +1117,5 @@
 
         int normIndex = PM_SUBTRACTION_INDEX_NORM(kernels); // Index for normalisation
-        calculatePenalty(sumMatrix, sumVector, kernels, sumMatrix->data.F64[normIndex][normIndex] / 10000.0);
+        calculatePenalty(sumMatrix, sumVector, kernels, sumMatrix->data.F64[normIndex][normIndex] / 100.0);
 #endif
 
Index: /branches/eam_branches/ipp-20100621/psModules/src/imcombine/pmSubtractionKernels.c
===================================================================
--- /branches/eam_branches/ipp-20100621/psModules/src/imcombine/pmSubtractionKernels.c	(revision 28893)
+++ /branches/eam_branches/ipp-20100621/psModules/src/imcombine/pmSubtractionKernels.c	(revision 28894)
@@ -218,5 +218,7 @@
         if (uOrder % 2 == 0 && vOrder % 2 == 0) {
             // Even functions: normalise to unit sum and subtract null pixel so that sum is zero
-            scale2D = 1.0 / fabs(sum);
+	    // Re-normalize 
+            // scale2D  = 1.0 / fabs(sum);
+            scale2D  = 1.0 / sqrt(sum2);
             zeroNull = true;
         } else {
@@ -252,19 +254,23 @@
 
     if (zeroNull) {
-        preCalc->kernel->kernel[0][0] -= 1.0;
-    }
-
-#if 0
+        // preCalc->kernel->kernel[0][0] -= 1.0;
+        preCalc->kernel->kernel[0][0] -= sum / sqrt (sum2);
+    }
+
+#if 1
     {
-        double sum = 0.0;   // Sum of kernel component
+        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++) {
-                sum += preCalc->kernel->kernel[v][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 mod: %lf, null: %f, min: %lf, max: %lf, scale: %f\n", index, sum, preCalc->kernel->kernel[0][0], min, max, scale2D);
+        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
Index: /branches/eam_branches/ipp-20100621/psModules/src/imcombine/pmSubtractionStamps.c
===================================================================
--- /branches/eam_branches/ipp-20100621/psModules/src/imcombine/pmSubtractionStamps.c	(revision 28893)
+++ /branches/eam_branches/ipp-20100621/psModules/src/imcombine/pmSubtractionStamps.c	(revision 28894)
@@ -891,16 +891,21 @@
     pmSubtractionGetFWHMs(&fwhm1, &fwhm2);
 
+    bool zeroNull = false;
+    int uOrder = kernels->u->data.S32[index];
+    int vOrder = kernels->v->data.S32[index];
+    if (uOrder % 2 == 0 && vOrder % 2 == 0) zeroNull = true;
+
     if (EMPIRICAL) {
 	psKernel *convolution1 = stamp->convolutions1->data[index];
-	penalty1 = pmSubtractionKernelPenaltySingle(convolution1);
+	penalty1 = pmSubtractionKernelPenaltySingle(convolution1, zeroNull);
 
 	psKernel *convolution2 = stamp->convolutions2->data[index];
-	penalty2 = pmSubtractionKernelPenaltySingle(convolution2);
+	penalty2 = pmSubtractionKernelPenaltySingle(convolution2, zeroNull);
     } else {
 	pmSubtractionKernelPreCalc *kernel = kernels->preCalc->data[index];
-	float M2 = pmSubtractionKernelPenaltySingle(kernel->kernel);
-
-	penalty1 = M2 * PS_SQR(fwhm1 / 2.35); // rescale the unconvolved second-moment by the image second moment 
-	penalty2 = M2 * PS_SQR(fwhm2 / 2.35); // rescale the unconvolved second-moment by the image second moment 
+	float M2 = pmSubtractionKernelPenaltySingle(kernel->kernel, zeroNull);
+
+	penalty1 = M2 + PS_SQR(fwhm1 / 2.35); // rescale the unconvolved second-moment by the image second moment 
+	penalty2 = M2 + PS_SQR(fwhm2 / 2.35); // rescale the unconvolved second-moment by the image second moment 
     }
     kernels->penalties1->data.F32[index] = kernels->penalty * penalty1;
@@ -915,5 +920,5 @@
 }
 
-float pmSubtractionKernelPenaltySingle(psKernel *kernel)
+float pmSubtractionKernelPenaltySingle(psKernel *kernel, bool zeroNull)
 {
     // Calculate moments
@@ -924,4 +929,7 @@
 	for (int u = kernel->xMin; u <= kernel->xMax; u++) {
             double value = kernel->kernel[v][u];
+	    if (false && zeroNull && (u == 0) && (v == 0)) {
+		value += 1.0;
+	    }
             double value2 = PS_SQR(value);
             sum += value;
@@ -934,6 +942,6 @@
     penalty *= 1.0 / sum2;
 
-    if (0) {
-	fprintf(stderr, "min: %lf, max: %lf, moment: %lf\n", min, max, penalty);
+    if (1) {
+	fprintf(stderr, "min: %lf, max: %lf, moment: %lf, flux^2: %lf\n", min, max, penalty, sum2);
 	// psTrace("psModules.imcombine", 7, "Kernel %d: %f %d %d %f\n", index, fwhm, uOrder, vOrder, penalty);
     }
Index: /branches/eam_branches/ipp-20100621/psModules/src/imcombine/pmSubtractionStamps.h
===================================================================
--- /branches/eam_branches/ipp-20100621/psModules/src/imcombine/pmSubtractionStamps.h	(revision 28893)
+++ /branches/eam_branches/ipp-20100621/psModules/src/imcombine/pmSubtractionStamps.h	(revision 28894)
@@ -200,5 +200,5 @@
 bool pmSubtractionKernelPenaltiesStamp(pmSubtractionStamp *stamp, pmSubtractionKernels *kernels);
 bool pmSubtractionKernelPenalties(pmSubtractionStamp *stamp, pmSubtractionKernels *kernels, int index);
-float pmSubtractionKernelPenaltySingle(psKernel *kernel);
+float pmSubtractionKernelPenaltySingle(psKernel *kernel, bool zeroNull);
 
 #endif
Index: /branches/eam_branches/ipp-20100621/psModules/src/imcombine/pmSubtractionVisual.c
===================================================================
--- /branches/eam_branches/ipp-20100621/psModules/src/imcombine/pmSubtractionVisual.c	(revision 28893)
+++ /branches/eam_branches/ipp-20100621/psModules/src/imcombine/pmSubtractionVisual.c	(revision 28894)
@@ -234,5 +234,5 @@
 bool pmSubtractionVisualShowKernels(pmSubtractionKernels *kernels) {
 
-    if (!pmVisualTestLevel("ppsub.kern.final", 1)) return true;
+    if (!pmVisualTestLevel("ppsub.kernels.final", 1)) return true;
 
     if (!pmVisualInitWindow (&kapa1, "PPSub:Images")) {