Index: /branches/pap/psModules/src/imcombine/pmSubtraction.c
===================================================================
--- /branches/pap/psModules/src/imcombine/pmSubtraction.c	(revision 25897)
+++ /branches/pap/psModules/src/imcombine/pmSubtraction.c	(revision 25898)
@@ -65,4 +65,40 @@
     return out;
 }
+
+// Contribute to an image of the solved kernel component for ISIS
+static void solvedKernelISIS(psKernel *kernel, // Kernel, updated
+                             const pmSubtractionKernels *kernels, // Kernel basis functions
+                             float value,                         // Normalisation value for basis function
+                             int index                  // Index of basis function of interest
+    )
+{
+    int size = kernels->size;           // Kernel half-size
+    psArray *preCalc = kernels->preCalc->data[index]; // Precalculated values
+#if 0
+    psVector *xKernel = preCalc->data[0]; // Kernel in x
+    psVector *yKernel = preCalc->data[1]; // Kernel in y
+    // Iterating over the kernel
+    for (int y = 0, v = -size; v <= size; y++, v++) {
+        float yValue = value * yKernel->data.F32[y];
+        for (int x = 0, u = -size; u <= size; x++, u++) {
+            kernel->kernel[v][u] +=  yValue * xKernel->data.F32[x];
+        }
+    }
+    // Photometric scaling for even kernels only
+    if (kernels->u->data.S32[i] % 2 == 0 && kernels->v->data.S32[i] % 2 == 0) {
+        kernel->kernel[0][0] -= value;
+    }
+#else
+    psKernel *k = preCalc->data[2]; // Kernel image
+    for (int v = -size; v <= size; v++) {
+        for (int u = -size; u <= size; u++) {
+            kernel->kernel[v][u] +=  value * k->kernel[v][u];
+        }
+    }
+#endif
+
+    return;
+}
+
 
 // Generate an image of the solved kernel
@@ -116,13 +152,5 @@
           case PM_SUBTRACTION_KERNEL_GUNK: {
               if (i < kernels->inner) {
-                  // Using pre-calculated function
-                  psKernel *preCalc = kernels->preCalc->data[i]; // Precalculated values
-                  // Iterating over the kernel
-                  for (int v = -size; v <= size; v++) {
-                      for (int u = -size; u <= size; u++) {
-                          kernel->kernel[v][u] += preCalc->kernel[v][u] * value;
-                          // Photometric scaling is built into the preCalc kernel --- no subtraction!
-                      }
-                  }
+                  solvedKernelISIS(kernel, kernels, value, i);
               } else {
                   // Using delta function
@@ -135,27 +163,5 @@
           }
           case PM_SUBTRACTION_KERNEL_ISIS: {
-              psArray *preCalc = kernels->preCalc->data[i]; // Precalculated values
-#if 0
-              psVector *xKernel = preCalc->data[0]; // Kernel in x
-              psVector *yKernel = preCalc->data[1]; // Kernel in y
-              // Iterating over the kernel
-              for (int y = 0, v = -size; v <= size; y++, v++) {
-                  float yValue = value * yKernel->data.F32[y];
-                  for (int x = 0, u = -size; u <= size; x++, u++) {
-                      kernel->kernel[v][u] +=  yValue * xKernel->data.F32[x];
-                  }
-              }
-              // Photometric scaling for even kernels only
-              if (kernels->u->data.S32[i] % 2 == 0 && kernels->v->data.S32[i] % 2 == 0) {
-                  kernel->kernel[0][0] -= value;
-              }
-#else
-              psKernel *k = preCalc->data[2]; // Kernel image
-              for (int v = -size; v <= size; v++) {
-                  for (int u = -size; u <= size; u++) {
-                      kernel->kernel[v][u] +=  value * k->kernel[v][u];
-                  }
-              }
-#endif
+              solvedKernelISIS(kernel, kernels, value, i);
               break;
           }
@@ -446,4 +452,63 @@
 
     return sys;
+}
+
+// Convolve a stamp using an ISIS kernel basis function
+static psKernel *convolveStampISIS(const psKernel *image, // Image to convolve
+                                   const pmSubtractionKernels *kernels, // Kernel basis functions
+                                   int index,                            // Index of basis function of interest
+                                   int footprint                         // Half-size of stamp
+    )
+{
+    psArray *preCalc = kernels->preCalc->data[index]; // Precalculated data
+#if 1
+    // Convolving using separable convolution
+    psVector *xKernel = preCalc->data[0]; // Kernel in x
+    psVector *yKernel = preCalc->data[1]; // Kernel in y
+    int size = kernels->size;     // Size of kernel
+
+    // Convolve in x
+    // Need to convolve a bit more than the footprint, for the y convolution
+    int yMin = -size - footprint, yMax = size + footprint; // Range for y
+    psKernel *temp = psKernelAlloc(yMin, yMax,
+                                   -footprint, footprint); // Temporary convolution; NOTE: wrong way!
+    for (int y = yMin; y <= yMax; y++) {
+        for (int x = -footprint; x <= footprint; x++) {
+            float value = 0.0;    // Value of convolved pixel
+            int uMin = x - size, uMax = x + size; // Range for u
+            psF32 *xKernelData = &xKernel->data.F32[xKernel->n - 1]; // Kernel values
+            psF32 *imageData = &image->kernel[y][uMin]; // Image values
+            for (int u = uMin; u <= uMax; u++, xKernelData--, imageData++) {
+                value += *xKernelData * *imageData;
+            }
+            temp->kernel[x][y] = value; // NOTE: putting in wrong way!
+        }
+    }
+
+    // Convolve in y
+    psKernel *convolved = psKernelAlloc(-footprint, footprint, -footprint, footprint);// Convolved image
+    for (int x = -footprint; x <= footprint; x++) {
+        for (int y = -footprint; y <= footprint; y++) {
+            float value = 0.0;    // Value of convolved pixel
+            int vMin = y - size, vMax = y + size; // Range for v
+            psF32 *yKernelData = &yKernel->data.F32[yKernel->n - 1]; // Kernel values
+            psF32 *imageData = &temp->kernel[x][vMin]; // Image values; NOTE: wrong way!
+            for (int v = vMin; v <= vMax; v++, yKernelData--, imageData++) {
+                value += *yKernelData * *imageData;
+            }
+            convolved->kernel[y][x] = value;
+        }
+    }
+    psFree(temp);
+
+    // Photometric scaling for even kernels only
+    if (kernels->u->data.S32[index] % 2 == 0 && kernels->v->data.S32[index] % 2 == 0) {
+        convolveSub(convolved, image, footprint);
+    }
+    return convolved;
+#else
+    // Convolving using precalculated kernel
+    return p_pmSubtractionConvolveStampPrecalc(image, preCalc->data[2]);
+#endif
 }
 
@@ -599,5 +664,5 @@
           if (index < kernels->inner) {
               // Photometric scaling is already built in to the precalculated kernel
-              return p_pmSubtractionConvolveStampPrecalc(image, kernels->preCalc->data[index]);
+              return convolveStampISIS(image, kernels, index, footprint);
           }
           // Using delta function
@@ -609,53 +674,5 @@
       }
       case PM_SUBTRACTION_KERNEL_ISIS: {
-          psArray *preCalc = kernels->preCalc->data[index]; // Precalculated values
-#if 1
-          psVector *xKernel = preCalc->data[0]; // Kernel in x
-          psVector *yKernel = preCalc->data[1]; // Kernel in y
-          int size = kernels->size;     // Size of kernel
-
-          // Convolve in x
-          // Need to convolve a bit more than the footprint, for the y convolution
-          int yMin = -size - footprint, yMax = size + footprint; // Range for y
-          psKernel *temp = psKernelAlloc(yMin, yMax,
-                                         -footprint, footprint); // Temporary convolution; NOTE: wrong way!
-          for (int y = yMin; y <= yMax; y++) {
-              for (int x = -footprint; x <= footprint; x++) {
-                  float value = 0.0;    // Value of convolved pixel
-                  int uMin = x - size, uMax = x + size; // Range for u
-                  psF32 *xKernelData = &xKernel->data.F32[xKernel->n - 1]; // Kernel values
-                  psF32 *imageData = &image->kernel[y][uMin]; // Image values
-                  for (int u = uMin; u <= uMax; u++, xKernelData--, imageData++) {
-                      value += *xKernelData * *imageData;
-                  }
-                  temp->kernel[x][y] = value; // NOTE: putting in wrong way!
-              }
-          }
-
-          // Convolve in y
-          psKernel *convolved = psKernelAlloc(-footprint, footprint, -footprint, footprint);// Convolved image
-          for (int x = -footprint; x <= footprint; x++) {
-              for (int y = -footprint; y <= footprint; y++) {
-                  float value = 0.0;    // Value of convolved pixel
-                  int vMin = y - size, vMax = y + size; // Range for v
-                  psF32 *yKernelData = &yKernel->data.F32[yKernel->n - 1]; // Kernel values
-                  psF32 *imageData = &temp->kernel[x][vMin]; // Image values; NOTE: wrong way!
-                  for (int v = vMin; v <= vMax; v++, yKernelData--, imageData++) {
-                      value += *yKernelData * *imageData;
-                  }
-                  convolved->kernel[y][x] = value;
-              }
-          }
-          psFree(temp);
-
-          // Photometric scaling for even kernels only
-          if (kernels->u->data.S32[index] % 2 == 0 && kernels->v->data.S32[index] % 2 == 0) {
-              convolveSub(convolved, image, footprint);
-          }
-          return convolved;
-#else
-          return p_pmSubtractionConvolveStampPrecalc(image, preCalc->data[2]);
-#endif
-
+          return convolveStampISIS(image, kernels, index, footprint);
       }
       case PM_SUBTRACTION_KERNEL_RINGS: {
Index: /branches/pap/psModules/src/imcombine/pmSubtractionKernels.c
===================================================================
--- /branches/pap/psModules/src/imcombine/pmSubtractionKernels.c	(revision 25897)
+++ /branches/pap/psModules/src/imcombine/pmSubtractionKernels.c	(revision 25898)
@@ -81,4 +81,5 @@
     kernels->penalties = psVectorRealloc(kernels->penalties, start + numNew);
     kernels->inner = start;
+    kernels->num += numNew;
 
     // Generate a set of kernels for each (u,v)
@@ -99,4 +100,10 @@
         }
     }
+
+    kernels->widths->n = start + numNew;
+    kernels->u->n = start + numNew;
+    kernels->v->n = start + numNew;
+    kernels->preCalc->n = start + numNew;
+    kernels->penalties->n = start + numNew;
 
     return true;
@@ -471,8 +478,7 @@
     PS_ASSERT_INT_LESS_THAN(inner, size, NULL);
 
-    // XXX GUNK doesn't seem to work --- doesn't add the POIS components, or at least, they're not noticed
-
     pmSubtractionKernels *kernels = p_pmSubtractionKernelsRawISIS(size, spatialOrder, fwhms, orders,
                                                                   penalty, mode); // Kernels
+    kernels->type = PM_SUBTRACTION_KERNEL_GUNK;
     psStringPrepend(&kernels->description, "GUNK=");
     psStringAppend(&kernels->description, "+POIS(%d,%d)", inner, spatialOrder);