Index: branches/eam_branches/20091201/psModules/src/imcombine/pmSubtractionEquation.c
===================================================================
--- branches/eam_branches/20091201/psModules/src/imcombine/pmSubtractionEquation.c	(revision 26701)
+++ branches/eam_branches/20091201/psModules/src/imcombine/pmSubtractionEquation.c	(revision 26702)
@@ -156,7 +156,5 @@
             if (mode & PM_SUBTRACTION_EQUATION_KERNELS) {
                 vector->data.F64[iIndex] = sumIC * poly[iTerm];
-                // XXX TEST for Hermitians: do not calculate A - norm*B - \sum(k x B),
-                // instead, calculate A - \sum(k x B), with full hermitians
-                if (0 && !(mode & PM_SUBTRACTION_EQUATION_NORM)) {
+                if (!(mode & PM_SUBTRACTION_EQUATION_NORM)) {
                     // subtract norm * sumRC * poly[iTerm]
                     psAssert (kernels->solution1, "programming error: define solution first!");
@@ -232,5 +230,6 @@
                                       const pmSubtractionKernels *kernels, // Kernels
                                       const psImage *polyValues, // Spatial polynomial values
-                                      int footprint // (Half-)Size of stamp
+                                      int footprint, // (Half-)Size of stamp
+				      const pmSubtractionEquationCalculationMode mode
                                       )
 {
@@ -272,4 +271,12 @@
     }
 
+
+    // initialize the matrix and vector for NOP on all coeffs.  we only fill in the coeffs we
+    // choose to calculate
+    psImageInit(matrix, 0.0);
+    psVectorInit(vector, 1.0);
+    for (int i = 0; i < matrix->numCols; i++) {
+        matrix->data.F64[i][i] = 1.0;
+    }
 
     for (int i = 0; i < numKernels; i++) {
@@ -306,21 +313,23 @@
             }
 
-            // Spatial variation
-            for (int iTerm = 0, iIndex = i; iTerm < numPoly; iTerm++, iIndex += numKernels) {
-                for (int jTerm = 0, jIndex = j; jTerm < numPoly; jTerm++, jIndex += numKernels) {
-                    double aa = sumAA * poly2[iTerm][jTerm];
-                    double bb = sumBB * poly2[iTerm][jTerm];
-                    double ab = sumAB * poly2[iTerm][jTerm];
-
-                    matrix->data.F64[iIndex][jIndex] = aa;
-                    matrix->data.F64[jIndex][iIndex] = aa;
-
-                    matrix->data.F64[iIndex + numParams][jIndex + numParams] = bb;
-                    matrix->data.F64[jIndex + numParams][iIndex + numParams] = bb;
-
-                    matrix->data.F64[iIndex][jIndex + numParams] = ab;
-                    matrix->data.F64[jIndex + numParams][iIndex] = ab;
-                }
-            }
+            // Spatial variation of kernel coeffs
+            if (mode & PM_SUBTRACTION_EQUATION_KERNELS) {
+		for (int iTerm = 0, iIndex = i; iTerm < numPoly; iTerm++, iIndex += numKernels) {
+		    for (int jTerm = 0, jIndex = j; jTerm < numPoly; jTerm++, jIndex += numKernels) {
+			double aa = sumAA * poly2[iTerm][jTerm];
+			double bb = sumBB * poly2[iTerm][jTerm];
+			double ab = sumAB * poly2[iTerm][jTerm];
+
+			matrix->data.F64[iIndex][jIndex] = aa;
+			matrix->data.F64[jIndex][iIndex] = aa;
+
+			matrix->data.F64[iIndex + numParams][jIndex + numParams] = bb;
+			matrix->data.F64[jIndex + numParams][iIndex + numParams] = bb;
+
+			matrix->data.F64[iIndex][jIndex + numParams] = ab;
+			matrix->data.F64[jIndex + numParams][iIndex] = ab;
+		    }
+		}
+	    }
         }
         for (int j = 0; j < i; j++) {
@@ -340,12 +349,14 @@
             }
 
-            // Spatial variation
-            for (int iTerm = 0, iIndex = i; iTerm < numPoly; iTerm++, iIndex += numKernels) {
-                for (int jTerm = 0, jIndex = j; jTerm < numPoly; jTerm++, jIndex += numKernels) {
-                    double ab = sumAB * poly2[iTerm][jTerm];
-                    matrix->data.F64[iIndex][jIndex + numParams] = ab;
-                    matrix->data.F64[jIndex + numParams][iIndex] = ab;
-                }
-            }
+            // Spatial variation of kernel coeffs
+            if (mode & PM_SUBTRACTION_EQUATION_KERNELS) {
+		for (int iTerm = 0, iIndex = i; iTerm < numPoly; iTerm++, iIndex += numKernels) {
+		    for (int jTerm = 0, jIndex = j; jTerm < numPoly; jTerm++, jIndex += numKernels) {
+			double ab = sumAB * poly2[iTerm][jTerm];
+			matrix->data.F64[iIndex][jIndex + numParams] = ab;
+			matrix->data.F64[jIndex + numParams][iIndex] = ab;
+		    }
+		}
+	    }
         }
 
@@ -403,18 +414,32 @@
             double bi2 = sumBI2 * poly[iTerm];
             double ai1 = sumAI1 * poly[iTerm];
-            double a = sumA * poly[iTerm];
+            double a   = sumA * poly[iTerm];
             double bi1 = sumBI1 * poly[iTerm];
-            double b = sumB * poly[iTerm];
-
-            matrix->data.F64[iIndex][normIndex] = ai1;
-            matrix->data.F64[normIndex][iIndex] = ai1;
-            matrix->data.F64[iIndex][bgIndex] = a;
-            matrix->data.F64[bgIndex][iIndex] = a;
-            matrix->data.F64[iIndex + numParams][normIndex] = bi1;
-            matrix->data.F64[normIndex][iIndex + numParams] = bi1;
-            matrix->data.F64[iIndex + numParams][bgIndex] = b;
-            matrix->data.F64[bgIndex][iIndex + numParams] = b;
-            vector->data.F64[iIndex] = ai2;
-            vector->data.F64[iIndex + numParams] = bi2;
+            double b   = sumB * poly[iTerm];
+
+            if ((mode & PM_SUBTRACTION_EQUATION_NORM) && (mode & PM_SUBTRACTION_EQUATION_KERNELS)) {
+		matrix->data.F64[iIndex][normIndex] = ai1;
+		matrix->data.F64[normIndex][iIndex] = ai1;
+		matrix->data.F64[iIndex + numParams][normIndex] = bi1;
+		matrix->data.F64[normIndex][iIndex + numParams] = bi1;
+	    }
+            if ((mode & PM_SUBTRACTION_EQUATION_BG) && (mode & PM_SUBTRACTION_EQUATION_KERNELS)) {
+		matrix->data.F64[iIndex][bgIndex] = a;
+		matrix->data.F64[bgIndex][iIndex] = a;
+		matrix->data.F64[iIndex + numParams][bgIndex] = b;
+		matrix->data.F64[bgIndex][iIndex + numParams] = b;
+	    }
+            if (mode & PM_SUBTRACTION_EQUATION_KERNELS) {
+		vector->data.F64[iIndex] = ai2;
+		vector->data.F64[iIndex + numParams] = bi2;
+                if (!(mode & PM_SUBTRACTION_EQUATION_NORM)) {
+                    // subtract norm * sumRC * poly[iTerm]
+                    psAssert (kernels->solution1, "programming error: define solution first!");
+                    int normIndex = PM_SUBTRACTION_INDEX_NORM(kernels); // Index for normalisation
+                    double norm = fabs(kernels->solution1->data.F64[normIndex]);  // Normalisation
+                    vector->data.F64[iIndex] -= norm * ai1;
+                    vector->data.F64[iIndex + numParams] -= norm * bi1;
+                }
+	    }
         }
     }
@@ -457,11 +482,16 @@
         }
     }
-    matrix->data.F64[bgIndex][normIndex] = sumI1;
-    matrix->data.F64[normIndex][bgIndex] = sumI1;
-    matrix->data.F64[normIndex][normIndex] = sumI1I1;
-    matrix->data.F64[bgIndex][bgIndex] = sum1;
-    vector->data.F64[bgIndex] = sumI2;
-    vector->data.F64[normIndex] = sumI1I2;
-
+    if (mode & PM_SUBTRACTION_EQUATION_NORM) {
+	matrix->data.F64[normIndex][normIndex] = sumI1I1;
+	vector->data.F64[normIndex] = sumI1I2;
+    }
+    if (mode & PM_SUBTRACTION_EQUATION_BG) {
+	matrix->data.F64[bgIndex][bgIndex] = sum1;
+	vector->data.F64[bgIndex] = sumI2;
+    }
+    if ((mode & PM_SUBTRACTION_EQUATION_NORM) && (mode & PM_SUBTRACTION_EQUATION_BG)) {
+	matrix->data.F64[bgIndex][normIndex] = sumI1;
+	matrix->data.F64[normIndex][bgIndex] = sumI1;
+    }
     return true;
 }
@@ -680,5 +710,5 @@
         status = calculateDualMatrixVector(stamp->matrix, stamp->vector,stamp->image1, stamp->image2,
                                            weight, window, stamp->convolutions1, stamp->convolutions2,
-                                           kernels, polyValues, footprint);
+                                           kernels, polyValues, footprint, mode);
         break;
       default:
@@ -1298,8 +1328,29 @@
         if (!kernels->solution1) {
             kernels->solution1 = psVectorAlloc(numSolution1, PS_TYPE_F64);
+            psVectorInit (kernels->solution1, 0.0);
         }
         if (!kernels->solution2) {
             kernels->solution2 = psVectorAlloc(numSolution2, PS_TYPE_F64);
-        }
+            psVectorInit (kernels->solution2, 0.0);
+        }
+
+        // only update the solutions that we chose to calculate:
+        if (mode & PM_SUBTRACTION_EQUATION_NORM) {
+            int normIndex = PM_SUBTRACTION_INDEX_NORM(kernels); // Index for normalisation
+            kernels->solution1->data.F64[normIndex] = solution->data.F64[normIndex];
+        }
+        if (mode & PM_SUBTRACTION_EQUATION_BG) {
+            int bgIndex = PM_SUBTRACTION_INDEX_BG(kernels); // Index in matrix for background
+            kernels->solution1->data.F64[bgIndex] = solution->data.F64[bgIndex];
+        }
+        if (mode & PM_SUBTRACTION_EQUATION_KERNELS) {
+            int numKernels = kernels->num;
+            for (int i = 0; i < numKernels * numSpatial; i++) {
+		// XXX fprintf (stderr, "keep\n");
+		kernels->solution1->data.F64[i] = solution->data.F64[i];
+		kernels->solution2->data.F64[i] = solution->data.F64[i + numSolution1];
+            }
+        }
+
 
         memcpy(kernels->solution1->data.F64, solution->data.F64,
Index: branches/eam_branches/20091201/psModules/src/imcombine/pmSubtractionMatch.c
===================================================================
--- branches/eam_branches/20091201/psModules/src/imcombine/pmSubtractionMatch.c	(revision 26701)
+++ branches/eam_branches/20091201/psModules/src/imcombine/pmSubtractionMatch.c	(revision 26702)
@@ -28,10 +28,9 @@
 static bool useFFT = true;              // Do convolutions using FFT
 
-# define SEPARATE 0
+# define SEPARATE 1
 # if (SEPARATE)
 # define SUBMODE PM_SUBTRACTION_EQUATION_NORM
 # else
 # define SUBMODE PM_SUBTRACTION_EQUATION_ALL
-// # define SUBMODE PM_SUBTRACTION_EQUATION_KERNELS
 # endif