Index: trunk/psphot/src/psphotEnsemblePSF.c
===================================================================
--- trunk/psphot/src/psphotEnsemblePSF.c	(revision 6379)
+++ trunk/psphot/src/psphotEnsemblePSF.c	(revision 6427)
@@ -3,5 +3,5 @@
 // 2006.02.07 : no leaks!
 // fit all reasonable sources with the linear PSF model
-bool psphotEnsemblePSF (pmReadout *readout, psMetadata *config, psArray *sources, pmPSF *psf, psStats *sky) { 
+bool psphotEnsemblePSF (pmReadout *readout, psMetadata *config, psArray *sources, pmPSF *psf) { 
 
     bool  status;
@@ -15,11 +15,10 @@
     sources = psArraySort (sources, psphotSortByY);
 
+    // radius for source moments if star is saturated
+    // XXX we already do this in psphotSourceStats, is it still needed?
     float INNER_RADIUS     = psMetadataLookupF32 (&status, config, "SKY_INNER_RADIUS");
-    float PSF_FIT_NSIGMA   = psMetadataLookupF32 (&status, config, "PSF_FIT_NSIGMA");
-    float PSF_FIT_PADDING  = psMetadataLookupF32 (&status, config, "PSF_FIT_PADDING");
 
-    // set the object surface-brightness limit for fitted pixels
-    float FLUX_LIMIT  = PSF_FIT_NSIGMA * sky->sampleStdev;
-    psLogMsg ("psphot.ensemble", 4, "fitting pixels with at least %f object counts\n", FLUX_LIMIT);
+    // setup the PSF fit details
+    psphotInitRadiusPSF (config, psf->type);
 
     // pre-calculate all model pixels
@@ -29,4 +28,5 @@
     index->n = 0;
 
+    // option to limit analysis to a specific region
     bool UseAnalysisRegion = false;
     psRegion AnalysisRegion = {0, 0, 0, 0};
@@ -98,13 +98,6 @@
 	y = model->params->data.F32[3];
 
-	// get function which specifies the radius at which the model hits the given flux
-	pmModelRadius modelRadius = pmModelRadius_GetFunction (psf->type);
-
 	// set the fit radius based on the object flux limit and the model
-	model->radius = modelRadius (model->params, FLUX_LIMIT) + PSF_FIT_PADDING;
-	if (isnan(model->radius)) psAbort ("psphotEnsemblePSF", "error in radius");
-
-	// if needed, ask for more object pixels
-	psphotRedefinePixels (inSource, readout, x, y, model->radius);
+	psphotCheckRadiusPSF (readout, inSource, model);
 
 	// make temporary copies of the image pixels and mask
@@ -116,5 +109,5 @@
 	psImage *mask = otSource->mask;
 
-	// set model to unit peak, zero sky (we assume sky is constant)
+	// set model to unit peak, zero sky (we assume sky is subtracted)
 	model->params->data.F32[0] = 0.0;
 	model->params->data.F32[1] = 1.0;
@@ -165,5 +158,5 @@
     psLogMsg ("psphot.emsemble", 4, "built matrix: %f (%d elements)\n", psTimerMark ("psphot"), sparse->Nelem);
 
-    // solve for normalization terms (need include local sky)
+    // solve for normalization terms (need include local sky?)
     psSparseResort (sparse);
     psVector *norm = psSparseSolve (NULL, sparse, 3);