Index: /branches/eam_branches/20091201/psModules/src/imcombine/pmSubtractionEquation.c
===================================================================
--- /branches/eam_branches/20091201/psModules/src/imcombine/pmSubtractionEquation.c	(revision 26339)
+++ /branches/eam_branches/20091201/psModules/src/imcombine/pmSubtractionEquation.c	(revision 26340)
@@ -104,5 +104,5 @@
 
 	    // Spatial variation of kernel coeffs
-	    if (mode | PM_SUBTRACTION_EQUATION_KERNELS) {
+	    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) {
@@ -147,14 +147,21 @@
 	    double normTerm = sumRC * poly[iTerm];
 	    double bgTerm = sumC * poly[iTerm];
-	    if ((mode | PM_SUBTRACTION_EQUATION_NORM) && (mode | PM_SUBTRACTION_EQUATION_KERNELS)) {
+	    if ((mode & PM_SUBTRACTION_EQUATION_NORM) && (mode & PM_SUBTRACTION_EQUATION_KERNELS)) {
 		matrix->data.F64[iIndex][normIndex] = normTerm;
 		matrix->data.F64[normIndex][iIndex] = normTerm;
 	    }
-	    if ((mode | PM_SUBTRACTION_EQUATION_BG) && (mode | PM_SUBTRACTION_EQUATION_KERNELS)) {
+	    if ((mode & PM_SUBTRACTION_EQUATION_BG) && (mode & PM_SUBTRACTION_EQUATION_KERNELS)) {
 		matrix->data.F64[iIndex][bgIndex] = bgTerm;
 		matrix->data.F64[bgIndex][iIndex] = bgTerm;
 	    }
-	    if (mode | PM_SUBTRACTION_EQUATION_KERNELS) {
+	    if (mode & PM_SUBTRACTION_EQUATION_KERNELS) {
 		vector->data.F64[iIndex] = sumIC * poly[iTerm];
+		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 * normTerm;
+		}
 	    }
 	}
@@ -196,13 +203,14 @@
         }
     }
-    if (mode | PM_SUBTRACTION_EQUATION_NORM) {
+    if (mode & PM_SUBTRACTION_EQUATION_NORM) {
 	matrix->data.F64[normIndex][normIndex] = sumRR;
 	vector->data.F64[normIndex] = sumIR;
-    }
-    if (mode | PM_SUBTRACTION_EQUATION_BG) {
+	// subtract sum over kernels * kernel solution
+    }
+    if (mode & PM_SUBTRACTION_EQUATION_BG) {
 	matrix->data.F64[bgIndex][bgIndex] = sum1;
 	vector->data.F64[bgIndex] = sumI;
     }
-    if ((mode | PM_SUBTRACTION_EQUATION_NORM) && (mode | PM_SUBTRACTION_EQUATION_BG)) {
+    if ((mode & PM_SUBTRACTION_EQUATION_NORM) && (mode & PM_SUBTRACTION_EQUATION_BG)) {
 	matrix->data.F64[normIndex][bgIndex] = sumR;
 	matrix->data.F64[bgIndex][normIndex] = sumR;
@@ -953,4 +961,16 @@
 #endif
 
+	// XXX test: save the matrix A and vector b:
+	{
+	    psFits *fits = psFitsOpen("matrix.fits", "w");
+            psFitsWriteImage(fits, NULL, sumMatrix, 0, NULL);
+            psFitsClose(fits);
+
+	    FILE *f = fopen ("vector.dat", "w");
+	    int fd = fileno(f);
+	    p_psVectorPrint (fd, sumVector, "B");
+	    fclose (f);
+	}
+
         psVector *permutation = NULL;       // Permutation vector, required for LU decomposition
         psImage *luMatrix = psMatrixLUDecomposition(NULL, &permutation, sumMatrix);
@@ -979,13 +999,13 @@
 
 	// only update the solutions that we chose to calculate:
-	if (mode | PM_SUBTRACTION_EQUATION_NORM) {
+	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) {
+	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) {
+	if (mode & PM_SUBTRACTION_EQUATION_KERNELS) {
 	    int numKernels = kernels->num;
 	    int spatialOrder = kernels->spatialOrder;       // Order of spatial variation
@@ -1225,5 +1245,5 @@
      }
 
-    pmSubtractionVisualPlotLeastSquares((pmSubtractionStampList *) stamps); //casting away const
+    // pmSubtractionVisualPlotLeastSquares((pmSubtractionStampList *) stamps); //casting away const
     return true;
 }
@@ -1241,5 +1261,4 @@
     double devNorm = 1.0 / (double)numPixels; // Normalisation for deviations
     int numKernels = kernels->num;      // Number of kernels
-    double norm = NAN;
 
     psImage *polyValues = NULL;         // Polynomial values
@@ -1258,5 +1277,5 @@
         // Calculate coefficients of the kernel basis functions
         polyValues = p_pmSubtractionPolynomial(polyValues, kernels->spatialOrder, stamp->xNorm, stamp->yNorm);
-        norm = p_pmSubtractionSolutionNorm(kernels); // Normalisation
+        double norm = p_pmSubtractionSolutionNorm(kernels); // Normalisation
         double background = p_pmSubtractionSolutionBackground(kernels, polyValues);// Difference in background
 
@@ -1404,11 +1423,19 @@
     
     }
+
+    // calculate and report the normalization and background for the image center
+    { 
+	polyValues = p_pmSubtractionPolynomial(polyValues, kernels->spatialOrder, 0.0, 0.0);
+	double norm = p_pmSubtractionSolutionNorm(kernels); // Normalisation
+	double background = p_pmSubtractionSolutionBackground(kernels, polyValues);// Difference in background
+	psLogMsg("psModules.imcombine", PS_LOG_INFO, "normalization: %f, background: %f", norm, background);
+    }
+
+    pmSubtractionVisualShowFit();
+    pmSubtractionVisualPlotFit(kernels);
+
     psFree(residual);
     psFree(polyValues);
 
-    psLogMsg("psModules.imcombine", PS_LOG_INFO, "normalization: %f", norm);
-    pmSubtractionVisualShowFit();
-    pmSubtractionVisualPlotFit(kernels);
-
     return deviations;
 }