Index: trunk/psModules/src/imcombine/pmSubtractionKernels.c
===================================================================
--- trunk/psModules/src/imcombine/pmSubtractionKernels.c	(revision 14680)
+++ trunk/psModules/src/imcombine/pmSubtractionKernels.c	(revision 14690)
@@ -150,16 +150,4 @@
     kernels->subIndex = -1;
 
-    if (psTraceGetLevel("psModules.imcombine.kernel") >= 10) {
-        for (int i = 0; i < num; i++) {
-            psKernel *kernel = kernels->preCalc->data[i]; // Kernel of interest
-            psString kernelName = NULL;
-            psStringAppend(&kernelName, "kernel%03d.fits", i);
-            psFits *kernelFile = psFitsOpen(kernelName, "w");
-            psFree(kernelName);
-            psFitsWriteImage(kernelFile, NULL, kernel->image, 0, NULL);
-            psFitsClose(kernelFile);
-        }
-    }
-
     return kernels;
 }
@@ -228,8 +216,28 @@
 
     // Subtract a reference kernel from all the others to maintain flux scaling
-    psKernel *subtract = kernels->preCalc->data[0]; // Kernel to subtract from all others
-    for (int i = 1; i < kernels->num; i++) {
-        psKernel *kernel = kernels->preCalc->data[i]; // Kernel to subtract
-        psBinaryOp(kernel->image, kernel->image, "-", subtract->image);
+    int subIndex = 0;                   // Index of kernel to subtract from others
+    psKernel *subtract = kernels->preCalc->data[subIndex]; // Kernel to subtract from others
+    int numGaussians = fwhms->n;       // Number of Gaussians
+    for (int i = 0, index = 0; i < numGaussians; i++) {
+        for (int uOrder = 0; uOrder <= orders->data.S32[i]; uOrder++) {
+            for (int vOrder = 0; vOrder <= orders->data.S32[i] - uOrder; vOrder++, index++) {
+                if (uOrder % 2 == 0 && vOrder % 2 == 0 && index != subIndex) {
+                    psKernel *kernel = kernels->preCalc->data[index]; // Kernel to subtract
+                    psBinaryOp(kernel->image, kernel->image, "-", subtract->image);
+                }
+            }
+        }
+    }
+
+    if (psTraceGetLevel("psModules.imcombine.kernel") >= 10) {
+        for (int i = 0; i < kernels->num; i++) {
+            psKernel *kernel = kernels->preCalc->data[i]; // Kernel of interest
+            psString kernelName = NULL;
+            psStringAppend(&kernelName, "kernel%03d.fits", i);
+            psFits *kernelFile = psFitsOpen(kernelName, "w");
+            psFree(kernelName);
+            psFitsWriteImage(kernelFile, NULL, kernel->image, 0, NULL);
+            psFitsClose(kernelFile);
+        }
     }
 
@@ -443,10 +451,17 @@
 
     // Subtract unity from the kernels to maintain photometric flux scaling
-    for (int i = 0; i < kernels->num; i++) {
-        psKernel *kernel = kernels->preCalc->data[i]; // Kernel of interest
-        kernel->kernel[0][0] -= 1.0;
-    }
-
-    int numGaussians = kernels->num;    // Number of ISIS kernels
+    int numGaussians = fwhms->n;       // Number of Gaussians
+    for (int i = 0, index = 0; i < numGaussians; i++) {
+        for (int uOrder = 0; uOrder <= orders->data.S32[i]; uOrder++) {
+            for (int vOrder = 0; vOrder <= orders->data.S32[i] - uOrder; vOrder++, index++) {
+                if (uOrder % 2 == 0 && vOrder % 2 == 0) {
+                    psKernel *kernel = kernels->preCalc->data[index]; // Kernel to subtract
+                    kernel->kernel[0][0] -= 1.0;
+                }
+            }
+        }
+    }
+
+    int numISIS = kernels->num;         // Number of ISIS kernels
     int numInner = PS_SQR(2 * inner + 1); // Number of inner kernel elements
 
@@ -455,5 +470,5 @@
     }
 
-    kernels->subIndex = numInner/2 + numGaussians;
+    kernels->subIndex = numInner/2 + numISIS;
     assert(kernels->u->data.S32[kernels->subIndex] == 0 &&
            kernels->v->data.S32[kernels->subIndex] == 0);