Index: trunk/psphot/src/psphotChoosePSF.c
===================================================================
--- trunk/psphot/src/psphotChoosePSF.c	(revision 6964)
+++ trunk/psphot/src/psphotChoosePSF.c	(revision 7758)
@@ -29,5 +29,4 @@
 
     stars = psArrayAlloc (sources->n);
-    // DROP stars->n = 0;
 
     // select the candidate PSF stars (pointers to original sources)
@@ -105,7 +104,4 @@
     psLogMsg ("psphot.pspsf", 3, "selected psf model %s, ApResid: %f +/- %f\n", modelName, psf->ApResid, psf->dApResid);
 
-    char *breakPt = psMetadataLookupStr (&status, recipe, "BREAK_POINT");
-    if (!strcasecmp (breakPt, "PSFFIT")) exit (0);
-
     return (psf);
 }
@@ -151,2 +147,44 @@
     return true;
 }
+
+bool psphotMomentsStats (pmReadout *readout, psMetadata *recipe, psArray *sources) {
+
+    // without the PSF model, we can only rely on the source->moments 
+    // as a measure of the FWHM values
+
+    double FWHM_X, FWHM_Y, FWHM_T;
+    int FWHM_N;
+
+    psEllipseMoments moments;
+    psEllipseAxes axes;
+
+    FWHM_N = 0;
+    FWHM_X = FWHM_Y = 0.0;
+
+    for (int i = 0; i < sources->n; i++) {
+	pmSource *source = sources->data[i];
+	if (!(source->mode & PM_SOURCE_MODE_PSFSTAR)) continue;
+	
+	// moments->Sx,Sy are roughly sigma_x,y
+	moments.x2 = PS_SQR(source->moments->Sx);
+	moments.y2 = PS_SQR(source->moments->Sy);
+	moments.xy = source->moments->Sxy;
+
+	axes = psEllipseMomentsToAxes (moments);
+
+	FWHM_X += axes.major * 2.35;
+	FWHM_Y += axes.minor * 2.35;
+	FWHM_T += axes.theta;
+	FWHM_N ++;
+    }
+
+    FWHM_X /= FWHM_N;
+    FWHM_Y /= FWHM_N;
+    FWHM_T /= FWHM_N;
+
+    psMetadataAdd (recipe, PS_LIST_TAIL, "FWHM_X", PS_DATA_F32 | PS_META_REPLACE, "PSF FWHM Major axis", FWHM_X);
+    psMetadataAdd (recipe, PS_LIST_TAIL, "FWHM_Y", PS_DATA_F32 | PS_META_REPLACE, "PSF FWHM Minor axis", FWHM_Y);
+    psMetadataAdd (recipe, PS_LIST_TAIL, "ANGLE",  PS_DATA_F32 | PS_META_REPLACE, "PSF angle",           FWHM_T);
+    
+    return true;
+}