Index: branches/tap_branches/psphot/src/psphotFitSourcesLinear.c
===================================================================
--- branches/tap_branches/psphot/src/psphotFitSourcesLinear.c	(revision 25900)
+++ branches/tap_branches/psphot/src/psphotFitSourcesLinear.c	(revision 27838)
@@ -12,5 +12,44 @@
 static bool SetBorderMatrixElements (psSparseBorder *border, pmReadout *readout, psArray *sources, bool constant_weights, int SKY_FIT_ORDER, psImageMaskType markVal);
 
-bool psphotFitSourcesLinear (pmReadout *readout, psArray *sources, const psMetadata *recipe, const pmPSF *psf, bool final) {
+// for now, let's store the detections on the readout->analysis for each readout
+bool psphotFitSourcesLinear (pmConfig *config, const pmFPAview *view, bool final)
+{
+    bool status = true;
+
+    // select the appropriate recipe information
+    psMetadata *recipe  = psMetadataLookupPtr (&status, config->recipes, PSPHOT_RECIPE);
+    assert (recipe);
+
+    int num = psMetadataLookupS32 (&status, config->arguments, "PSPHOT.INPUT.NUM");
+    psAssert (status, "programming error: must define PSPHOT.INPUT.NUM");
+
+    // loop over the available readouts
+    for (int i = 0; i < num; i++) {
+
+	// find the currently selected readout
+	pmFPAfile *file = pmFPAfileSelectSingle(config->files, "PSPHOT.INPUT", i); // File of interest
+	psAssert (file, "missing file?");
+
+	pmReadout *readout = pmFPAviewThisReadout(view, file->fpa);
+	psAssert (readout, "missing readout?");
+
+	pmDetections *detections = psMetadataLookupPtr (&status, readout->analysis, "PSPHOT.DETECTIONS");
+	psAssert (detections, "missing detections?");
+
+	psArray *sources = detections->allSources;
+	psAssert (sources, "missing sources?");
+
+	pmPSF *psf = psMetadataLookupPtr (&status, readout->analysis, "PSPHOT.PSF");
+	psAssert (psf, "missing psf?");
+
+	if (!psphotFitSourcesLinearReadout (recipe, readout, sources, psf, final)) {
+            psError (PSPHOT_ERR_CONFIG, false, "failed to fit sources (linear) for PSPHOT.INPUT entry %d", i);
+	    return false;
+	}
+    }
+    return true;
+}
+
+bool psphotFitSourcesLinearReadout (psMetadata *recipe, pmReadout *readout, psArray *sources, pmPSF *psf, bool final) {
 
     bool status;
@@ -18,5 +57,9 @@
     float y;
     float f;
-    // float r;
+
+    if (!sources->n) {
+	psLogMsg ("psphot", PS_LOG_INFO, "no sources, skipping linear fit");
+	return true;
+    }
 
     psTimerStart ("psphot.linear");
@@ -59,4 +102,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 +170,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 +183,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 +189,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 +215,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 +269,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,5 +289,10 @@
 
     psphotVisualShowResidualImage (readout);
+    psphotVisualPlotChisq (sources);
     // psphotVisualShowFlags (sources);
+
+    // We have to place this visualization here because the models are not realized until
+    // psphotGuessModels or fitted until psphotFitSourcesLinear.
+    psphotVisualShowPSFStars (recipe, psf, sources);
 
     return true;
@@ -244,4 +300,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