Index: trunk/ppSub/src/ppSubReadoutSubtract.c
===================================================================
--- trunk/ppSub/src/ppSubReadoutSubtract.c	(revision 26982)
+++ trunk/ppSub/src/ppSubReadoutSubtract.c	(revision 28006)
@@ -48,14 +48,33 @@
     pmReadout *outRO = pmFPAfileThisReadout(config->files, view, "PPSUB.OUTPUT");
     outRO->image = (psImage*)psBinaryOp(outRO->image, minuend->image, "-", subtrahend->image);
-    if (minuend->variance && subtrahend->variance) {
-        outRO->variance = (psImage*)psBinaryOp(outRO->variance, minuend->variance, "+", subtrahend->variance);
-    }
     outRO->mask = (psImage*)psBinaryOp(outRO->mask, minuend->mask, "|", subtrahend->mask);
+    outRO->variance = (psImage*)psBinaryOp(outRO->variance, minuend->variance, "+", subtrahend->variance);
 
+    // Measure the variance scales
+    psStats *varStats = psStatsAlloc(PS_STAT_ROBUST_MEDIAN); // Statistics for variance images
+    psRandom *rng = psRandomAlloc(PS_RANDOM_TAUS);           // Random number generator
+    psImageMaskType maskVal = pmConfigMaskGet("MASK.VALUE", config) |
+        pmConfigMaskGet("BLANK", config); // Bits to mask in inputs
+    psImageBackground(varStats, NULL, minuend->variance, minuend->mask, maskVal, rng);
+    float minuendVar = varStats->robustMedian; // Mean variance for minuend
+    psImageBackground(varStats, NULL, subtrahend->variance, subtrahend->mask, maskVal, rng);
+    float subtrahendVar = varStats->robustMedian; // Mean variance for subtrahend
+    psFree(varStats);
+    psFree(rng);
+
+    // 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);
-    outRO->covariance = psImageCovarianceAverage(covars);
+    covarWeights->data.F32[0] = minuendVar;
+    covarWeights->data.F32[1] = subtrahendVar;
+    outRO->covariance = psImageCovarianceAverageWeighted(covars, covarWeights);
     psFree(covars);
+    psFree(covarWeights);
+
+    psImageCovarianceTransfer(outRO->variance, outRO->covariance);
 
     outRO->data_exists = true;