Index: trunk/psModules/src/imcombine/pmSubtractionEquation.c
===================================================================
--- trunk/psModules/src/imcombine/pmSubtractionEquation.c	(revision 26031)
+++ trunk/psModules/src/imcombine/pmSubtractionEquation.c	(revision 26035)
@@ -17,5 +17,5 @@
 //#define TESTING                         // TESTING output for debugging; may not work with threads!
 
-#define USE_VARIANCE                    // Include variance in equation?
+#define USE_WEIGHT                      // Include weight (1/variance) in equation?
 
 
@@ -29,5 +29,5 @@
                                   const psKernel *input, // Input image (target)
                                   const psKernel *reference, // Reference image (convolution source)
-                                  const psKernel *variance,  // Variance image
+                                  const psKernel *weight,  // Weight image
                                   const psArray *convolutions,         // Convolutions for each kernel
                                   const pmSubtractionKernels *kernels, // Kernels
@@ -74,6 +74,6 @@
                 for (int x = - footprint; x <= footprint; x++) {
                     double cc = iConv->kernel[y][x] * jConv->kernel[y][x];
-#ifdef USE_VARIANCE
-                    cc /= variance->kernel[y][x];
+#ifdef USE_WEIGHT
+                    cc *= weight->kernel[y][x];
 #endif
                     sumCC += cc;
@@ -102,9 +102,9 @@
                 double rc = ref * conv;
                 double c = conv;
-#ifdef USE_VARIANCE
-                float varVal = 1.0 / variance->kernel[y][x];
-                ic *= varVal;
-                rc *= varVal;
-                c *= varVal;
+#ifdef USE_WEIGHT
+                float wtVal = weight->kernel[y][x];
+                ic *= wtVal;
+                rc *= wtVal;
+                c *= wtVal;
 #endif
                 sumIC += ic;
@@ -137,11 +137,11 @@
             double rr = PS_SQR(ref);
             double one = 1.0;
-#ifdef USE_VARIANCE
-            float varVal = 1.0 / variance->kernel[y][x];
-            rr *= varVal;
-            ir *= varVal;
-            in *= varVal;
-            ref *= varVal;
-            one *= varVal;
+#ifdef USE_WEIGHT
+            float wtVal = weight->kernel[y][x];
+            rr *= wtVal;
+            ir *= wtVal;
+            in *= wtVal;
+            ref *= wtVal;
+            one *= wtVal;
 #endif
             sumRR += rr;
@@ -169,5 +169,5 @@
                                       const psKernel *image1, // Image 1
                                       const psKernel *image2, // Image 2
-                                      const psKernel *variance,  // Variance image
+                                      const psKernel *weight,  // Weight image
                                       const psArray *convolutions1, // Convolutions of image 1 for each kernel
                                       const psArray *convolutions2, // Convolutions of image 2 for each kernel
@@ -227,8 +227,9 @@
                     double bb = iConv2->kernel[y][x] * jConv2->kernel[y][x];
                     double ab = iConv1->kernel[y][x] * jConv2->kernel[y][x];
-#ifdef USE_VARIANCE
-                    aa /= variance->kernel[y][x];
-                    bb /= variance->kernel[y][x];
-                    ab /= variance->kernel[y][x];
+#ifdef USE_WEIGHT
+                    float wtVal = weight->kernel[y][x];
+                    aa *= wtVal;
+                    bb *= wtVal;
+                    ab *= wtVal;
 #endif
                     sumAA += aa;
@@ -258,6 +259,6 @@
                 for (int x = - footprint; x <= footprint; x++) {
                     double ab = iConv1->kernel[y][x] * jConv2->kernel[y][x];
-#ifdef USE_VARIANCE
-                    ab /= variance->kernel[y][x];
+#ifdef USE_WEIGHT
+                    ab *= weight->kernel[y][x];
 #endif
                     sumAB += ab;
@@ -295,14 +296,14 @@
                 double i1i2 = i1 * i2;
 
-#ifdef USE_VARIANCE
-                float varVal = 1.0 / variance->kernel[y][x];
-                ai2 *= varVal;
-                bi2 *= varVal;
-                ai1 *= varVal;
-                bi1 *= varVal;
-                i1i2 *= varVal;
-                a *= varVal;
-                b *= varVal;
-                i2 *= varVal;
+#ifdef USE_WEIGHT
+                float wtVal = weight->kernel[y][x];
+                ai2 *= wtVal;
+                bi2 *= wtVal;
+                ai1 *= wtVal;
+                bi1 *= wtVal;
+                i1i2 *= wtVal;
+                a *= wtVal;
+                b *= wtVal;
+                i2 *= wtVal;
 #endif
 
@@ -351,11 +352,11 @@
             double i1i2 = i1 * i2;
 
-#ifdef USE_VARIANCE
-            float varVal = 1.0 / variance->kernel[y][x];
-            i1 *= varVal;
-            i1i1 *= varVal;
-            one *= varVal;
-            i2 *= varVal;
-            i1i2 *= varVal;
+#ifdef USE_WEIGHT
+            float wtVal = weight->kernel[y][x];
+            i1 *= wtVal;
+            i1i1 *= wtVal;
+            one *= wtVal;
+            i2 *= wtVal;
+            i1i2 *= wtVal;
 #endif
 
@@ -619,10 +620,10 @@
       case PM_SUBTRACTION_MODE_1:
         status = calculateMatrixVector(stamp->matrix1, stamp->vector1, stamp->image2, stamp->image1,
-                                       stamp->variance, stamp->convolutions1, kernels, polyValues,
+                                       stamp->weight, stamp->convolutions1, kernels, polyValues,
                                        footprint);
         break;
       case PM_SUBTRACTION_MODE_2:
         status = calculateMatrixVector(stamp->matrix1, stamp->vector1, stamp->image1, stamp->image2,
-                                       stamp->variance, stamp->convolutions2, kernels, polyValues,
+                                       stamp->weight, stamp->convolutions2, kernels, polyValues,
                                        footprint);
         break;
@@ -639,5 +640,5 @@
 #endif
         status = calculateDualMatrixVector(stamp->matrix1, stamp->vector1, stamp->matrix2, stamp->vector2,
-                                           stamp->matrixX, stamp->image1, stamp->image2, stamp->variance,
+                                           stamp->matrixX, stamp->image1, stamp->image2, stamp->weight,
                                            stamp->convolutions1, stamp->convolutions2, kernels, polyValues,
                                            footprint);
@@ -1163,5 +1164,5 @@
 
         // Calculate residuals
-        psKernel *variance = stamp->variance; // Variance postage stamp
+        psKernel *weight = stamp->weight; // Weight postage stamp
         psImageInit(residual->image, 0.0);
         if (kernels->mode != PM_SUBTRACTION_MODE_DUAL) {
@@ -1247,5 +1248,5 @@
         for (int y = - footprint; y <= footprint; y++) {
             for (int x = - footprint; x <= footprint; x++) {
-                double dev = PS_SQR(residual->kernel[y][x]) / variance->kernel[y][x];
+                double dev = PS_SQR(residual->kernel[y][x]) * weight->kernel[y][x];
                 deviation += dev;
 #ifdef TESTING
@@ -1290,10 +1291,10 @@
             psFitsClose(fits);
         }
-        if (stamp->variance) {
+        if (stamp->weight) {
             psString filename = NULL;
-            psStringAppend(&filename, "stamp_variance_%03d.fits", i);
+            psStringAppend(&filename, "stamp_weight_%03d.fits", i);
             psFits *fits = psFitsOpen(filename, "w");
             psFree(filename);
-            psFitsWriteImage(fits, NULL, stamp->variance->image, 0, NULL);
+            psFitsWriteImage(fits, NULL, stamp->weight->image, 0, NULL);
             psFitsClose(fits);
         }