Index: /trunk/psphot/src/psphotChoosePSF.c
===================================================================
--- /trunk/psphot/src/psphotChoosePSF.c	(revision 10201)
+++ /trunk/psphot/src/psphotChoosePSF.c	(revision 10202)
@@ -26,4 +26,10 @@
     }
 
+    float PSF_MIN_DS = psMetadataLookupF32 (&status, recipe, "PSF_MIN_DS");
+    if (!status) {
+	PSF_MIN_DS = 0.01;
+        psWarning("PSF_MIN_DS is not set --- defaulting to %f\n", PSF_MIN_DS);
+    }
+
     // supply the measured sky variance for optional constant errors (non-poissonian)
     float SKY_STDEV = psMetadataLookupF32 (&status, recipe, "SKY_STDEV");
@@ -39,4 +45,11 @@
     }
     pmSourceFitModelInit (15, 0.01, PS_SQR(SKY_STDEV), POISSON_ERRORS);
+
+    // use poissonian errors or local-sky errors
+    bool PSF_PARAM_WEIGHTS = psMetadataLookupBool (&status, recipe, "PSF_PARAM_WEIGHTS");
+    if (!status) {
+        PSF_PARAM_WEIGHTS = false;
+        psWarning("PSF_PARAM_WEIGHTS is not set in the recipe --- defaulting to false.\n");
+    }
 
     // how to model the PSF variations across the field
@@ -98,5 +111,5 @@
         psMetadataItem *item = psListGetAndIncrement (iter);
         modelName = item->data.V;
-        models->data[i] = pmPSFtryModel (stars, modelName, RADIUS, POISSON_ERRORS, psfTrendMask);
+        models->data[i] = pmPSFtryModel (stars, modelName, RADIUS, POISSON_ERRORS, psfTrendMask, PSF_PARAM_WEIGHTS);
     }
 
@@ -132,4 +145,70 @@
     try = models->data[bestN];
 
+    // measure and save the median value of dparams[PM_PAR_SXX],dparams[PM_PAR_SYY]
+    // these are used by psphotEvalPSF to identify extended sources
+    psVector *Sx = psVectorAllocEmpty (try->sources->n, PS_TYPE_F32);
+    psVector *Sy = psVectorAllocEmpty (try->sources->n, PS_TYPE_F32);
+    psVector *dSx = psVectorAllocEmpty (try->sources->n, PS_TYPE_F32);
+    psVector *dSy = psVectorAllocEmpty (try->sources->n, PS_TYPE_F32);
+    psVector *dSN = psVectorAllocEmpty (try->sources->n, PS_TYPE_F32);
+    for (int i = 0; i < try->sources->n; i++) {
+	// masked for: bad model fit, outlier in parameters
+	if (try->mask->data.U8[i] & PSFTRY_MASK_ALL)
+	    continue;
+
+	pmSource *source = try->sources->data[i];
+	Sx->data.F32[Sx->n] = source->modelPSF->params->data.F32[PM_PAR_SXX];
+	Sy->data.F32[Sy->n] = source->modelPSF->params->data.F32[PM_PAR_SYY];
+	dSN->data.F32[dSN->n] = source->modelPSF->dparams->data.F32[PM_PAR_I0] / source->modelPSF->params->data.F32[PM_PAR_I0];
+	dSx->data.F32[dSx->n] = source->modelPSF->dparams->data.F32[PM_PAR_SXX];
+	dSy->data.F32[dSy->n] = source->modelPSF->dparams->data.F32[PM_PAR_SYY];
+	Sx->n ++;
+	Sy->n ++;
+	dSN->n ++;
+	dSx->n ++;
+	dSy->n ++;
+    }
+    psStats *stats = psStatsAlloc (PS_STAT_SAMPLE_MEAN | PS_STAT_SAMPLE_MEDIAN | PS_STAT_SAMPLE_STDEV);
+
+    psVectorStats (stats, dSx, NULL, NULL, 0);
+    float dSxo = stats->sampleMean;
+    psLogMsg ("psphot.choosePSF", 3, "dSx: mean: %f, median: %f, stdev: %f, npts: %ld", stats->sampleMean, stats->sampleMedian, stats->sampleStdev, dSx->n);
+
+    psVectorStats (stats, dSy, NULL, NULL, 0);
+    float dSyo = stats->sampleMean;
+    psLogMsg ("psphot.choosePSF", 3, "dSy: mean: %f, median: %f, stdev: %f, npts: %ld", stats->sampleMean, stats->sampleMedian, stats->sampleStdev, dSy->n);
+
+    for (int i = 0; i < dSN->n; i++) {
+	dSx->data.F32[i] = (dSx->data.F32[i] - dSxo) / PS_MAX (PSF_MIN_DS, Sx->data.F32[i]*dSN->data.F32[i]);
+	dSy->data.F32[i] = (dSy->data.F32[i] - dSyo) / PS_MAX (PSF_MIN_DS, Sy->data.F32[i]*dSN->data.F32[i]);
+    }
+
+    psVectorStats (stats, dSx, NULL, NULL, 0);
+    float nSx = stats->sampleStdev;
+    psLogMsg ("psphot.choosePSF", 3, "dSx: mean: %f, median: %f, stdev: %f, npts: %ld", stats->sampleMean, stats->sampleMedian, stats->sampleStdev, dSx->n);
+    psVectorStats (stats, dSy, NULL, NULL, 0);
+    float nSy = stats->sampleStdev;
+    psLogMsg ("psphot.choosePSF", 3, "dSy: mean: %f, median: %f, stdev: %f, npts: %ld", stats->sampleMean, stats->sampleMedian, stats->sampleStdev, dSy->n);
+
+    psFree (Sx);
+    psFree (Sy);
+    psFree (dSx);
+    psFree (dSy);
+    psFree (dSN);
+    psFree (stats);
+
+    psMetadataAddF32 (recipe, PS_LIST_TAIL, "DSX_MEAN", PS_META_REPLACE, "mean offset of dSXX", dSxo);
+    psMetadataAddF32 (recipe, PS_LIST_TAIL, "DSY_MEAN", PS_META_REPLACE, "mean offset of dSYY", dSyo);
+    psMetadataAddF32 (recipe, PS_LIST_TAIL, "DSX_STDV", PS_META_REPLACE, "stdev of mean-corrected dSXX", nSx);
+    psMetadataAddF32 (recipe, PS_LIST_TAIL, "DSY_STDV", PS_META_REPLACE, "stdev of mean-corrected dSYY", nSy);
+
+    // unset the PSFSTAR flag for stars not used for PSF model
+    for (int i = 0; i < try->sources->n; i++) {
+        pmSource *source = try->sources->data[i];
+        if (try->mask->data.U8[i]) {
+            source->mode &= ~PM_SOURCE_MODE_PSFSTAR;
+        }
+    }
+
     // XXX test dump of psf star data and psf-subtracted image
     if (psTraceGetLevel("psphot.psfstars") > 5) { 
@@ -159,12 +238,4 @@
     }
 
-    // unset the PSFSTAR flag for stars not used for PSF model
-    for (int i = 0; i < try->sources->n; i++) {
-        pmSource *source = try->sources->data[i];
-        if (try->mask->data.U8[i] & PSFTRY_MASK_EXT_FAIL) {
-            source->mode &= ~PM_SOURCE_MODE_PSFSTAR;
-        }
-    }
-
     // save only the best model;
     // XXX we are not saving the fitted sources