Index: trunk/psModules/src/imcombine/pmSubtractionKernels.c
===================================================================
--- trunk/psModules/src/imcombine/pmSubtractionKernels.c	(revision 14199)
+++ trunk/psModules/src/imcombine/pmSubtractionKernels.c	(revision 14305)
@@ -58,4 +58,5 @@
     kernels->preCalc = NULL;
     kernels->size = size;
+    kernels->inner = 0;
     kernels->spatialOrder = spatialOrder;
 
@@ -141,12 +142,4 @@
                 }
 
-#if 0
-                // Subtract a particular kernel in order to preserve photometric calibration across image
-                if (spatialOrder > 0 && index != kernels->subIndex) {
-                    psKernel *subKernel = kernels->preCalc->data[kernels->subIndex]; // Kernel to subtract
-                    (void)psBinaryOp(preCalc->image, preCalc->image, "-", subKernel->image);
-                }
-#endif
-
                 // Iterate over spatial order.  This loop creates the terms for
                 // x^xOrder * y^yOrder  such that (xOrder+yOrder) <= spatialOrder.
@@ -390,4 +383,112 @@
 }
 
+// Grid United with Normal Kernel
+pmSubtractionKernels *pmSubtractionKernelsGUNK(int size, int spatialOrder, const psVector *sigmas,
+                                               const psVector *orders, int inner)
+{
+    PS_ASSERT_INT_POSITIVE(size, NULL);
+    PS_ASSERT_INT_NONNEGATIVE(spatialOrder, NULL);
+    PS_ASSERT_VECTOR_NON_NULL(sigmas, NULL);
+    PS_ASSERT_VECTOR_TYPE(sigmas, PS_TYPE_F32, NULL);
+    PS_ASSERT_VECTOR_NON_NULL(orders, NULL);
+    PS_ASSERT_VECTOR_TYPE(orders, PS_TYPE_S32, NULL);
+    PS_ASSERT_VECTORS_SIZE_EQUAL(sigmas, orders, NULL);
+    PS_ASSERT_INT_NONNEGATIVE(inner, NULL);
+    PS_ASSERT_INT_LESS_THAN(inner, size, NULL);
+
+    int numGaussians = sigmas->n;       // Number of Gaussians
+    int numGaussianVars = 0;            // Number of Gaussian variant functions in the kernel
+    for (int i = 0; i < numGaussians; i++) {
+        int gaussOrder = orders->data.S32[i]; // Polynomial order to apply to Gaussian
+        numGaussianVars += (gaussOrder + 1) * (gaussOrder + 2) / 2;
+    }
+
+    int numInner = 2 * PS_SQR(inner) + 1; // Number of inner kernel elements
+    int numSpatial = (spatialOrder + 1) * (spatialOrder + 2) / 2; // Number of spatial variations of a kernel
+
+    int num = (numGaussianVars + numInner) * numSpatial; // Total number of basis functions
+
+    pmSubtractionKernels *kernels = pmSubtractionKernelsAlloc(num, PM_SUBTRACTION_KERNEL_GUNK,
+                                                              size, spatialOrder); // The kernels
+    kernels->preCalc = psArrayAlloc(numGaussianVars * numSpatial);
+    kernels->inner = numGaussianVars * numSpatial;
+
+    // Set the Gaussian kernel parameters
+    for (int i = 0, index = 0; i < numGaussians; i++) {
+        // Iterate over (u,v) order
+        for (int uOrder = 0; uOrder <= orders->data.S32[i]; uOrder++) {
+            for (int vOrder = 0; vOrder <= orders->data.S32[i] - uOrder; vOrder++) {
+
+
+                // Set the pre-calculated kernel
+                psKernel *preCalc = psKernelAlloc(-size, size, -size, size);
+                for (int v = -size; v <= size; v++) {
+                    for (int u = -size; u <= size; u++) {
+                        preCalc->kernel[v][u] = power(u, uOrder) * power(v, vOrder) *
+                            expf(-0.5 * (PS_SQR(u) + PS_SQR(v)) / PS_SQR(sigmas->data.F32[i]));
+                    }
+                }
+
+                // Iterate over spatial order.  This loop creates the terms for
+                // x^xOrder * y^yOrder  such that (xOrder+yOrder) <= spatialOrder.
+                for (int xOrder = 0; xOrder <= spatialOrder; xOrder++) {
+                    for (int yOrder = 0; yOrder <= spatialOrder - xOrder; yOrder++, index++) {
+                        kernels->sigma->data.F32[index] = sigmas->data.F32[i];
+                        kernels->u->data.S32[index] = uOrder;
+                        kernels->v->data.S32[index] = vOrder;
+                        kernels->xOrder->data.S32[index] = xOrder;
+                        kernels->yOrder->data.S32[index] = yOrder;
+                        kernels->preCalc->data[index] = psMemIncrRefCounter(preCalc);
+
+                        psTrace("psModules.imcombine", 7, "Kernel %d: %f %d %d %d %d\n", index,
+                                sigmas->data.F32[i], uOrder, vOrder, xOrder, yOrder);
+                    }
+                }
+
+                psFree(preCalc);        // Drop reference
+            }
+        }
+    }
+
+
+    // Generate a grid set of kernels for each (u,v)
+    for (int v = - inner, index = kernels->inner; v <= inner; v++) {
+        for (int u = - inner; u <= inner; u++) {
+            // Iterate over spatial order.  This loop creates the terms for
+            // x^xOrder * y^yOrder  such that (xOrder+yOrder) <= spatialOrder.
+            for (int xOrder = 0; xOrder <= spatialOrder; xOrder++) {
+                for (int yOrder = 0; yOrder <= spatialOrder - xOrder; yOrder++, index++) {
+                    kernels->u->data.S32[index] = u;
+                    kernels->v->data.S32[index] = v;
+                    kernels->xOrder->data.S32[index] = xOrder;
+                    kernels->yOrder->data.S32[index] = yOrder;
+
+                    psTrace("psModules.imcombine", 7, "Kernel %d: %d %d %d %d\n", index,
+                            u, v, xOrder, yOrder);
+                }
+            }
+        }
+    }
+
+    kernels->subIndex = kernels->inner + numInner / 2;
+    assert(kernels->u->data.S32[kernels->subIndex] == 0 &&
+           kernels->v->data.S32[kernels->subIndex] == 0 &&
+           kernels->xOrder->data.S32[kernels->subIndex] == 0 &&
+           kernels->yOrder->data.S32[kernels->subIndex] == 0);
+
+    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;
+}
 
 pmSubtractionKernels *pmSubtractionKernelsGenerate(pmSubtractionKernelsType type, int size, int spatialOrder,
@@ -404,4 +505,6 @@
       case PM_SUBTRACTION_KERNEL_FRIES:
         return pmSubtractionKernelsFRIES(size, spatialOrder, inner);
+      case PM_SUBTRACTION_KERNEL_GUNK:
+        return pmSubtractionKernelsGUNK(size, spatialOrder, sigmas, orders, inner);
       default:
         psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Unknown kernel type: %x", type);