Index: /trunk/psphot/src/psphotFitSourcesLinear.c
===================================================================
--- /trunk/psphot/src/psphotFitSourcesLinear.c	(revision 25989)
+++ /trunk/psphot/src/psphotFitSourcesLinear.c	(revision 25990)
@@ -59,4 +59,11 @@
     }
 
+    // XXX test: choose a larger-than expected radius:
+    float covarFactor = psImageCovarianceFactorForAperture(readout->covariance, 10.0); // Covariance matrix
+    psLogMsg ("psphot.ensemble", PS_LOG_MINUTIA, "covariance factor: %f\n", covarFactor);
+
+    // XXX do not apply covarFactor for the moment...
+    // covarFactor = 1.0;
+
     // select the sources which will be used for the fitting analysis
     for (int i = 0; i < sources->n; i++) {
@@ -120,10 +127,10 @@
 
         // diagonal elements of the sparse matrix (auto-cross-product)
-        f = pmSourceModelDotModel (SRCi, SRCi, CONSTANT_PHOTOMETRIC_WEIGHTS);
+        f = pmSourceModelDotModel (SRCi, SRCi, CONSTANT_PHOTOMETRIC_WEIGHTS, covarFactor);
         psSparseMatrixElement (sparse, i, i, f);
 
         // the formal error depends on the weighting scheme
         if (CONSTANT_PHOTOMETRIC_WEIGHTS) {
-            float var = pmSourceModelDotModel (SRCi, SRCi, false);
+            float var = pmSourceModelDotModel (SRCi, SRCi, false, covarFactor);
             errors->data.F32[i] = 1.0 / sqrt(var);
         } else {
@@ -133,5 +140,5 @@
 
         // find the image x model value
-        f = pmSourceDataDotModel (SRCi, SRCi, CONSTANT_PHOTOMETRIC_WEIGHTS);
+        f = pmSourceDataDotModel (SRCi, SRCi, CONSTANT_PHOTOMETRIC_WEIGHTS, covarFactor);
         psSparseVectorElement (sparse, i, f);
 
@@ -139,11 +146,11 @@
         switch (SKY_FIT_ORDER) {
           case 1:
-            f = pmSourceModelWeight (SRCi, 1, CONSTANT_PHOTOMETRIC_WEIGHTS);
+            f = pmSourceModelWeight (SRCi, 1, CONSTANT_PHOTOMETRIC_WEIGHTS, covarFactor);
             psSparseBorderElementB (border, i, 1, f);
-            f = pmSourceModelWeight (SRCi, 2, CONSTANT_PHOTOMETRIC_WEIGHTS);
+            f = pmSourceModelWeight (SRCi, 2, CONSTANT_PHOTOMETRIC_WEIGHTS, covarFactor);
             psSparseBorderElementB (border, i, 2, f);
 
           case 0:
-            f = pmSourceModelWeight (SRCi, 0, CONSTANT_PHOTOMETRIC_WEIGHTS);
+            f = pmSourceModelWeight (SRCi, 0, CONSTANT_PHOTOMETRIC_WEIGHTS, covarFactor);
             psSparseBorderElementB (border, i, 0, f);
             break;
@@ -165,5 +172,5 @@
 
             // got an overlap; calculate cross-product and add to output array
-            f = pmSourceModelDotModel (SRCi, SRCj, CONSTANT_PHOTOMETRIC_WEIGHTS);
+            f = pmSourceModelDotModel (SRCi, SRCj, CONSTANT_PHOTOMETRIC_WEIGHTS, covarFactor);
             psSparseMatrixElement (sparse, j, i, f);
         }
@@ -219,9 +226,10 @@
 
     // measure chisq for each source
-    for (int i = 0; final && (i < fitSources->n); i++) {
+    // for (int i = 0; final && (i < fitSources->n); i++) {
+    for (int i = 0; i < fitSources->n; i++) {
         pmSource *source = fitSources->data[i];
         if (source->mode & PM_SOURCE_MODE_NONLINEAR_FIT) continue;
         pmModel *model = pmSourceGetModel (NULL, source);
-        pmSourceChisq (model, source->pixels, source->maskObj, source->variance, maskVal);
+        pmSourceChisq (model, source->pixels, source->maskObj, source->variance, maskVal, covarFactor);
     }
     psLogMsg ("psphot.ensemble", PS_LOG_MINUTIA, "get chisqs: %f sec (%d elements)\n", psTimerMark ("psphot.linear"), sparse->Nelem);
@@ -238,4 +246,5 @@
 
     psphotVisualShowResidualImage (readout);
+    psphotVisualPlotChisq (sources);
     // psphotVisualShowFlags (sources);
 
@@ -244,4 +253,5 @@
 
 // XXX do we need this?
+// XXX disallow the simultaneous sky fit and remove this code...
 
 // Calculate the weight terms for the sky fit component of the matrix.  This function operates