Index: branches/eam_branches/20091201/psModules/src/imcombine/pmSubtractionEquation.c
===================================================================
--- branches/eam_branches/20091201/psModules/src/imcombine/pmSubtractionEquation.c	(revision 26575)
+++ branches/eam_branches/20091201/psModules/src/imcombine/pmSubtractionEquation.c	(revision 26576)
@@ -17,6 +17,6 @@
 // #define TESTING                         // TESTING output for debugging; may not work with threads!
 
-// #define USE_WEIGHT                      // Include weight (1/variance) in equation?
-// #define USE_WINDOW                      // Include weight (1/variance) in equation?
+#define USE_WEIGHT                      // Include weight (1/variance) in equation?
+#define USE_WINDOW                      // Include weight (1/variance) in equation?
 
 
@@ -36,5 +36,5 @@
                                   const psImage *polyValues, // Spatial polynomial values
                                   int footprint, // (Half-)Size of stamp
-                                  const pmSubtractionEquationCalculationMode mode
+				  const pmSubtractionEquationCalculationMode mode
                                   )
 {
@@ -68,10 +68,10 @@
     }
 
-    // initialize the matrix and vector for NOP on all coeffs.  we only fill in the coeffs we
+    // 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;
+	matrix->data.F64[i][i] = 1.0;
     }
 
@@ -84,87 +84,87 @@
 
     for (int i = 0; i < numKernels; i++) {
-        psKernel *iConv = convolutions->data[i]; // Convolution for index i
-        for (int j = i; j < numKernels; j++) {
-            psKernel *jConv = convolutions->data[j]; // Convolution for index j
-
-            double sumCC = 0.0;         // Sum of convolution products
-            for (int y = - footprint; y <= footprint; y++) {
-                for (int x = - footprint; x <= footprint; x++) {
-                    double cc = iConv->kernel[y][x] * jConv->kernel[y][x];
-                    if (weight) {
-                        cc *= weight->kernel[y][x];
-                    }
-                    if (window) {
-                        cc *= window->kernel[y][x];
-                    }
-                    sumCC += cc;
-                }
-            }
-
-            // 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 value = sumCC * poly2[iTerm][jTerm];
-                        matrix->data.F64[iIndex][jIndex] = value;
-                        matrix->data.F64[jIndex][iIndex] = value;
-                    }
-                }
-            }
-        }
-
-        double sumRC = 0.0;             // Sum of the reference-convolution products
-        double sumIC = 0.0;             // Sum of the input-convolution products
-        double sumC = 0.0;              // Sum of the convolution
-        for (int y = - footprint; y <= footprint; y++) {
-            for (int x = - footprint; x <= footprint; x++) {
-                float conv = iConv->kernel[y][x];
-                float in = input->kernel[y][x];
-                float ref = reference->kernel[y][x];
-                double ic = in * conv;
-                double rc = ref * conv;
-                double c = conv;
-                if (weight) {
-                    float wtVal = weight->kernel[y][x];
-                    ic *= wtVal;
-                    rc *= wtVal;
-                    c *= wtVal;
-                }
-                if (window) {
-                    float winVal = window->kernel[y][x];
-                    ic *= winVal;
-                    rc *= winVal;
-                    c  *= winVal;
-                }
-                sumIC += ic;
-                sumRC += rc;
-                sumC += c;
-            }
-        }
-        // Spatial variation
-        for (int iTerm = 0, iIndex = i; iTerm < numPoly; iTerm++, iIndex += numKernels) {
-            double normTerm = sumRC * poly[iTerm];
-            double bgTerm = sumC * poly[iTerm];
-            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)) {
-                matrix->data.F64[iIndex][bgIndex] = bgTerm;
-                matrix->data.F64[bgIndex][iIndex] = bgTerm;
-            }
-            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)) {
-                    // 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;
-                }
-            }
-        }
+	psKernel *iConv = convolutions->data[i]; // Convolution for index i
+	for (int j = i; j < numKernels; j++) {
+	    psKernel *jConv = convolutions->data[j]; // Convolution for index j
+
+	    double sumCC = 0.0;         // Sum of convolution products
+	    for (int y = - footprint; y <= footprint; y++) {
+		for (int x = - footprint; x <= footprint; x++) {
+		    double cc = iConv->kernel[y][x] * jConv->kernel[y][x];
+		    if (weight) {
+			cc *= weight->kernel[y][x];
+		    }
+		    if (window) {
+			cc *= window->kernel[y][x];
+		    }
+		    sumCC += cc;
+		}
+	    }
+
+	    // 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 value = sumCC * poly2[iTerm][jTerm];
+			matrix->data.F64[iIndex][jIndex] = value;
+			matrix->data.F64[jIndex][iIndex] = value;
+		    }
+		}
+	    }
+	}
+
+	double sumRC = 0.0;             // Sum of the reference-convolution products
+	double sumIC = 0.0;             // Sum of the input-convolution products
+	double sumC = 0.0;              // Sum of the convolution
+	for (int y = - footprint; y <= footprint; y++) {
+	    for (int x = - footprint; x <= footprint; x++) {
+		float conv = iConv->kernel[y][x];
+		float in = input->kernel[y][x];
+		float ref = reference->kernel[y][x];
+		double ic = in * conv;
+		double rc = ref * conv;
+		double c = conv;
+		if (weight) {
+		    float wtVal = weight->kernel[y][x];
+		    ic *= wtVal;
+		    rc *= wtVal;
+		    c *= wtVal;
+		}
+		if (window) {
+		    float winVal = window->kernel[y][x];
+		    ic *= winVal;
+		    rc *= winVal;
+		    c  *= winVal;
+		}
+		sumIC += ic;
+		sumRC += rc;
+		sumC += c;
+	    }
+	}
+	// Spatial variation
+	for (int iTerm = 0, iIndex = i; iTerm < numPoly; iTerm++, iIndex += numKernels) {
+	    double normTerm = sumRC * poly[iTerm];
+	    double bgTerm = sumC * poly[iTerm];
+	    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)) {
+		matrix->data.F64[iIndex][bgIndex] = bgTerm;
+		matrix->data.F64[bgIndex][iIndex] = bgTerm;
+	    }
+	    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)) {
+		    // 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;
+		}
+	    }
+	}
     }
 
@@ -181,20 +181,20 @@
             double rr = PS_SQR(ref);
             double one = 1.0;
-            if (weight) {
-                float wtVal = weight->kernel[y][x];
-                rr *= wtVal;
-                ir *= wtVal;
-                in *= wtVal;
-                ref *= wtVal;
-                one *= wtVal;
-            }
-            if (window) {
-                float  winVal = window->kernel[y][x];
-                rr      *= winVal;
-                ir      *= winVal;
-                in      *= winVal;
-                ref *= winVal;
-                one *= winVal;
-            }
+	    if (weight) {
+		float wtVal = weight->kernel[y][x];
+		rr *= wtVal;
+		ir *= wtVal;
+		in *= wtVal;
+		ref *= wtVal;
+		one *= wtVal;
+	    }
+	    if (window) {
+		float  winVal = window->kernel[y][x];
+		rr 	*= winVal;
+		ir 	*= winVal;
+		in 	*= winVal;
+		ref *= winVal;
+		one *= winVal;
+	    }
             sumRR += rr;
             sumIR += ir;
@@ -205,15 +205,15 @@
     }
     if (mode & PM_SUBTRACTION_EQUATION_NORM) {
-        matrix->data.F64[normIndex][normIndex] = sumRR;
-        vector->data.F64[normIndex] = sumIR;
-        // subtract sum over kernels * kernel solution
+	matrix->data.F64[normIndex][normIndex] = sumRR;
+	vector->data.F64[normIndex] = sumIR;
+	// subtract sum over kernels * kernel solution
     }
     if (mode & PM_SUBTRACTION_EQUATION_BG) {
-        matrix->data.F64[bgIndex][bgIndex] = sum1;
-        vector->data.F64[bgIndex] = sumI;
+	matrix->data.F64[bgIndex][bgIndex] = sum1;
+	vector->data.F64[bgIndex] = sumI;
     }
     if ((mode & PM_SUBTRACTION_EQUATION_NORM) && (mode & PM_SUBTRACTION_EQUATION_BG)) {
-        matrix->data.F64[normIndex][bgIndex] = sumR;
-        matrix->data.F64[bgIndex][normIndex] = sumR;
+	matrix->data.F64[normIndex][bgIndex] = sumR;
+	matrix->data.F64[bgIndex][normIndex] = sumR;
     }
     return true;
@@ -288,16 +288,16 @@
                     double bb = iConv2->kernel[y][x] * jConv2->kernel[y][x];
                     double ab = iConv1->kernel[y][x] * jConv2->kernel[y][x];
-                    if (weight) {
-                        float wtVal = weight->kernel[y][x];
-                        aa *= wtVal;
-                        bb *= wtVal;
-                        ab *= wtVal;
-                    }
-                    if (window) {
-                        float wtVal = window->kernel[y][x];
-                        aa *= wtVal;
-                        bb *= wtVal;
-                        ab *= wtVal;
-                    }
+		    if (weight) {
+			float wtVal = weight->kernel[y][x];
+			aa *= wtVal;
+			bb *= wtVal;
+			ab *= wtVal;
+		    }
+		    if (window) {
+			float wtVal = window->kernel[y][x];
+			aa *= wtVal;
+			bb *= wtVal;
+			ab *= wtVal;
+		    }
                     sumAA += aa;
                     sumBB += bb;
@@ -330,10 +330,10 @@
                 for (int x = - footprint; x <= footprint; x++) {
                     double ab = iConv1->kernel[y][x] * jConv2->kernel[y][x];
-                    if (weight) {
-                        ab *= weight->kernel[y][x];
-                    }
-                    if (window) {
-                        ab *= window->kernel[y][x];
-                    }
+		    if (weight) {
+			ab *= weight->kernel[y][x];
+		    }
+		    if (window) {
+			ab *= window->kernel[y][x];
+		    }
                     sumAB += ab;
                 }
@@ -369,24 +369,24 @@
                 double bi1 = b * i1;
 
-                if (weight) {
-                    float wtVal = weight->kernel[y][x];
-                    ai2 *= wtVal;
-                    bi2 *= wtVal;
-                    ai1 *= wtVal;
-                    bi1 *= wtVal;
-                    a *= wtVal;
-                    b *= wtVal;
-                    i2 *= wtVal;
-                }
-                if (window) {
-                    float wtVal = window->kernel[y][x];
-                    ai2 *= wtVal;
-                    bi2 *= wtVal;
-                    ai1 *= wtVal;
-                    bi1 *= wtVal;
-                    a *= wtVal;
-                    b *= wtVal;
-                    i2 *= wtVal;
-                }
+		if (weight) {
+		    float wtVal = weight->kernel[y][x];
+		    ai2 *= wtVal;
+		    bi2 *= wtVal;
+		    ai1 *= wtVal;
+		    bi1 *= wtVal;
+		    a *= wtVal;
+		    b *= wtVal;
+		    i2 *= wtVal;
+		}
+		if (window) {
+		    float wtVal = window->kernel[y][x];
+		    ai2 *= wtVal;
+		    bi2 *= wtVal;
+		    ai1 *= wtVal;
+		    bi1 *= wtVal;
+		    a *= wtVal;
+		    b *= wtVal;
+		    i2 *= wtVal;
+		}
                 sumAI2 += ai2;
                 sumBI2 += bi2;
@@ -434,20 +434,20 @@
             double i1i2 = i1 * i2;
 
-            if (weight) {
-                float wtVal = weight->kernel[y][x];
-                i1 *= wtVal;
-                i1i1 *= wtVal;
-                one *= wtVal;
-                i2 *= wtVal;
-                i1i2 *= wtVal;
-            }
-            if (window) {
-                float wtVal = window->kernel[y][x];
-                i1 *= wtVal;
-                i1i1 *= wtVal;
-                one *= wtVal;
-                i2 *= wtVal;
-                i1i2 *= wtVal;
-            }
+	    if (weight) {
+		float wtVal = weight->kernel[y][x];
+		i1 *= wtVal;
+		i1i1 *= wtVal;
+		one *= wtVal;
+		i2 *= wtVal;
+		i1i2 *= wtVal;
+	    }
+	    if (window) {
+		float wtVal = window->kernel[y][x];
+		i1 *= wtVal;
+		i1i1 *= wtVal;
+		one *= wtVal;
+		i2 *= wtVal;
+		i1i2 *= wtVal;
+	    }
             sumI1 += i1;
             sumI1I1 += i1i1;
@@ -486,5 +486,5 @@
             for (int xOrder = 0; xOrder <= spatialOrder - yOrder; xOrder++, index += numKernels) {
                 // Contribution to chi^2: a_i^2 P_i
-                psAssert(isfinite(penalties->data.F32[i]), "Invalid penalty");
+		psAssert(isfinite(penalties->data.F32[i]), "Invalid penalty");
                 matrix->data.F64[index][index] -= norm * penalties->data.F32[i];
             }
@@ -494,6 +494,5 @@
     return true;
 }
-#endif
-
+# endif
 
 //////////////////////////////////////////////////////////////////////////////////////////////////////////////
@@ -596,5 +595,5 @@
     int numBackground = PM_SUBTRACTION_POLYTERMS(kernels->bgOrder); // Number of background terms
 
-    // numKernels is the number of unique kernel images (one for each Gaussian modified by a specific polynomial).
+    // numKernels is the number of unique kernel images (one for each Gaussian modified by a specific polynomial).  
     // = \sum_i^N_Gaussians [(order + 1) * (order + 2) / 2], eg for 1 Gauss and 1st order, numKernels = 3
 
@@ -671,10 +670,10 @@
         status = calculateMatrixVector(stamp->matrix, stamp->vector, stamp->image2, stamp->image1,
                                        weight, window, stamp->convolutions1, kernels,
-                                       polyValues, footprint, mode);
+				       polyValues, footprint, mode);
         break;
       case PM_SUBTRACTION_MODE_2:
         status = calculateMatrixVector(stamp->matrix, stamp->vector, stamp->image1, stamp->image2,
                                        weight, window, stamp->convolutions2, kernels,
-                                       polyValues, footprint, mode);
+				       polyValues, footprint, mode);
         break;
       case PM_SUBTRACTION_MODE_DUAL:
@@ -737,10 +736,10 @@
         }
 
-        if ((stamp->x <= 0.0) && (stamp->y <= 0.0)) {
-            psAbort ("bad stamp");
-        }
-        if (!isfinite(stamp->x) && !isfinite(stamp->y)) {
-            psAbort ("bad stamp");
-        }
+	if ((stamp->x <= 0.0) && (stamp->y <= 0.0)) {
+	    psAbort ("bad stamp");
+	}
+	if (!isfinite(stamp->x) && !isfinite(stamp->y)) {
+	    psAbort ("bad stamp");
+	}
 
         if (pmSubtractionThreaded()) {
@@ -797,7 +796,7 @@
 double p_pmSubSolve_ChiSquare (pmSubtractionKernels *kernels, const pmSubtractionStampList *stamps);
 
-bool pmSubtractionSolveEquation(pmSubtractionKernels *kernels,
-                                const pmSubtractionStampList *stamps,
-                                const pmSubtractionEquationCalculationMode mode)
+bool pmSubtractionSolveEquation(pmSubtractionKernels *kernels, 
+				const pmSubtractionStampList *stamps, 
+				const pmSubtractionEquationCalculationMode mode)
 {
     PM_ASSERT_SUBTRACTION_KERNELS_NON_NULL(kernels, false);
@@ -914,261 +913,261 @@
 # define SVD_TOL 0.0
 
-        psVector *solution = NULL;
-        psVector *solutionErr = NULL;
+	psVector *solution = NULL;
+	psVector *solutionErr = NULL;
 
 #ifdef TESTING
-        psFitsWriteImageSimple("A.fits", sumMatrix, NULL);
-        psVectorWriteFile("B.dat", sumVector);
-#endif
-
-        // Use SVD to determine the kernel coeffs (and validate)
-        if (SVD_ANALYSIS) {
-
-            // We have sumVector and sumMatrix.  we are trying to solve the following equation:
-            // sumMatrix * x = sumVector.
-
-            // we can use any standard matrix inversion to solve this.  However, the basis
-            // functions in general have substantial correlation, so that the solution may be
-            // somewhat poorly determined or unstable.  If not numerically ill-conditioned, the
-            // system of equations may be statistically ill-conditioned.  Noise in the image
-            // will drive insignificant, but correlated, terms in the solution.  To avoid these
-            // problems, we can use SVD to identify numerically unconstrained values and to
-            // avoid statistically badly determined value.
-
-            // A = sumMatrix, B = sumVector
-            // SVD: A = U w V^T  -> A^{-1} = V (1/w) U^T
-            // x = V (1/w) (U^T B)
-            // \sigma_x = sqrt(diag(A^{-1}))
-            // solve for x and A^{-1} to get x & dx
-            // identify the elements of (1/w) that are nan (1/0.0) -> set to 0.0
-            // identify the elements of x that are insignificant (x / dx < 1.0? < 0.5?) -> set to 0.0
-
-            // If I use the SVD trick to re-condition the matrix, I need to break out the
-            // kernel and normalization terms from the background term.
-            // XXX is this true?  or was this due to an error in the analysis?
-
-            int bgIndex = PM_SUBTRACTION_INDEX_BG(kernels); // Index in matrix for background
-
-            // now pull out the kernel elements into their own square matrix
-            psImage  *kernelMatrix = psImageAlloc  (sumMatrix->numCols - 1, sumMatrix->numRows - 1, PS_TYPE_F64);
-            psVector *kernelVector = psVectorAlloc (sumMatrix->numCols - 1, PS_TYPE_F64);
-
-            for (int ix = 0, kx = 0; ix < sumMatrix->numCols; ix++) {
-                if (ix == bgIndex) continue;
-                for (int iy = 0, ky = 0; iy < sumMatrix->numRows; iy++) {
-                    if (iy == bgIndex) continue;
-                    kernelMatrix->data.F64[ky][kx] = sumMatrix->data.F64[iy][ix];
-                    ky++;
-                }
-                kernelVector->data.F64[kx] = sumVector->data.F64[ix];
-                kx++;
-            }
-
-            psImage *U = NULL;
-            psImage *V = NULL;
-            psVector *w = NULL;
-            if (!psMatrixSVD (&U, &w, &V, kernelMatrix)) {
-                psError(PS_ERR_UNKNOWN, false, "failed to perform SVD on sumMatrix\n");
-                return NULL;
-            }
-
-            // calculate A_inverse:
-            // Ainv = V * w * U^T
-            psImage *wUt  = p_pmSubSolve_wUt (w, U);
-            psImage *Ainv = p_pmSubSolve_VwUt (V, wUt);
-            psImage *Xvar = NULL;
-            psFree (wUt);
+	psFitsWriteImageSimple("A.fits", sumMatrix, NULL);
+	psVectorWriteFile ("B.dat", sumVector);
+#endif
+
+	// Use SVD to determine the kernel coeffs (and validate)
+	if (SVD_ANALYSIS) {
+
+	    // We have sumVector and sumMatrix.  we are trying to solve the following equation:
+	    // sumMatrix * x = sumVector.
+
+	    // we can use any standard matrix inversion to solve this.  However, the basis
+	    // functions in general have substantial correlation, so that the solution may be
+	    // somewhat poorly determined or unstable.  If not numerically ill-conditioned, the
+	    // system of equations may be statistically ill-conditioned.  Noise in the image
+	    // will drive insignificant, but correlated, terms in the solution.  To avoid these
+	    // problems, we can use SVD to identify numerically unconstrained values and to
+	    // avoid statistically badly determined value.
+
+	    // A = sumMatrix, B = sumVector
+	    // SVD: A = U w V^T  -> A^{-1} = V (1/w) U^T
+	    // x = V (1/w) (U^T B)
+	    // \sigma_x = sqrt(diag(A^{-1}))
+	    // solve for x and A^{-1} to get x & dx
+	    // identify the elements of (1/w) that are nan (1/0.0) -> set to 0.0
+	    // identify the elements of x that are insignificant (x / dx < 1.0? < 0.5?) -> set to 0.0
+
+	    // If I use the SVD trick to re-condition the matrix, I need to break out the
+	    // kernel and normalization terms from the background term.
+	    // XXX is this true?  or was this due to an error in the analysis?
+
+	    int bgIndex = PM_SUBTRACTION_INDEX_BG(kernels); // Index in matrix for background
+
+	    // now pull out the kernel elements into their own square matrix
+	    psImage  *kernelMatrix = psImageAlloc  (sumMatrix->numCols - 1, sumMatrix->numRows - 1, PS_TYPE_F64);
+	    psVector *kernelVector = psVectorAlloc (sumMatrix->numCols - 1, PS_TYPE_F64);
+
+	    for (int ix = 0, kx = 0; ix < sumMatrix->numCols; ix++) {
+		if (ix == bgIndex) continue;
+		for (int iy = 0, ky = 0; iy < sumMatrix->numRows; iy++) {
+		    if (iy == bgIndex) continue;
+		    kernelMatrix->data.F64[ky][kx] = sumMatrix->data.F64[iy][ix];
+		    ky++;
+		}
+		kernelVector->data.F64[kx] = sumVector->data.F64[ix];
+		kx++;
+	    }
+
+	    psImage *U = NULL;
+	    psImage *V = NULL;
+	    psVector *w = NULL;
+	    if (!psMatrixSVD (&U, &w, &V, kernelMatrix)) {
+		psError(PS_ERR_UNKNOWN, false, "failed to perform SVD on sumMatrix\n");
+		return NULL;
+	    }
+
+	    // calculate A_inverse:
+	    // Ainv = V * w * U^T
+	    psImage *wUt  = p_pmSubSolve_wUt (w, U);
+	    psImage *Ainv = p_pmSubSolve_VwUt (V, wUt);
+	    psImage *Xvar = NULL;
+	    psFree (wUt);
 
 # ifdef TESTING
-            // kernel terms:
-            for (int i = 0; i < w->n; i++) {
-                fprintf (stderr, "w: %f\n", w->data.F64[i]);
-            }
+	    // kernel terms:
+	    for (int i = 0; i < w->n; i++) {
+		fprintf (stderr, "w: %f\n", w->data.F64[i]);
+	    }
 # endif
-            // loop over w adding in more and more of the values until chisquare is no longer
-            // dropping significantly.
-            // XXX this does not seem to work very well: we seem to need all terms even for
-            // simple cases...
-
-            psVector *Xsvd = NULL;
-            {
-                psVector *Ax = NULL;
-                psVector *UtB = NULL;
-                psVector *wUtB = NULL;
-
-                psVector *wApply = psVectorAlloc(w->n, PS_TYPE_F64);
-                psVector *wMask = psVectorAlloc(w->n, PS_TYPE_U8);
-                psVectorInit (wMask, 1); // start by masking everything
-
-                double chiSquareLast = NAN;
-                int maxWeight = 0;
-
-                double Axx, Bx, y2;
-
-                // XXX this is an attempt to exclude insignificant modes.
-                // it was not successful with the ISIS kernel set: removing even
-                // the least significant mode leaves additional ringing / noise
-                // because the terms are so coupled.
-                for (int k = 0; false && (k < w->n); k++) {
-
-                    // unmask the k-th weight
-                    wMask->data.U8[k] = 0;
-                    p_pmSubSolve_SetWeights(wApply, w, wMask);
-
-                    // solve for x:
-                    // x = V * w * (U^T * B)
-                    p_pmSubSolve_UtB (&UtB, U, kernelVector);
-                    p_pmSubSolve_wUtB (&wUtB, wApply, UtB);
-                    p_pmSubSolve_VwUtB (&Xsvd, V, wUtB);
-
-                    // chi-square for this system of equations:
-                    // chi-square = sum over terms of: (Ax - B)*x - b*x - y^2
-                    // y^2 = \sum_stamps \sum_pixels input->kernel[y][x]^2
-                    p_pmSubSolve_Ax (&Ax, kernelMatrix, Xsvd);
-                    p_pmSubSolve_VdV (&Axx, Ax, Xsvd);
-                    p_pmSubSolve_VdV (&Bx, kernelVector, Xsvd);
-                    p_pmSubSolve_y2 (&y2, kernels, stamps);
-
-                    // apparently, this works (compare with the brute force value below
-                    double chiSquare = Axx - 2.0*Bx + y2;
-                    double deltaChi = (k == 0) ? chiSquare : chiSquareLast - chiSquare;
-                    chiSquareLast = chiSquare;
-
-                    // fprintf (stderr, "chi square = %f, delta: %f\n", chiSquare, deltaChi);
-                    if (k && !maxWeight && (deltaChi < 1.0)) {
-                        maxWeight = k;
-                    }
-                }
-
-                // keep all terms or we get extra ringing
-                maxWeight = w->n;
-                psVectorInit (wMask, 1);
-                for (int k = 0; k < maxWeight; k++) {
-                    wMask->data.U8[k] = 0;
-                }
-                p_pmSubSolve_SetWeights(wApply, w, wMask);
-
-                // solve for x:
-                // x = V * w * (U^T * B)
-                p_pmSubSolve_UtB (&UtB, U, kernelVector);
-                p_pmSubSolve_wUtB (&wUtB, wApply, UtB);
-                p_pmSubSolve_VwUtB (&Xsvd, V, wUtB);
-
-                // chi-square for this system of equations:
-                // chi-square = sum over terms of: (Ax - B)*x - b*x - y^2
-                // y^2 = \sum_stamps \sum_pixels input->kernel[y][x]^2
-                p_pmSubSolve_Ax (&Ax, kernelMatrix, Xsvd);
-                p_pmSubSolve_VdV (&Axx, Ax, Xsvd);
-                p_pmSubSolve_VdV (&Bx, kernelVector, Xsvd);
-                p_pmSubSolve_y2 (&y2, kernels, stamps);
-
-                // apparently, this works (compare with the brute force value below
-                double chiSquare = Axx - 2.0*Bx + y2;
-                psLogMsg ("psModules.imcombine", PS_LOG_INFO, "model kernel with %d terms; chi square = %f\n", maxWeight, chiSquare);
-
-                // re-calculate A^{-1} to get new variances:
-                // Ainv = V * w * U^T
-                // XXX since we keep all terms, this is identical to Ainv
-                psImage *wUt  = p_pmSubSolve_wUt (wApply, U);
-                Xvar = p_pmSubSolve_VwUt (V, wUt);
-                psFree (wUt);
-
-                psFree (Ax);
-                psFree (UtB);
-                psFree (wUtB);
-                psFree (wApply);
-                psFree (wMask);
-            }
-
-            // copy the kernel solutions to the full solution vector:
-            solution = psVectorAlloc(sumVector->n, PS_TYPE_F64);
-            solutionErr = psVectorAlloc(sumVector->n, PS_TYPE_F64);
-
-            for (int ix = 0, kx = 0; ix < sumVector->n; ix++) {
-                if (ix == bgIndex) {
-                    solution->data.F64[ix] = 0;
-                    solutionErr->data.F64[ix] = 0.001;
-                    continue;
-                }
-                solutionErr->data.F64[ix] = sqrt(Ainv->data.F64[kx][kx]);
-                solution->data.F64[ix] = Xsvd->data.F64[kx];
-                kx++;
-            }
-
-            psFree (kernelMatrix);
-            psFree (kernelVector);
-
-            psFree (U);
-            psFree (V);
-            psFree (w);
-
-            psFree (Ainv);
-            psFree (Xsvd);
-        } else {
-            psVector *permutation = NULL;       // Permutation vector, required for LU decomposition
-            psImage *luMatrix = psMatrixLUDecomposition(NULL, &permutation, sumMatrix);
-            if (!luMatrix) {
-                psError(PS_ERR_UNKNOWN, true, "LU Decomposition of least-squares matrix failed.\n");
-                psFree(solution);
-                psFree(sumVector);
-                psFree(sumMatrix);
-                psFree(luMatrix);
-                psFree(permutation);
-                return NULL;
-            }
-
-            solution = psMatrixLUSolution(NULL, luMatrix, sumVector, permutation);
-            psFree(luMatrix);
-            psFree(permutation);
-            if (!solution) {
-                psError(PS_ERR_UNKNOWN, true, "Failed to solve the least-squares system.\n");
-                psFree(solution);
-                psFree(sumVector);
-                psFree(sumMatrix);
-                return NULL;
-            }
-
-            // XXX LUD does not provide A^{-1}?  fake the error for now
-            solutionErr = psVectorAlloc(sumVector->n, PS_TYPE_F64);
-            for (int ix = 0; ix < sumVector->n; ix++) {
-                solutionErr->data.F64[ix] = 0.1*solution->data.F64[ix];
-            }
-        }
-
-        if (!kernels->solution1) {
-            kernels->solution1 = psVectorAlloc (sumVector->n, PS_TYPE_F64);
-            psVectorInit (kernels->solution1, 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;
-            int spatialOrder = kernels->spatialOrder;       // Order of spatial variation
-            int numPoly = PM_SUBTRACTION_POLYTERMS(spatialOrder); // Number of polynomial terms
-            for (int i = 0; i < numKernels * numPoly; i++) {
-                // XXX fprintf (stderr, "%f +/- %f (%f) -> ", solution->data.F64[i], solutionErr->data.F64[i], fabs(solution->data.F64[i]/solutionErr->data.F64[i]));
-                if (fabs(solution->data.F64[i] / solutionErr->data.F64[i]) < COEFF_SIG) {
-                    // XXX fprintf (stderr, "drop\n");
-                    kernels->solution1->data.F64[i] = 0.0;
-                } else {
-                    // XXX fprintf (stderr, "keep\n");
-                    kernels->solution1->data.F64[i] = solution->data.F64[i];
-                }
-            }
-        }
-        // double chiSquare = p_pmSubSolve_ChiSquare (kernels, stamps);
-        // fprintf (stderr, "chi square Brute = %f\n", chiSquare);
-
-        psFree(solution);
-        psFree(sumVector);
-        psFree(sumMatrix);
+	    // loop over w adding in more and more of the values until chisquare is no longer
+	    // dropping significantly.
+	    // XXX this does not seem to work very well: we seem to need all terms even for
+	    // simple cases...
+
+	    psVector *Xsvd = NULL;
+	    { 
+		psVector *Ax = NULL;
+		psVector *UtB = NULL;
+		psVector *wUtB = NULL;
+
+		psVector *wApply = psVectorAlloc(w->n, PS_TYPE_F64);
+		psVector *wMask = psVectorAlloc(w->n, PS_TYPE_U8);
+		psVectorInit (wMask, 1); // start by masking everything
+
+		double chiSquareLast = NAN;
+		int maxWeight = 0;
+
+		double Axx, Bx, y2;
+
+		// XXX this is an attempt to exclude insignificant modes.
+		// it was not successful with the ISIS kernel set: removing even 
+		// the least significant mode leaves additional ringing / noise
+		// because the terms are so coupled.
+		for (int k = 0; false && (k < w->n); k++) {
+
+		    // unmask the k-th weight
+		    wMask->data.U8[k] = 0;
+		    p_pmSubSolve_SetWeights(wApply, w, wMask);
+
+		    // solve for x:
+		    // x = V * w * (U^T * B)
+		    p_pmSubSolve_UtB (&UtB, U, kernelVector);
+		    p_pmSubSolve_wUtB (&wUtB, wApply, UtB);
+		    p_pmSubSolve_VwUtB (&Xsvd, V, wUtB);
+
+		    // chi-square for this system of equations:
+		    // chi-square = sum over terms of: (Ax - B)*x - b*x - y^2
+		    // y^2 = \sum_stamps \sum_pixels input->kernel[y][x]^2
+		    p_pmSubSolve_Ax (&Ax, kernelMatrix, Xsvd);
+		    p_pmSubSolve_VdV (&Axx, Ax, Xsvd);
+		    p_pmSubSolve_VdV (&Bx, kernelVector, Xsvd);
+		    p_pmSubSolve_y2 (&y2, kernels, stamps);
+
+		    // apparently, this works (compare with the brute force value below
+		    double chiSquare = Axx - 2.0*Bx + y2;
+		    double deltaChi = (k == 0) ? chiSquare : chiSquareLast - chiSquare;
+		    chiSquareLast = chiSquare;
+
+		    // fprintf (stderr, "chi square = %f, delta: %f\n", chiSquare, deltaChi);
+		    if (k && !maxWeight && (deltaChi < 1.0)) {
+			maxWeight = k;
+		    }
+		}
+
+		// keep all terms or we get extra ringing
+		maxWeight = w->n;
+		psVectorInit (wMask, 1);
+		for (int k = 0; k < maxWeight; k++) {
+		    wMask->data.U8[k] = 0;
+		}
+		p_pmSubSolve_SetWeights(wApply, w, wMask);
+
+		// solve for x:
+		// x = V * w * (U^T * B)
+		p_pmSubSolve_UtB (&UtB, U, kernelVector);
+		p_pmSubSolve_wUtB (&wUtB, wApply, UtB);
+		p_pmSubSolve_VwUtB (&Xsvd, V, wUtB);
+
+		// chi-square for this system of equations:
+		// chi-square = sum over terms of: (Ax - B)*x - b*x - y^2
+		// y^2 = \sum_stamps \sum_pixels input->kernel[y][x]^2
+		p_pmSubSolve_Ax (&Ax, kernelMatrix, Xsvd);
+		p_pmSubSolve_VdV (&Axx, Ax, Xsvd);
+		p_pmSubSolve_VdV (&Bx, kernelVector, Xsvd);
+		p_pmSubSolve_y2 (&y2, kernels, stamps);
+
+		// apparently, this works (compare with the brute force value below
+		double chiSquare = Axx - 2.0*Bx + y2;
+		psLogMsg ("psModules.imcombine", PS_LOG_INFO, "model kernel with %d terms; chi square = %f\n", maxWeight, chiSquare);
+
+		// re-calculate A^{-1} to get new variances:
+		// Ainv = V * w * U^T
+		// XXX since we keep all terms, this is identical to Ainv
+		psImage *wUt  = p_pmSubSolve_wUt (wApply, U);
+		Xvar = p_pmSubSolve_VwUt (V, wUt);
+		psFree (wUt);
+
+		psFree (Ax);
+		psFree (UtB);
+		psFree (wUtB);
+		psFree (wApply);
+		psFree (wMask);
+	    }
+
+	    // copy the kernel solutions to the full solution vector:
+	    solution = psVectorAlloc(sumVector->n, PS_TYPE_F64);
+	    solutionErr = psVectorAlloc(sumVector->n, PS_TYPE_F64);
+
+	    for (int ix = 0, kx = 0; ix < sumVector->n; ix++) {
+		if (ix == bgIndex) {
+		    solution->data.F64[ix] = 0;
+		    solutionErr->data.F64[ix] = 0.001;
+		    continue;
+		}
+		solutionErr->data.F64[ix] = sqrt(Ainv->data.F64[kx][kx]);
+		solution->data.F64[ix] = Xsvd->data.F64[kx];
+		kx++;
+	    }
+
+	    psFree (kernelMatrix);
+	    psFree (kernelVector);
+
+	    psFree (U);
+	    psFree (V);
+	    psFree (w);
+
+	    psFree (Ainv);
+	    psFree (Xsvd);
+	} else {
+	    psVector *permutation = NULL;       // Permutation vector, required for LU decomposition
+	    psImage *luMatrix = psMatrixLUDecomposition(NULL, &permutation, sumMatrix);
+	    if (!luMatrix) {
+		psError(PS_ERR_UNKNOWN, true, "LU Decomposition of least-squares matrix failed.\n");
+		psFree(solution);
+		psFree(sumVector);
+		psFree(sumMatrix);
+		psFree(luMatrix);
+		psFree(permutation);
+		return NULL;
+	    }
+
+	    solution = psMatrixLUSolution(NULL, luMatrix, sumVector, permutation);
+	    psFree(luMatrix);
+	    psFree(permutation);
+	    if (!solution) {
+		psError(PS_ERR_UNKNOWN, true, "Failed to solve the least-squares system.\n");
+		psFree(solution);
+		psFree(sumVector);
+		psFree(sumMatrix);
+		return NULL;
+	    }
+
+	    // XXX LUD does not provide A^{-1}?  fake the error for now 
+	    solutionErr = psVectorAlloc(sumVector->n, PS_TYPE_F64);
+	    for (int ix = 0; ix < sumVector->n; ix++) {
+		solutionErr->data.F64[ix] = 0.1*solution->data.F64[ix];
+	    }
+	}
+
+	if (!kernels->solution1) {
+	    kernels->solution1 = psVectorAlloc (sumVector->n, PS_TYPE_F64);
+	    psVectorInit (kernels->solution1, 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;
+	    int spatialOrder = kernels->spatialOrder;       // Order of spatial variation
+	    int numPoly = PM_SUBTRACTION_POLYTERMS(spatialOrder); // Number of polynomial terms
+	    for (int i = 0; i < numKernels * numPoly; i++) {
+		// XXX fprintf (stderr, "%f +/- %f (%f) -> ", solution->data.F64[i], solutionErr->data.F64[i], fabs(solution->data.F64[i]/solutionErr->data.F64[i]));
+		if (fabs(solution->data.F64[i] / solutionErr->data.F64[i]) < COEFF_SIG) {
+		    // XXX fprintf (stderr, "drop\n");
+		    kernels->solution1->data.F64[i] = 0.0;
+		} else {
+		    // XXX fprintf (stderr, "keep\n");
+		    kernels->solution1->data.F64[i] = solution->data.F64[i];
+		}
+	    }
+	}
+	// double chiSquare = p_pmSubSolve_ChiSquare (kernels, stamps);
+	// fprintf (stderr, "chi square Brute = %f\n", chiSquare);
+
+	psFree(solution);
+	psFree(sumVector);
+	psFree(sumMatrix);
 
 #ifdef TESTING
@@ -1201,5 +1200,4 @@
             }
         }
-
 
 #ifdef TESTING
@@ -1319,8 +1317,8 @@
             psVectorWriteFile("sumVectorFix.dat", sumVector);
 #endif
-
-        }
-#endif
-
+        }
+#endif
+    }
+#endif
         solution = psMatrixSolveSVD(solution, sumMatrix, sumVector, NAN);
 
@@ -1369,5 +1367,5 @@
 
 psVector *pmSubtractionCalculateDeviations(pmSubtractionStampList *stamps,
-                                           const pmSubtractionKernels *kernels)
+                                           pmSubtractionKernels *kernels)
 {
     PM_ASSERT_SUBTRACTION_STAMP_LIST_NON_NULL(stamps, NULL);
@@ -1387,177 +1385,244 @@
     pmSubtractionVisualShowFitInit (stamps);
 
+    psVector *fSigRes = psVectorAllocEmpty(stamps->num, PS_TYPE_F32); 
+    psVector *fMinRes = psVectorAllocEmpty(stamps->num, PS_TYPE_F32); 
+    psVector *fMaxRes = psVectorAllocEmpty(stamps->num, PS_TYPE_F32); 
+
     for (int i = 0; i < stamps->num; i++) {
-        pmSubtractionStamp *stamp = stamps->stamps->data[i]; // The stamp of interest
-        if (stamp->status != PM_SUBTRACTION_STAMP_USED) {
-            deviations->data.F32[i] = NAN;
-            continue;
-        }
-
-        // Calculate coefficients of the kernel basis functions
-        polyValues = p_pmSubtractionPolynomial(polyValues, kernels->spatialOrder, stamp->xNorm, stamp->yNorm);
-        double norm = p_pmSubtractionSolutionNorm(kernels); // Normalisation
-        double background = p_pmSubtractionSolutionBackground(kernels, polyValues);// Difference in background
-
-        // Calculate residuals
-        psKernel *weight = stamp->weight; // Weight postage stamp
-        psImageInit(residual->image, 0.0);
-        if (kernels->mode != PM_SUBTRACTION_MODE_DUAL) {
-            psKernel *target;           // Target postage stamp
-            psKernel *source;           // Source postage stamp
-            psArray *convolutions;      // Convolution postage stamps for each kernel basis function
-            switch (kernels->mode) {
-              case PM_SUBTRACTION_MODE_1:
-                target = stamp->image2;
-                source = stamp->image1;
-                convolutions = stamp->convolutions1;
-
-                // Having convolved image1 and changed its normalisation, we need to renormalise the residual
-                // so that it is on the scale of image1.
-                psImage *image = pmSubtractionKernelImage(kernels, stamp->xNorm, stamp->yNorm,
-                                                          false); // Kernel image
-                if (!image) {
-                    psError(PS_ERR_UNKNOWN, false, "Unable to generate image of kernel.");
-                    return false;
-                }
-                double sumKernel = 0;   // Sum of kernel, for normalising residual
-                int size = kernels->size; // Half-size of kernel
-                int fullSize = 2 * size + 1; // Full size of kernel
-                for (int y = 0; y < fullSize; y++) {
-                    for (int x = 0; x < fullSize; x++) {
-                        sumKernel += image->data.F32[y][x];
-                    }
-                }
-                psFree(image);
-                devNorm = 1.0 / sumKernel / numPixels;
-                break;
-              case PM_SUBTRACTION_MODE_2:
-                target = stamp->image1;
-                source = stamp->image2;
-                convolutions = stamp->convolutions2;
-                break;
-              default:
-                psAbort("Unsupported subtraction mode: %x", kernels->mode);
-            }
-
-            for (int j = 0; j < numKernels; j++) {
-                psKernel *convolution = convolutions->data[j]; // Convolution
-                double coefficient = p_pmSubtractionSolutionCoeff(kernels, polyValues, j,
-                                                                  false); // Coefficient
-                for (int y = - footprint; y <= footprint; y++) {
-                    for (int x = - footprint; x <= footprint; x++) {
-                        residual->kernel[y][x] += convolution->kernel[y][x] * coefficient;
-                    }
-                }
-            }
-
-            // XXX visualize the target, source, convolution and residual
-            pmSubtractionVisualShowFitAddStamp (target, source, residual, background, norm, i);
-
-            for (int y = - footprint; y <= footprint; y++) {
-                for (int x = - footprint; x <= footprint; x++) {
-                    residual->kernel[y][x] += background + source->kernel[y][x] * norm - target->kernel[y][x];
-                }
-            }
-        } else {
-            // Dual convolution
-            psArray *convolutions1 = stamp->convolutions1; // Convolutions of the first image
-            psArray *convolutions2 = stamp->convolutions2; // Convolutions of the second image
-            psKernel *image1 = stamp->image1; // The first image
-            psKernel *image2 = stamp->image2; // The second image
-
-            for (int j = 0; j < numKernels; j++) {
-                psKernel *conv1 = convolutions1->data[j]; // Convolution of first image
-                psKernel *conv2 = convolutions2->data[j]; // Convolution of second image
-                double coeff1 = p_pmSubtractionSolutionCoeff(kernels, polyValues, j, false); // Coefficient 1
-                double coeff2 = p_pmSubtractionSolutionCoeff(kernels, polyValues, j, true); // Coefficient 2
-
-                for (int y = - footprint; y <= footprint; y++) {
-                    for (int x = - footprint; x <= footprint; x++) {
-                        residual->kernel[y][x] += conv2->kernel[y][x] * coeff2 + conv1->kernel[y][x] * coeff1;
-                    }
-                }
-            }
-
-            // XXX visualize the target, source, convolution and residual
-            pmSubtractionVisualShowFitAddStamp (image2, image1, residual, background, norm, i);
-
-            for (int y = - footprint; y <= footprint; y++) {
-                for (int x = - footprint; x <= footprint; x++) {
-                    residual->kernel[y][x] += background + image1->kernel[y][x] * norm - image2->kernel[y][x];
-                }
-            }
-        }
-
-        double deviation = 0.0;         // Sum of differences
-        for (int y = - footprint; y <= footprint; y++) {
-            for (int x = - footprint; x <= footprint; x++) {
-                double dev = PS_SQR(residual->kernel[y][x]) * weight->kernel[y][x];
-                deviation += dev;
+	pmSubtractionStamp *stamp = stamps->stamps->data[i]; // The stamp of interest
+	if (stamp->status != PM_SUBTRACTION_STAMP_USED) {
+	    deviations->data.F32[i] = NAN;
+	    continue;
+	}
+
+	// Calculate coefficients of the kernel basis functions
+	polyValues = p_pmSubtractionPolynomial(polyValues, kernels->spatialOrder, stamp->xNorm, stamp->yNorm);
+	double norm = p_pmSubtractionSolutionNorm(kernels); // Normalisation
+	double background = p_pmSubtractionSolutionBackground(kernels, polyValues);// Difference in background
+
+	// Calculate residuals
+	psKernel *weight = stamp->weight; // Weight postage stamp
+	psImageInit(residual->image, 0.0);
+	if (kernels->mode != PM_SUBTRACTION_MODE_DUAL) {
+	    psKernel *target;           // Target postage stamp
+	    psKernel *source;           // Source postage stamp
+	    psArray *convolutions;      // Convolution postage stamps for each kernel basis function
+	    switch (kernels->mode) {
+	      case PM_SUBTRACTION_MODE_1:
+		target = stamp->image2;
+		source = stamp->image1;
+		convolutions = stamp->convolutions1;
+
+		// Having convolved image1 and changed its normalisation, we need to renormalise the residual
+		// so that it is on the scale of image1.
+		psImage *image = pmSubtractionKernelImage(kernels, stamp->xNorm, stamp->yNorm,
+							  false); // Kernel image
+		if (!image) {
+		    psError(PS_ERR_UNKNOWN, false, "Unable to generate image of kernel.");
+		    return false;
+		}
+		double sumKernel = 0;   // Sum of kernel, for normalising residual
+		int size = kernels->size; // Half-size of kernel
+		int fullSize = 2 * size + 1; // Full size of kernel
+		for (int y = 0; y < fullSize; y++) {
+		    for (int x = 0; x < fullSize; x++) {
+			sumKernel += image->data.F32[y][x];
+		    }
+		}
+		psFree(image);
+		devNorm = 1.0 / sumKernel / numPixels;
+		break;
+	      case PM_SUBTRACTION_MODE_2:
+		target = stamp->image1;
+		source = stamp->image2;
+		convolutions = stamp->convolutions2;
+		break;
+	      default:
+		psAbort("Unsupported subtraction mode: %x", kernels->mode);
+	    }
+
+	    for (int j = 0; j < numKernels; j++) {
+		psKernel *convolution = convolutions->data[j]; // Convolution
+		double coefficient = p_pmSubtractionSolutionCoeff(kernels, polyValues, j,
+								  false); // Coefficient
+		for (int y = - footprint; y <= footprint; y++) {
+		    for (int x = - footprint; x <= footprint; x++) {
+			residual->kernel[y][x] += convolution->kernel[y][x] * coefficient;
+		    }
+		}
+	    }
+
+	    // XXX visualize the target, source, convolution and residual
+	    pmSubtractionVisualShowFitAddStamp (target, source, residual, background, norm, i);
+
+	    for (int y = - footprint; y <= footprint; y++) {
+		for (int x = - footprint; x <= footprint; x++) {
+		    residual->kernel[y][x] += background + source->kernel[y][x] * norm - target->kernel[y][x];
+		}
+	    }
+
+	    // XXX measure some useful stats on the residuals
+	    float sum = 0.0;
+	    float peak = 0.0;
+	    for (int y = - footprint; y <= footprint; y++) {
+		for (int x = - footprint; x <= footprint; x++) {
+		    sum += 0.5*(target->kernel[y][x] + source->kernel[y][x] * norm);
+		    peak = PS_MAX(peak, 0.5*(target->kernel[y][x] + source->kernel[y][x] * norm));
+		}
+	    }
+
+	    // only count pixels with more than X% of the source flux
+	    // calculate stdev(dflux)
+	    float dflux1 = 0.0;
+	    float dflux2 = 0.0;
+	    int npix = 0;
+
+	    float dmax = 0.0;
+	    float dmin = 0.0;
+
+	    for (int y = - footprint; y <= footprint; y++) {
+		for (int x = - footprint; x <= footprint; x++) {
+		    float dflux = 0.5*(target->kernel[y][x] + source->kernel[y][x] * norm);
+		    if (dflux < 0.02*sum) continue;
+		    dflux1 += residual->kernel[y][x];
+		    dflux2 += PS_SQR(residual->kernel[y][x]);
+		    dmax = PS_MAX(residual->kernel[y][x], dmax);
+		    dmin = PS_MIN(residual->kernel[y][x], dmin);
+		    npix ++;
+		}
+	    }
+	    float sigma = sqrt(dflux2 / npix - PS_SQR(dflux1/npix));
+	    // fprintf (stderr, "sum: %f, peak: %f, sigma: %f, fsigma: %f, fmax: %f, fmin: %f\n", sum, peak, sigma, sigma/sum, dmax/peak, dmin/peak);
+	    psVectorAppend(fSigRes, sigma/sum);
+	    psVectorAppend(fMaxRes, dmax/peak);
+	    psVectorAppend(fMinRes, dmin/peak);
+	} else {
+	    // Dual convolution
+	    psArray *convolutions1 = stamp->convolutions1; // Convolutions of the first image
+	    psArray *convolutions2 = stamp->convolutions2; // Convolutions of the second image
+	    psKernel *image1 = stamp->image1; // The first image
+	    psKernel *image2 = stamp->image2; // The second image
+
+	    for (int j = 0; j < numKernels; j++) {
+		psKernel *conv1 = convolutions1->data[j]; // Convolution of first image
+		psKernel *conv2 = convolutions2->data[j]; // Convolution of second image
+		double coeff1 = p_pmSubtractionSolutionCoeff(kernels, polyValues, j, false); // Coefficient 1
+		double coeff2 = p_pmSubtractionSolutionCoeff(kernels, polyValues, j, true); // Coefficient 2
+
+		for (int y = - footprint; y <= footprint; y++) {
+		    for (int x = - footprint; x <= footprint; x++) {
+			residual->kernel[y][x] += conv2->kernel[y][x] * coeff2 + conv1->kernel[y][x] * coeff1;
+		    }
+		}
+	    }
+
+	    // XXX visualize the target, source, convolution and residual
+	    pmSubtractionVisualShowFitAddStamp (image2, image1, residual, background, norm, i);
+
+	    for (int y = - footprint; y <= footprint; y++) {
+		for (int x = - footprint; x <= footprint; x++) {
+		    residual->kernel[y][x] += background + image1->kernel[y][x] * norm - image2->kernel[y][x];
+		}
+	    }
+	}
+
+	double deviation = 0.0;         // Sum of differences
+	for (int y = - footprint; y <= footprint; y++) {
+	    for (int x = - footprint; x <= footprint; x++) {
+		double dev = PS_SQR(residual->kernel[y][x]) * weight->kernel[y][x];
+		deviation += dev;
 #ifdef TESTING
-                residual->kernel[y][x] = dev;
-#endif
-            }
-        }
-        deviations->data.F32[i] = devNorm * deviation;
-        psTrace("psModules.imcombine", 5, "Deviation for stamp %d (%d,%d): %f\n",
-                i, (int)(stamp->x + 0.5), (int)(stamp->y + 0.5), deviations->data.F32[i]);
-        if (!isfinite(deviations->data.F32[i])) {
-            stamp->status = PM_SUBTRACTION_STAMP_REJECTED;
-            psTrace("psModules.imcombine", 5,
-                    "Rejecting stamp %d (%d,%d) because of non-finite deviation\n",
-                    i, (int)(stamp->x + 0.5), (int)(stamp->y + 0.5));
-            continue;
-        }
+		residual->kernel[y][x] = dev;
+#endif
+	    }
+	}
+	deviations->data.F32[i] = devNorm * deviation;
+	psTrace("psModules.imcombine", 5, "Deviation for stamp %d (%d,%d): %f\n",
+		i, (int)(stamp->x + 0.5), (int)(stamp->y + 0.5), deviations->data.F32[i]);
+	if (!isfinite(deviations->data.F32[i])) {
+	    stamp->status = PM_SUBTRACTION_STAMP_REJECTED;
+	    psTrace("psModules.imcombine", 5,
+		    "Rejecting stamp %d (%d,%d) because of non-finite deviation\n",
+		    i, (int)(stamp->x + 0.5), (int)(stamp->y + 0.5));
+	    continue;
+	}
 
 #ifdef TESTING
-        {
-            psString filename = NULL;
-            psStringAppend(&filename, "resid_%03d.fits", i);
-            psFits *fits = psFitsOpen(filename, "w");
-            psFree(filename);
-            psFitsWriteImage(fits, NULL, residual->image, 0, NULL);
-            psFitsClose(fits);
-        }
-        if (stamp->image1) {
-            psString filename = NULL;
-            psStringAppend(&filename, "stamp_image1_%03d.fits", i);
-            psFits *fits = psFitsOpen(filename, "w");
-            psFree(filename);
-            psFitsWriteImage(fits, NULL, stamp->image1->image, 0, NULL);
-            psFitsClose(fits);
-        }
-        if (stamp->image2) {
-            psString filename = NULL;
-            psStringAppend(&filename, "stamp_image2_%03d.fits", i);
-            psFits *fits = psFitsOpen(filename, "w");
-            psFree(filename);
-            psFitsWriteImage(fits, NULL, stamp->image2->image, 0, NULL);
-            psFitsClose(fits);
-        }
-        if (stamp->weight) {
-            psString filename = NULL;
-            psStringAppend(&filename, "stamp_weight_%03d.fits", i);
-            psFits *fits = psFitsOpen(filename, "w");
-            psFree(filename);
-            psFitsWriteImage(fits, NULL, stamp->weight->image, 0, NULL);
-            psFitsClose(fits);
-        }
-#endif
-
+	{
+	    psString filename = NULL;
+	    psStringAppend(&filename, "resid_%03d.fits", i);
+	    psFits *fits = psFitsOpen(filename, "w");
+	    psFree(filename);
+	    psFitsWriteImage(fits, NULL, residual->image, 0, NULL);
+	    psFitsClose(fits);
+	}
+	if (stamp->image1) {
+	    psString filename = NULL;
+	    psStringAppend(&filename, "stamp_image1_%03d.fits", i);
+	    psFits *fits = psFitsOpen(filename, "w");
+	    psFree(filename);
+	    psFitsWriteImage(fits, NULL, stamp->image1->image, 0, NULL);
+	    psFitsClose(fits);
+	}
+	if (stamp->image2) {
+	    psString filename = NULL;
+	    psStringAppend(&filename, "stamp_image2_%03d.fits", i);
+	    psFits *fits = psFitsOpen(filename, "w");
+	    psFree(filename);
+	    psFitsWriteImage(fits, NULL, stamp->image2->image, 0, NULL);
+	    psFitsClose(fits);
+	}
+	if (stamp->weight) {
+	    psString filename = NULL;
+	    psStringAppend(&filename, "stamp_weight_%03d.fits", i);
+	    psFits *fits = psFitsOpen(filename, "w");
+	    psFree(filename);
+	    psFitsWriteImage(fits, NULL, stamp->weight->image, 0, NULL);
+	    psFitsClose(fits);
+	}
+#endif
+    
     }
 
     // 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(norm);
-        pmSubtractionVisualPlotFit(kernels);
+    { 
+	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(norm);
+	pmSubtractionVisualPlotFit(kernels);
+
+	psStats *stats = psStatsAlloc(PS_STAT_ROBUST_MEDIAN | PS_STAT_ROBUST_STDEV);
+	psVectorStats (stats, fSigRes, NULL, NULL, 0);
+	kernels->fSigResMean = stats->robustMedian;
+	kernels->fSigResStdev = stats->robustStdev;
+
+	psStatsInit (stats);
+	psVectorStats (stats, fMaxRes, NULL, NULL, 0);
+	kernels->fMaxResMean = stats->robustMedian;
+	kernels->fMaxResStdev = stats->robustStdev;
+
+	psStatsInit (stats);
+	psVectorStats (stats, fMinRes, NULL, NULL, 0);
+	kernels->fMinResMean = stats->robustMedian;
+	kernels->fMinResStdev = stats->robustStdev;
+
+	// XXX save these values somewhere
+	psLogMsg("psModules.imcombine", PS_LOG_INFO, "fSigma: %f +/- %f, fMaxRes: %f +/- %f, fMinRes: %f +/- %f", 
+		 kernels->fSigResMean, kernels->fSigResStdev, 
+		 kernels->fMaxResMean, kernels->fMaxResStdev, 
+		 kernels->fMinResMean, kernels->fMinResStdev);
+
+	psFree (fSigRes);
+	psFree (fMaxRes);
+	psFree (fMinRes);
+	psFree (stats);
     }
 
     psFree(residual);
     psFree(polyValues);
+
 
     return deviations;
@@ -1574,9 +1639,9 @@
 
     for (int i = 0; i < wUt->numCols; i++) {
-        for (int j = 0; j < wUt->numRows; j++) {
-            if (!isfinite(w->data.F64[j])) continue;
-            if (w->data.F64[j] == 0.0) continue;
-            wUt->data.F64[j][i] = U->data.F64[i][j] / w->data.F64[j];
-        }
+	for (int j = 0; j < wUt->numRows; j++) {
+	    if (!isfinite(w->data.F64[j])) continue;
+	    if (w->data.F64[j] == 0.0) continue;
+	    wUt->data.F64[j][i] = U->data.F64[i][j] / w->data.F64[j];
+	}
     }
     return wUt;
@@ -1591,11 +1656,11 @@
 
     for (int i = 0; i < Ainv->numCols; i++) {
-        for (int j = 0; j < Ainv->numRows; j++) {
-            double sum = 0.0;
-            for (int k = 0; k < V->numCols; k++) {
-                sum += V->data.F64[j][k] * wUt->data.F64[k][i];
-            }
-            Ainv->data.F64[j][i] = sum;
-        }
+	for (int j = 0; j < Ainv->numRows; j++) {
+	    double sum = 0.0;
+	    for (int k = 0; k < V->numCols; k++) {
+		sum += V->data.F64[j][k] * wUt->data.F64[k][i];
+	    }
+	    Ainv->data.F64[j][i] = sum;
+	}
     }
     return Ainv;
@@ -1610,9 +1675,9 @@
 
     for (int i = 0; i < U->numCols; i++) {
-        double sum = 0.0;
-        for (int j = 0; j < U->numRows; j++) {
-            sum += B->data.F64[j] * U->data.F64[j][i];
-        }
-        UtB[0]->data.F64[i] = sum;
+	double sum = 0.0;
+	for (int j = 0; j < U->numRows; j++) {
+	    sum += B->data.F64[j] * U->data.F64[j][i];
+	}
+	UtB[0]->data.F64[i] = sum;
     }
     return true;
@@ -1629,7 +1694,7 @@
 
     for (int i = 0; i < w->n; i++) {
-        if (!isfinite(w->data.F64[i])) continue;
-        if (w->data.F64[i] == 0.0) continue;
-        wUtB[0]->data.F64[i] = UtB->data.F64[i] / w->data.F64[i];
+	if (!isfinite(w->data.F64[i])) continue;
+	if (w->data.F64[i] == 0.0) continue;
+	wUtB[0]->data.F64[i] = UtB->data.F64[i] / w->data.F64[i];
     }
     return true;
@@ -1644,9 +1709,9 @@
 
     for (int j = 0; j < V->numRows; j++) {
-        double sum = 0.0;
-        for (int i = 0; i < V->numCols; i++) {
-            sum += V->data.F64[j][i] * wUtB->data.F64[i];
-        }
-        VwUtB[0]->data.F64[j] = sum;
+	double sum = 0.0;
+	for (int i = 0; i < V->numCols; i++) {
+	    sum += V->data.F64[j][i] * wUtB->data.F64[i];
+	}
+	VwUtB[0]->data.F64[j] = sum;
     }
     return true;
@@ -1661,9 +1726,9 @@
 
     for (int j = 0; j < A->numRows; j++) {
-        double sum = 0.0;
-        for (int i = 0; i < A->numCols; i++) {
-            sum += A->data.F64[j][i] * x->data.F64[i];
-        }
-        B[0]->data.F64[j] = sum;
+	double sum = 0.0;
+	for (int i = 0; i < A->numCols; i++) {
+	    sum += A->data.F64[j][i] * x->data.F64[i];
+	}
+	B[0]->data.F64[j] = sum;
     }
     return true;
@@ -1677,5 +1742,5 @@
     double sum = 0.0;
     for (int i = 0; i < x->n; i++) {
-        sum += x->data.F64[i] * y->data.F64[i];
+	sum += x->data.F64[i] * y->data.F64[i];
     }
     *value = sum;
@@ -1690,45 +1755,45 @@
     for (int i = 0; i < stamps->num; i++) {
 
-        pmSubtractionStamp *stamp = stamps->stamps->data[i];
-        if (stamp->status != PM_SUBTRACTION_STAMP_USED) continue;
-
-        psKernel *weight = NULL;
-        psKernel *window = NULL;
-        psKernel *input = NULL;
+	pmSubtractionStamp *stamp = stamps->stamps->data[i];
+	if (stamp->status != PM_SUBTRACTION_STAMP_USED) continue;
+
+	psKernel *weight = NULL;
+	psKernel *window = NULL;
+	psKernel *input = NULL;
 
 #ifdef USE_WEIGHT
-        weight = stamp->weight;
+	weight = stamp->weight;
 #endif
 #ifdef USE_WINDOW
-        window = stamps->window;
-#endif
-
-        switch (kernels->mode) {
-            // MODE_1 : convolve image 1 to match image 2 (and vice versa)
-          case PM_SUBTRACTION_MODE_1:
-            input = stamp->image2;
-            break;
-          case PM_SUBTRACTION_MODE_2:
-            input = stamp->image1;
-            break;
-          default:
-            psAbort ("programming error");
-        }
-
-        for (int y = - footprint; y <= footprint; y++) {
-            for (int x = - footprint; x <= footprint; x++) {
-                double in = input->kernel[y][x];
-                double value = in*in;
-                if (weight) {
-                    float wtVal = weight->kernel[y][x];
-                    value *= wtVal;
-                }
-                if (window) {
-                    float  winVal = window->kernel[y][x];
-                    value *= winVal;
-                }
-                sum += value;
-            }
-        }
+	window = stamps->window;
+#endif
+
+	switch (kernels->mode) {
+	    // MODE_1 : convolve image 1 to match image 2 (and vice versa)
+	  case PM_SUBTRACTION_MODE_1:
+	    input = stamp->image2;
+	    break;
+	  case PM_SUBTRACTION_MODE_2:
+	    input = stamp->image1;
+	    break;
+	  default:
+	    psAbort ("programming error");
+	}
+
+	for (int y = - footprint; y <= footprint; y++) {
+	    for (int x = - footprint; x <= footprint; x++) {
+		double in = input->kernel[y][x];
+		double value = in*in;
+		if (weight) {
+		    float wtVal = weight->kernel[y][x];
+		    value *= wtVal;
+		}
+		if (window) {
+		    float  winVal = window->kernel[y][x];
+		    value *= winVal;
+		}
+		sum += value;
+	    }
+	}
     }
     *y2 = sum;
@@ -1750,68 +1815,68 @@
     for (int i = 0; i < stamps->num; i++) {
 
-        pmSubtractionStamp *stamp = stamps->stamps->data[i];
-        if (stamp->status != PM_SUBTRACTION_STAMP_USED) continue;
-
-        psKernel *weight = NULL;
-        psKernel *window = NULL;
-        psKernel *target = NULL;
-        psKernel *source = NULL;
-
-        psArray *convolutions = NULL;
+	pmSubtractionStamp *stamp = stamps->stamps->data[i];
+	if (stamp->status != PM_SUBTRACTION_STAMP_USED) continue;
+
+	psKernel *weight = NULL;
+	psKernel *window = NULL;
+	psKernel *target = NULL;
+	psKernel *source = NULL;
+
+	psArray *convolutions = NULL;
 
 #ifdef USE_WEIGHT
-        weight = stamp->weight;
+	weight = stamp->weight;
 #endif
 #ifdef USE_WINDOW
-        window = stamps->window;
-#endif
-
-        switch (kernels->mode) {
-            // MODE_1 : convolve image 1 to match image 2 (and vice versa)
-          case PM_SUBTRACTION_MODE_1:
-            target = stamp->image2;
-            source = stamp->image1;
-            convolutions = stamp->convolutions1;
-            break;
-          case PM_SUBTRACTION_MODE_2:
-            target = stamp->image1;
-            source = stamp->image2;
-            convolutions = stamp->convolutions2;
-            break;
-          default:
-            psAbort ("programming error");
-        }
-
-        // Calculate coefficients of the kernel basis functions
-        polyValues = p_pmSubtractionPolynomial(polyValues, kernels->spatialOrder, stamp->xNorm, stamp->yNorm);
-        double norm = p_pmSubtractionSolutionNorm(kernels); // Normalisation
-        double background = p_pmSubtractionSolutionBackground(kernels, polyValues);// Difference in background
-
-        psImageInit(residual->image, 0.0);
-        for (int j = 0; j < numKernels; j++) {
-            psKernel *convolution = convolutions->data[j]; // Convolution
-            double coefficient = p_pmSubtractionSolutionCoeff(kernels, polyValues, j, false); // Coefficient
-            for (int y = - footprint; y <= footprint; y++) {
-                for (int x = - footprint; x <= footprint; x++) {
-                    residual->kernel[y][x] -= convolution->kernel[y][x] * coefficient;
-                }
-            }
-        }
-
-        for (int y = - footprint; y <= footprint; y++) {
-            for (int x = - footprint; x <= footprint; x++) {
-                double resid = target->kernel[y][x] - background - source->kernel[y][x] * norm + residual->kernel[y][x];
-                double value = PS_SQR(resid);
-                if (weight) {
-                    float wtVal = weight->kernel[y][x];
-                    value *= wtVal;
-                }
-                if (window) {
-                    float  winVal = window->kernel[y][x];
-                    value *= winVal;
-                }
-                sum += value;
-            }
-        }
+	window = stamps->window;
+#endif
+
+	switch (kernels->mode) {
+	    // MODE_1 : convolve image 1 to match image 2 (and vice versa)
+	  case PM_SUBTRACTION_MODE_1:
+	    target = stamp->image2;
+	    source = stamp->image1;
+	    convolutions = stamp->convolutions1;
+	    break;
+	  case PM_SUBTRACTION_MODE_2:
+	    target = stamp->image1;
+	    source = stamp->image2;
+	    convolutions = stamp->convolutions2;
+	    break;
+	  default:
+	    psAbort ("programming error");
+	}
+
+	// Calculate coefficients of the kernel basis functions
+	polyValues = p_pmSubtractionPolynomial(polyValues, kernels->spatialOrder, stamp->xNorm, stamp->yNorm);
+	double norm = p_pmSubtractionSolutionNorm(kernels); // Normalisation
+	double background = p_pmSubtractionSolutionBackground(kernels, polyValues);// Difference in background
+
+	psImageInit(residual->image, 0.0);
+	for (int j = 0; j < numKernels; j++) {
+	    psKernel *convolution = convolutions->data[j]; // Convolution
+	    double coefficient = p_pmSubtractionSolutionCoeff(kernels, polyValues, j, false); // Coefficient
+	    for (int y = - footprint; y <= footprint; y++) {
+		for (int x = - footprint; x <= footprint; x++) {
+		    residual->kernel[y][x] -= convolution->kernel[y][x] * coefficient;
+		}
+	    }
+	}
+
+	for (int y = - footprint; y <= footprint; y++) {
+	    for (int x = - footprint; x <= footprint; x++) {
+		double resid = target->kernel[y][x] - background - source->kernel[y][x] * norm + residual->kernel[y][x];
+		double value = PS_SQR(resid);
+		if (weight) {
+		    float wtVal = weight->kernel[y][x];
+		    value *= wtVal;
+		}
+		if (window) {
+		    float  winVal = window->kernel[y][x];
+		    value *= winVal;
+		}
+		sum += value;
+	    }
+	}
     }
     psFree (polyValues);
@@ -1824,5 +1889,5 @@
 
     for (int i = 0; i < w->n; i++) {
-        wApply->data.F64[i] = wMask->data.U8[i] ? 0.0 : w->data.F64[i];
+	wApply->data.F64[i] = wMask->data.U8[i] ? 0.0 : w->data.F64[i];
     }
     return true;
@@ -1839,9 +1904,9 @@
     // generate Vn = V * w^{-1}
     for (int j = 0; j < Vn->numRows; j++) {
-        for (int i = 0; i < Vn->numCols; i++) {
-            if (!isfinite(w->data.F64[i])) continue;
-            if (w->data.F64[i] == 0.0) continue;
-            Vn->data.F64[j][i] = V->data.F64[j][i] / w->data.F64[i];
-        }
+	for (int i = 0; i < Vn->numCols; i++) {
+	    if (!isfinite(w->data.F64[i])) continue;
+	    if (w->data.F64[i] == 0.0) continue;
+	    Vn->data.F64[j][i] = V->data.F64[j][i] / w->data.F64[i];
+	}
     }
 
@@ -1851,11 +1916,11 @@
     // generate Xvar = Vn * Vn^T
     for (int j = 0; j < Vn->numRows; j++) {
-        for (int i = 0; i < Vn->numCols; i++) {
-            double sum = 0.0;
-            for (int k = 0; k < Vn->numCols; k++) {
-                sum += Vn->data.F64[k][i]*Vn->data.F64[k][j];
-            }
-            Xvar->data.F64[j][i] = sum;
-        }
+	for (int i = 0; i < Vn->numCols; i++) {
+	    double sum = 0.0;
+	    for (int k = 0; k < Vn->numCols; k++) {
+		sum += Vn->data.F64[k][i]*Vn->data.F64[k][j];
+	    }
+	    Xvar->data.F64[j][i] = sum;
+	}
     }
     return Xvar;
@@ -1872,5 +1937,5 @@
     psFitsWriteImage(fits, header, image, 0, NULL);
     psFitsClose(fits);
-
+    
     return true;
 }