Index: trunk/psModules/src/imcombine/pmStack.c
===================================================================
--- trunk/psModules/src/imcombine/pmStack.c	(revision 20487)
+++ trunk/psModules/src/imcombine/pmStack.c	(revision 20497)
@@ -8,6 +8,6 @@
  *  @author GLG, MHPCC
  *
- *  @version $Revision: 1.43 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2008-10-31 21:50:59 $
+ *  @version $Revision: 1.44 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2008-11-01 02:59:33 $
  *  Copyright 2004-2007 Institute for Astronomy, University of Hawaii
  *
@@ -219,4 +219,5 @@
                           int numIter, // Number of rejection iterations
                           float rej, // Number of standard deviations at which to reject
+                          float sys,    // Relative systematic error
                           bool useVariance, // Use variance for rejection when combining?
                           bool safe,    // Combine safely?
@@ -311,11 +312,17 @@
           if (useVariance && numIter > 0) {
               // Use variance to check that the two are consistent
-              float diff = pixelData->data.F32[0] - pixelData->data.F32[1];
+              float diff = 0.5 * (pixelData->data.F32[0] - pixelData->data.F32[1]); // Mean flux difference
+              float var1 = pixelVariances->data.F32[0]; // Variance of first
+              float var2 = pixelVariances->data.F32[1]; // Variance of second
 #ifdef VARIANCE_FACTORS
-              float sigma2 = pixelVariances->data.F32[0] * varFactors->data.F32[pixelSources->data.U16[0]] +
-                  pixelVariances->data.F32[1] * varFactors->data.F32[pixelSources->data.U16[1]];
-#else
-              float sigma2 = pixelVariances->data.F32[0] + pixelVariances->data.F32[1];
+              // Variance factor contributes to the rejection level
+              var1 *= varFactors->data.F32[pixelSources->data.U16[0]];
+              var2 *= varFactors->data.F32[pixelSources->data.U16[1]];
 #endif
+              // Systematic error contributes to the rejection level
+              var1 += PS_SQR(sys * pixelData->data.F32[0]);
+              var2 += PS_SQR(sys * pixelData->data.F32[1]);
+
+              float sigma2 = var1 + var2; //
               if (PS_SQR(diff) > PS_SQR(rej) * sigma2) {
                   // Not consistent: mark both for inspection
@@ -346,9 +353,11 @@
                   float rej2 = PS_SQR(rej); // Rejection level squared
                   for (int i = 0; i < num; i++) {
+                      pixelVariances->data.F32[i] *= rej2;
 #ifdef VARIANCE_FACTORS
-                      pixelVariances->data.F32[i] *= rej2 * varFactors->data.F32[pixelSources->data.U16[i]];
-#else
-                      pixelVariances->data.F32[i] *= rej2;
+                      // Variance factor contributes to the rejection level
+                      pixelVariances->data.F32[i] *= varFactors->data.F32[pixelSources->data.U16[i]];
 #endif
+                      // Systematic error contributes to the rejection level
+                      pixelVariances->data.F32[i] += PS_SQR(sys * pixelData->data.F32[i]);
                   }
               }
@@ -554,5 +563,6 @@
 /// Stack input images
 bool pmStackCombine(pmReadout *combined, psArray *input, psMaskType maskVal, psMaskType bad,
-                    int kernelSize, int numIter, float rej, bool entire, bool useVariance, bool safe)
+                    int kernelSize, int numIter, float rej, float sys,
+                    bool entire, bool useVariance, bool safe)
 {
     PS_ASSERT_PTR_NON_NULL(combined, false);
@@ -660,5 +670,5 @@
                                                      x, y); // Reject these images
                     combinePixels(combinedImage, combinedMask, combinedVariance, input, weights, varFactors,
-                                  reject, x, y, maskVal, bad, numIter, rej, useVariance, safe, buffer);
+                                  reject, x, y, maskVal, bad, numIter, rej, sys, useVariance, safe, buffer);
                 }
             }
@@ -674,5 +684,5 @@
                                                  x, y); // Reject these images
                 combinePixels(combinedImage, combinedMask, combinedVariance, input, weights, varFactors,
-                              reject, x, y, maskVal, bad, numIter, rej, useVariance, safe, buffer);
+                              reject, x, y, maskVal, bad, numIter, rej, sys, useVariance, safe, buffer);
             }
         }
@@ -701,5 +711,5 @@
             for (int x = minInputCols; x < maxInputCols; x++) {
                 combinePixels(combinedImage, combinedMask, combinedVariance, input, weights, varFactors,
-                              NULL, x, y, maskVal, bad, numIter, rej, useVariance, safe, buffer);
+                              NULL, x, y, maskVal, bad, numIter, rej, sys, useVariance, safe, buffer);
             }
         }