Index: trunk/psModules/src/imcombine/pmSubtractionEquation.c
===================================================================
--- trunk/psModules/src/imcombine/pmSubtractionEquation.c	(revision 19532)
+++ trunk/psModules/src/imcombine/pmSubtractionEquation.c	(revision 19765)
@@ -758,4 +758,25 @@
         calculatePenalty(sumVector, kernels);
 
+#ifdef TESTING
+        {
+            psImage *inverse = psMatrixInvert(NULL, sumMatrix, NULL);
+            psFits *fits = psFitsOpen("matrixInv.fits", "w");
+            psFitsWriteImage(fits, NULL, inverse, 0, NULL);
+            psFitsClose(fits);
+            psFree(inverse);
+        }
+        {
+            psImage *X = psMatrixInvert(NULL, sumMatrix, NULL);
+            psImage *Xt = psMatrixTranspose(NULL, X);
+            psImage *XtX = psMatrixMultiply(NULL, Xt, X);
+            psFits *fits = psFitsOpen("matrixErr.fits", "w");
+            psFitsWriteImage(fits, NULL, XtX, 0, NULL);
+            psFitsClose(fits);
+            psFree(X);
+            psFree(Xt);
+            psFree(XtX);
+        }
+#endif
+
         psVector *permutation = NULL;       // Permutation vector, required for LU decomposition
         psImage *luMatrix = psMatrixLUD(NULL, &permutation, sumMatrix);
@@ -769,4 +790,5 @@
         }
         kernels->solution1 = psMatrixLUSolve(kernels->solution1, luMatrix, sumVector, permutation);
+
         psFree(sumVector);
         psFree(luMatrix);
@@ -806,13 +828,4 @@
         calculatePenalty(sumVector2, kernels);
 
-#if 0
-        // Apply weighting to maximise the difference between solution coefficients for the same functions
-        double fudge = PS_SQR(2 * stamps->footprint + 1); // Fudge factor
-        for (int i = 0; i < kernels->num; i++) {
-            sumMatrix1->data.F64[i][i] -= fudge;
-            sumMatrix2->data.F64[i][i] += fudge;
-        }
-#endif
-
         // Pure matrix operations
 
@@ -941,4 +954,5 @@
 #endif
 
+#ifdef TESTING
         {
             psFits *fits = psFitsOpen("sumMatrix1.fits", "w");
@@ -961,5 +975,5 @@
             psFitsClose(fits);
         }
-
+#endif
 
         kernels->solution1 = a;