Index: trunk/psphot/src/psphotChoosePSF.c
===================================================================
--- trunk/psphot/src/psphotChoosePSF.c	(revision 20216)
+++ trunk/psphot/src/psphotChoosePSF.c	(revision 20234)
@@ -318,6 +318,4 @@
 bool psphotPSFstats (pmReadout *readout, psMetadata *recipe, pmPSF *psf) {
 
-    double FWHM_X, FWHM_Y;
-
     psEllipseShape shape;
     psEllipseAxes axes;
@@ -330,41 +328,56 @@
     PS_ASSERT_PTR_NON_NULL(image, false);
 
-    // XXXX A Serious hack: the psf might not be determined in the field center.
-    // for the moment, just try a few locations until it is defined!
-
-    pmModel *modelPSF = NULL;
-    for (int ix = -1; ix <= +1; ix ++) {
-	for (int iy = -1; iy <= +1; iy ++) {
+    // XXX a minor hack: measure the values for a grid of points across the image, determine mean, UQ, LQ, etc:
+    psVector *fwhmMajor = psVectorAllocEmpty (100, PS_DATA_F32);
+    psVector *fwhmMinor = psVectorAllocEmpty (100, PS_DATA_F32);
+
+    for (float ix = -0.4; ix <= +0.4; ix += 0.1) {
+	for (float iy = -0.4; iy <= +0.4; iy += 0.1) {
 
 	    // use the center of the center pixel of the image
-	    float xc = (0.5 + 0.3*ix)*image->numCols + image->col0 + 0.5;
-	    float yc = (0.5 + 0.3*ix)*image->numRows + image->row0 + 0.5;
+	    float xc = ix*image->numCols + 0.5*image->numCols + image->col0 + 0.5;
+	    float yc = iy*image->numRows + 0.5*image->numRows + image->row0 + 0.5;
 
 	    // create modelPSF from this model
-	    modelPSF = pmModelFromPSFforXY (psf, xc, yc, 1.0);
-	    if (modelPSF) goto got_model;
+	    pmModel *modelPSF = pmModelFromPSFforXY (psf, xc, yc, 1.0);
+	    if (!modelPSF) continue;
+
+	    // get the model full-width at half-max
+	    float FWHM_MAJOR = 2*modelPSF->modelRadius (modelPSF->params, 0.5);
+
+	    // XXX make sure this is consistent with the re-definition of PM_PAR_SXX
+	    shape.sx  = modelPSF->params->data.F32[PM_PAR_SXX];
+	    shape.sy  = modelPSF->params->data.F32[PM_PAR_SYY];
+	    shape.sxy = modelPSF->params->data.F32[PM_PAR_SXY];
+	    axes = psEllipseShapeToAxes (shape, 20.0);
+	    psFree (modelPSF);
+
+	    float FWHM_MINOR = FWHM_MAJOR * (axes.minor / axes.major);
+	    if (!isfinite(FWHM_MAJOR) || !isfinite(FWHM_MINOR)) continue;
+	    psVectorAppend (fwhmMajor, FWHM_MAJOR);
+	    psVectorAppend (fwhmMinor, FWHM_MINOR);
 	}
     }
-    psAssert (modelPSF, "Failed to estimate PSF model at image center (psf is invalid everywhere?)");
-
-got_model:
-    // get the model full-width at half-max
-    FWHM_X = 2*modelPSF->modelRadius (modelPSF->params, 0.5);
-
-    // XXX make sure this is consistent with the re-definition of PM_PAR_SXX
-    shape.sx  = modelPSF->params->data.F32[PM_PAR_SXX];
-    shape.sy  = modelPSF->params->data.F32[PM_PAR_SYY];
-    shape.sxy = modelPSF->params->data.F32[PM_PAR_SXY];
-    axes = psEllipseShapeToAxes (shape, 20.0);
-
-    FWHM_Y = FWHM_X * (axes.minor / axes.major);
-
-    psMetadataAddF32 (recipe, PS_LIST_TAIL, "FWHM_X",   PS_META_REPLACE, "PSF FWHM Major axis", FWHM_X);
-    psMetadataAddF32 (recipe, PS_LIST_TAIL, "FWHM_Y",   PS_META_REPLACE, "PSF FWHM Minor axis", FWHM_Y);
+
+    psStats *stats = psStatsAlloc (PS_STAT_SAMPLE_MEAN | PS_STAT_SAMPLE_STDEV | PS_STAT_SAMPLE_QUARTILE);
+    psVectorStats (stats, fwhmMajor, NULL, NULL, 0);
+    psMetadataAddF32 (recipe, PS_LIST_TAIL, "FWHM_MAJ",   PS_META_REPLACE, "PSF FWHM Major axis (mean)", stats->sampleMean);
+    psMetadataAddF32 (recipe, PS_LIST_TAIL, "FW_MJ_SG",   PS_META_REPLACE, "PSF FWHM Major axis (sigma)", stats->sampleStdev);
+    psMetadataAddF32 (recipe, PS_LIST_TAIL, "FW_MJ_LQ",   PS_META_REPLACE, "PSF FWHM Major axis (upper quartile)", stats->sampleLQ);
+    psMetadataAddF32 (recipe, PS_LIST_TAIL, "FW_MJ_UQ",   PS_META_REPLACE, "PSF FWHM Major axis (lower quartile)", stats->sampleUQ);
+
+    psVectorStats (stats, fwhmMinor, NULL, NULL, 0);
+    psMetadataAddF32 (recipe, PS_LIST_TAIL, "FWHM_MIN",   PS_META_REPLACE, "PSF FWHM Minor axis (mean)", stats->sampleMean);
+    psMetadataAddF32 (recipe, PS_LIST_TAIL, "FW_MN_SG",   PS_META_REPLACE, "PSF FWHM Minor axis (sigma)", stats->sampleStdev);
+    psMetadataAddF32 (recipe, PS_LIST_TAIL, "FW_MN_LQ",   PS_META_REPLACE, "PSF FWHM Minor axis (upper quartile)", stats->sampleLQ);
+    psMetadataAddF32 (recipe, PS_LIST_TAIL, "FW_MN_UQ",   PS_META_REPLACE, "PSF FWHM Minor axis (lower quartile)", stats->sampleUQ);
+
     psMetadataAddF32 (recipe, PS_LIST_TAIL, "ANGLE",    PS_META_REPLACE, "PSF angle",           axes.theta);
     psMetadataAddS32 (recipe, PS_LIST_TAIL, "NPSFSTAR", PS_META_REPLACE, "Number of stars used to make PSF", psf->nPSFstars);
     psMetadataAddBool(recipe, PS_LIST_TAIL, "PSFMODEL", PS_META_REPLACE, "Valid PSF Model?", true);
-    psFree (modelPSF);
-
+
+    psFree (fwhmMajor);
+    psFree (fwhmMinor);
+    psFree (stats);
     return true;
 }