Index: trunk/psphot/src/psphotChoosePSF.c
===================================================================
--- trunk/psphot/src/psphotChoosePSF.c	(revision 23989)
+++ trunk/psphot/src/psphotChoosePSF.c	(revision 25755)
@@ -1,18 +1,3 @@
 # include "psphotInternal.h"
-
-void psphotCountPSFStars (psArray *sources) {
-
-    int nPSF = 0;
-
-    // select the candidate PSF stars (pointers to original sources)
-    for (int i = 0; i < sources->n; i++) {
-        pmSource *source = sources->data[i];
-        if (source->mode & PM_SOURCE_MODE_PSFSTAR) {
-            nPSF ++;
-        }
-    }
-
-    fprintf (stderr, "N PSF: %d\n", nPSF);
-}
 
 // try PSF models and select best option
@@ -73,6 +58,18 @@
     // get the fixed PSF fit radius
     // XXX check that PSF_FIT_RADIUS < SKY_OUTER_RADIUS
-    options->radius = psMetadataLookupF32 (&status, recipe, "PSF_FIT_RADIUS");
-    assert (status);
+    // options->radius = psMetadataLookupF32 (&status, recipe, "PSF_FIT_RADIUS");
+    // assert (status);
+
+    // We have calculated a Gaussian window function, use that for both the PSF fit radius and
+    // the aperture radius (scaling SIGMA)
+    float gaussSigma = psMetadataLookupF32(&status, recipe, "MOMENTS_GAUSS_SIGMA");
+    float fitScale = psMetadataLookupF32(&status, recipe, "PSF_FIT_RADIUS_SCALE");
+    float apScale = psMetadataLookupF32(&status, recipe, "PSF_APERTURE_SCALE");
+    options->fitRadius = (int)(fitScale*gaussSigma);
+    options->apRadius = (int)(apScale*gaussSigma);
+
+    // XXX use the same radii for standard analysis as for the PSF creation
+    psMetadataAddF32(recipe, PS_LIST_TAIL, "PSF_FIT_RADIUS", PS_META_REPLACE, "fit radius", options->fitRadius);
+    psMetadataAddF32(recipe, PS_LIST_TAIL, "PSF_APERTURE", PS_META_REPLACE, "psf aperture", options->apRadius);
 
     // XXX ROBUST seems to be too agressive given the small numbers.
@@ -99,6 +96,4 @@
 
     psArray *stars = psArrayAllocEmpty (sources->n);
-
-    // psphotCountPSFStars (sources);
 
     // select the candidate PSF stars (pointers to original sources)
@@ -115,11 +110,7 @@
     }
 
-    // psphotCountPSFStars (sources);
-
     // check that the identified psf stars sufficiently cover the region; if not, extend the
     // limits somewhat
     psphotCheckStarDistribution (stars, sources, options);
-
-    // psphotCountPSFStars (sources);
 
     psLogMsg ("psphot.pspsf", PS_LOG_DETAIL, "selected candidate %ld PSF objects\n", stars->n);
@@ -289,10 +280,8 @@
     // XXX test dump of psf star data and psf-subtracted image
     if (psTraceGetLevel("psphot.psfstars") > 5) {
-        psphotDumpPSFStars (readout, try, options->radius, maskVal, markVal);
+        psphotDumpPSFStars (readout, try, options->fitRadius, maskVal, markVal);
     }
 
     // save only the best model;
-    // XXX we are not saving the fitted sources
-    // XXX do we want to keep them so we may optionally write them out?
     pmPSF *psf = psMemIncrRefCounter(try->psf);
     psFree (models);
@@ -305,6 +294,4 @@
     }
 
-    // psphotCountPSFStars (sources);
-
     char *modelName = pmModelClassGetName (psf->type);
     psLogMsg ("psphot.pspsf", PS_LOG_INFO, "select psf model: %f sec\n", psTimerMark ("psphot.choose.psf"));
@@ -341,5 +328,8 @@
             // create modelPSF from this model
             pmModel *modelPSF = pmModelFromPSFforXY (psf, xc, yc, 1.0);
-            if (!modelPSF) continue;
+            if (!modelPSF) {
+		fprintf (stderr, "?");
+		continue;
+	    }
 
             // get the model full-width at half-max
@@ -354,5 +344,8 @@
 
             float FWHM_MINOR = FWHM_MAJOR * (axes.minor / axes.major);
-            if (!isfinite(FWHM_MAJOR) || !isfinite(FWHM_MINOR)) continue;
+            if (!isfinite(FWHM_MAJOR) || !isfinite(FWHM_MINOR)) {
+		fprintf (stderr, "!");
+		continue;
+	    }
             psVectorAppend (fwhmMajor, FWHM_MAJOR);
             psVectorAppend (fwhmMinor, FWHM_MINOR);