Index: branches/pap/ppSub/src/ppSubReadoutSubtract.c
===================================================================
--- branches/pap/ppSub/src/ppSubReadoutSubtract.c	(revision 27709)
+++ branches/pap/ppSub/src/ppSubReadoutSubtract.c	(revision 27727)
@@ -49,11 +49,5 @@
     outRO->image = (psImage*)psBinaryOp(outRO->image, minuend->image, "-", subtrahend->image);
     outRO->mask = (psImage*)psBinaryOp(outRO->mask, minuend->mask, "|", subtrahend->mask);
-
-    // Combining the variance is not straight-forward.  The problem is, we're representing the variance by a
-    // variance map and a covariance pseudo-matrix, and we cannot simply add them and get the correct result,
-    // because they don't combine nicely.  It's like we want dD = aA + bB, but the upper- and lower-case
-    // variables can't be merged.  Our solution is to pull the normalisation of the covariance pseudo-matrices
-    // into the variance map, sum the scaled variance maps and (weighted) average the covariance
-    // pseudo-matrices.  It's not perfectly correct, but it's gotta be good enough.
+    outRO->variance = (psImage*)psBinaryOp(outRO->variance, minuend->variance, "+", subtrahend->variance);
 
     // Measure the variance scales
@@ -69,39 +63,17 @@
     psFree(rng);
 
-    // Measure and apply the covariance scales
-    float minuendSum = psImageCovarianceSum(minuend->covariance); // Sum of covariance for minuend
-    float subtrahendSum = psImageCovarianceSum(subtrahend->covariance); // Sum of covariance for subtrahend
-    psBinaryOp(minuend->covariance->image, minuend->covariance->image, "/",
-               psScalarAlloc(minuendSum, PS_TYPE_F32));
-    psBinaryOp(subtrahend->covariance->image, subtrahend->covariance->image, "/",
-               psScalarAlloc(subtrahendSum, PS_TYPE_F32));
-
-    // Combine the scaled covariances
+    // Combine the covariances
+    // These are weighted by the appropriate mean variance.  This is probably not perfectly correct, but it
+    // does seem to reproduce the correct magnitude limit in psphot.
     psArray *covars = psArrayAlloc(2);  // Covariance pseudo-matrices
     psVector *covarWeights = psVectorAlloc(2, PS_TYPE_F32); // Weights for covariances
     covars->data[0] = psMemIncrRefCounter(minuend->covariance);
     covars->data[1] = psMemIncrRefCounter(subtrahend->covariance);
-    covarWeights->data.F32[0] = minuendSum * minuendVar;
-    covarWeights->data.F32[1] = subtrahendSum * subtrahendVar;
+    covarWeights->data.F32[0] = minuendVar;
+    covarWeights->data.F32[1] = subtrahendVar;
     outRO->covariance = psImageCovarianceAverageWeighted(covars, covarWeights);
     psFree(covars);
     psFree(covarWeights);
 
-    // Restore the scales, in case someone cares to look at the noise in the convolved images
-    psBinaryOp(minuend->covariance->image, minuend->covariance->image, "*",
-               psScalarAlloc(minuendSum, PS_TYPE_F32));
-    psBinaryOp(subtrahend->covariance->image, subtrahend->covariance->image, "*",
-               psScalarAlloc(subtrahendSum, PS_TYPE_F32));
-
-    // Combine the variance maps
-    int numCols = outRO->image->numCols, numRows = outRO->image->numRows; // Size of image
-    psImage *outVar = outRO->variance = psImageAlloc(numCols, numRows, PS_TYPE_F32); // Output variance map
-    psImage *minuendVariance = minuend->variance, *subtrahendVariance = subtrahend->variance; // Variance maps
-    for (int y = 0; y < numRows; y++) {
-        for (int x = 0; x < numCols; x++) {
-            outVar->data.F32[y][x] = minuendVariance->data.F32[y][x] * minuendSum +
-                subtrahendVariance->data.F32[y][x] * subtrahendSum;
-        }
-    }
     psImageCovarianceTransfer(outRO->variance, outRO->covariance);