Index: trunk/psModules/src/imcombine/pmSubtractionKernels.c
===================================================================
--- trunk/psModules/src/imcombine/pmSubtractionKernels.c	(revision 25120)
+++ trunk/psModules/src/imcombine/pmSubtractionKernels.c	(revision 25999)
@@ -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)
@@ -94,9 +95,15 @@
             kernels->v->data.S32[index] = v;
             kernels->preCalc->data[index] = NULL;
-            kernels->penalties->data.F32[index] = kernels->penalty * (PS_SQR(u) + PS_SQR(v));
+            kernels->penalties->data.F32[index] = kernels->penalty * PS_SQR(PS_SQR(u) + PS_SQR(v));
 
             psTrace("psModules.imcombine", 7, "Kernel %d: %d %d\n", index, u, v);
         }
     }
+
+    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;
@@ -140,7 +147,8 @@
         for (int uOrder = 0; uOrder <= orders->data.S32[i]; uOrder++) {
             for (int vOrder = 0; vOrder <= orders->data.S32[i] - uOrder; vOrder++, index++) {
-                psArray *preCalc = psArrayAlloc(2); // Array to hold precalculated values
+                psArray *preCalc = psArrayAlloc(3); // Array to hold precalculated values
                 psVector *xKernel = preCalc->data[0] = subtractionKernelISIS(sigma, uOrder, size); // x Kernel
                 psVector *yKernel = preCalc->data[1] = subtractionKernelISIS(sigma, vOrder, size); // y Kernel
+                psKernel *kernel = preCalc->data[2] = psKernelAlloc(-size, size, -size, size);      // Kernel
 
                 // Calculate moments
@@ -149,5 +157,6 @@
                     for (int u = -size, x = 0; u <= size; u++, x++) {
                         double value = xKernel->data.F32[x] * yKernel->data.F32[y]; // Value of kernel
-                        moment += value * (PS_SQR(u) + PS_SQR(v));
+                        kernel->kernel[v][u] = value;
+                        moment += value * PS_SQR((PS_SQR(u) + PS_SQR(v)));
                     }
                 }
@@ -164,6 +173,19 @@
                     psBinaryOp(xKernel, xKernel, "*", psScalarAlloc(sum, PS_TYPE_F32));
                     psBinaryOp(yKernel, yKernel, "*", psScalarAlloc(sum, PS_TYPE_F32));
+                    psBinaryOp(kernel->image, kernel->image, "*", psScalarAlloc(PS_SQR(sum), PS_TYPE_F32));
+                    kernel->kernel[0][0] -= 1.0;
                     moment *= PS_SQR(sum);
                 }
+
+
+#if 0
+                double sum = 0.0;   // Sum of kernel component
+                for (int v = -size; v <= size; v++) {
+                    for (int u = -size; u <= size; u++) {
+                        sum += kernel->kernel[v][u];
+                    }
+                }
+                fprintf(stderr, "%d sum: %lf\n", index, sum);
+#endif
 
                 kernels->widths->data.F32[index] = fwhms->data.F32[i];
@@ -456,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);
@@ -577,5 +598,5 @@
                                     poly->data.F32[j] = polyVal;
                                     norm += polyVal;
-                                    moment += polyVal * (PS_SQR(u) + PS_SQR(v));
+                                    moment += polyVal * PS_SQR(PS_SQR(u) + PS_SQR(v));
 
                                     psVectorExtend(uCoords, RINGS_BUFFER, 1);
@@ -603,5 +624,5 @@
                         psBinaryOp(poly, poly, "*", psScalarAlloc(1.0 / norm, PS_TYPE_F32));
                     }
-//                    moment /= norm;
+                    moment /= norm;
                 }