Index: /branches/eam_branches/20091201/psModules/src/imcombine/pmSubtractionEquation.c
===================================================================
--- /branches/eam_branches/20091201/psModules/src/imcombine/pmSubtractionEquation.c	(revision 26608)
+++ /branches/eam_branches/20091201/psModules/src/imcombine/pmSubtractionEquation.c	(revision 26609)
@@ -15,5 +15,5 @@
 #include "pmSubtractionVisual.h"
 
-// #define TESTING                         // TESTING output for debugging; may not work with threads!
+//#define TESTING                         // TESTING output for debugging; may not work with threads!
 
 //#define USE_WEIGHT                      // Include weight (1/variance) in equation?
@@ -487,5 +487,5 @@
                 // Contribution to chi^2: a_i^2 P_i
                 psAssert(isfinite(penalties->data.F32[i]), "Invalid penalty");
-                matrix->data.F64[index][index] -= norm * penalties->data.F32[i];
+                matrix->data.F64[index][index] += norm * penalties->data.F32[i];
             }
         }
@@ -1337,8 +1337,8 @@
     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));
-	}
+        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));
+        }
     }
 
@@ -1353,13 +1353,13 @@
 
     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 ++;
-	}
+        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));
@@ -1398,29 +1398,29 @@
     psVector *keepStamps  = psVectorAlloc(stamps->num, PS_TYPE_S32);
     psVectorInit (keepStamps, 0);
-    { 
-	psVector *flux  = psVectorAlloc(stamps->num, PS_TYPE_F32);
-	psVectorInit (flux, 0.0);
-	
-	for (int i = 0; i < stamps->num; i++) {
-	    pmSubtractionStamp *stamp = stamps->stamps->data[i];
-	    if (!isfinite(stamp->flux)) continue;
-	    flux->data.F32[i] = stamp->flux;
-	}
-
-	psVector *index = psVectorSortIndex(NULL, flux);
-	for (int i = 0; (i < stamps->num) && (i < 9); i++) {
-	    int n = stamps->num - i - 1;
-	    keepStamps->data.S32[index->data.S32[n]] = 1;
-	    fprintf (stderr, "keeping %d (%d of %d)\n", index->data.S32[n], n, 9);
-	}
-	psFree (flux);
-	psFree (index);
-
-	// this function is called multiple times in the iteration, but 
-	// we only know after the interation is done if we will try again.
-	// therefore we must save the sample each time, and blow away the old one
-	// if it exists.
-	psFree (kernels->sampleStamps);
-	kernels->sampleStamps = psArrayAllocEmpty(9);
+    {
+        psVector *flux  = psVectorAlloc(stamps->num, PS_TYPE_F32);
+        psVectorInit (flux, 0.0);
+
+        for (int i = 0; i < stamps->num; i++) {
+            pmSubtractionStamp *stamp = stamps->stamps->data[i];
+            if (!isfinite(stamp->flux)) continue;
+            flux->data.F32[i] = stamp->flux;
+        }
+
+        psVector *index = psVectorSortIndex(NULL, flux);
+        for (int i = 0; (i < stamps->num) && (i < 9); i++) {
+            int n = stamps->num - i - 1;
+            keepStamps->data.S32[index->data.S32[n]] = 1;
+            fprintf (stderr, "keeping %d (%d of %d)\n", index->data.S32[n], n, 9);
+        }
+        psFree (flux);
+        psFree (index);
+
+        // this function is called multiple times in the iteration, but
+        // we only know after the interation is done if we will try again.
+        // therefore we must save the sample each time, and blow away the old one
+        // if it exists.
+        psFree (kernels->sampleStamps);
+        kernels->sampleStamps = psArrayAllocEmpty(9);
     }
 
@@ -1498,11 +1498,11 @@
             }
 
-	    if (keepStamps->data.S32[i]) {
-		psImage *sample = psImageCopy(NULL, residual->image, PS_TYPE_F32);
-		psArrayAdd (kernels->sampleStamps, 9, sample);
-		psFree (sample);
-	    }
-
-	    pmSubtractionResidualStats(fSigRes, fMaxRes, fMinRes, target, source, residual, norm, footprint);
+            if (keepStamps->data.S32[i]) {
+                psImage *sample = psImageCopy(NULL, residual->image, PS_TYPE_F32);
+                psArrayAdd (kernels->sampleStamps, 9, sample);
+                psFree (sample);
+            }
+
+            pmSubtractionResidualStats(fSigRes, fMaxRes, fMinRes, target, source, residual, norm, footprint);
 
         } else {
@@ -1534,11 +1534,11 @@
                 }
             }
-	    if (keepStamps->data.S32[i]) {
-		psImage *sample = psImageCopy(NULL, residual->image, PS_TYPE_F32);
-		psArrayAdd (kernels->sampleStamps, 9, sample);
-		psFree (sample);
-	    }
-
-	    pmSubtractionResidualStats(fSigRes, fMaxRes, fMinRes, image1, image2, residual, norm, footprint);
+            if (keepStamps->data.S32[i]) {
+                psImage *sample = psImageCopy(NULL, residual->image, PS_TYPE_F32);
+                psArrayAdd (kernels->sampleStamps, 9, sample);
+                psFree (sample);
+            }
+
+            pmSubtractionResidualStats(fSigRes, fMaxRes, fMinRes, image1, image2, residual, norm, footprint);
         }