Index: /branches/eam_branches/20090715/psphot/src/psphotReadoutKnownSources.c
===================================================================
--- /branches/eam_branches/20090715/psphot/src/psphotReadoutKnownSources.c	(revision 25527)
+++ /branches/eam_branches/20090715/psphot/src/psphotReadoutKnownSources.c	(revision 25528)
@@ -41,5 +41,5 @@
 
     // construct sources and measure basic stats
-    psArray *sources = psphotSourceStats (config, readout, detections);
+    psArray *sources = psphotSourceStats (config, readout, detections, true);
     if (!sources) return false;
 
Index: /branches/eam_branches/20090715/psphot/src/psphotReadoutMinimal.c
===================================================================
--- /branches/eam_branches/20090715/psphot/src/psphotReadoutMinimal.c	(revision 25527)
+++ /branches/eam_branches/20090715/psphot/src/psphotReadoutMinimal.c	(revision 25528)
@@ -54,5 +54,5 @@
 
     // construct sources and measure basic stats
-    psArray *sources = psphotSourceStats (config, readout, detections);
+    psArray *sources = psphotSourceStats (config, readout, detections, true);
     if (!sources) return false;
 
Index: /branches/eam_branches/20090715/psphot/src/psphotSourceStats.c
===================================================================
--- /branches/eam_branches/20090715/psphot/src/psphotSourceStats.c	(revision 25527)
+++ /branches/eam_branches/20090715/psphot/src/psphotSourceStats.c	(revision 25528)
@@ -1,5 +1,7 @@
 # include "psphotInternal.h"
 
-psArray *psphotSourceStats (pmConfig *config, pmReadout *readout, pmDetections *detections) {
+bool psphotSetMomentsWindow (psMetadata *recipe, psArray *sources);
+
+psArray *psphotSourceStats (pmConfig *config, pmReadout *readout, pmDetections *detections, bool setWindow) {
 
     bool status = false;
@@ -65,86 +67,9 @@
     }
 
-    // XXX *** this section is an attempt to iteratively determine the best value for sigma of the moments weighting Gaussian
-    {
-	float MIN_SN       = psMetadataLookupF32 (&status, recipe, "MOMENTS_SN_MIN");
-	if (!status) return false;
-
-	float sigma[4] = {1.0, 2.0, 3.0, 4.5};
-	float Sout[4];
-
-	// XXX check that this sorts by peak->SN
-	sources = psArraySort (sources, pmSourceSortBySN);
-
-	// loop over radii:
-	for (int i = 0; i < 4; i++) {
-
-	    for (int j = 0; (j < sources->n) && (j < 400); j++) {
- 
-		pmSource *source = sources->data[j];
-		psAssert (source->moments, "force moments to exist");
-		source->moments->nPixels = 0;
-
-		// skip faint sources for moments measurement
-		if (source->peak->SN < MIN_SN) {
-		    continue;
-		}
-
-		// measure basic source moments
-		status = pmSourceMoments (source, 4*sigma[i], sigma[i], 0.0);
-	    }
-
-	    // choose a grid scale that is a fixed fraction of the psf sigma^2
-	    psMetadataAddF32(recipe, PS_LIST_TAIL, "PSF_CLUMP_GRID_SCALE", PS_META_REPLACE, "clump grid", 0.1*PS_SQR(sigma[i]));
-	    psMetadataAddF32(recipe, PS_LIST_TAIL, "MOMENTS_SX_MAX", PS_META_REPLACE, "moments limit", 2.0*PS_SQR(sigma[i]));
-	    psMetadataAddF32(recipe, PS_LIST_TAIL, "MOMENTS_SY_MAX", PS_META_REPLACE, "moments limit", 2.0*PS_SQR(sigma[i]));
-
-	    // determine the PSF parameters from the source moment values
-	    pmPSFClump psfClump = pmSourcePSFClump (NULL, sources, recipe);
-	    psLogMsg ("psphot", 3, "radius %.1f, nStars: %d, nSigma: %5.2f, X,  Y: %f, %f (%f, %f)\n", sigma[i], psfClump.nStars, psfClump.nSigma, psfClump.X, psfClump.Y, sqrt(psfClump.X) / sigma[i], sqrt(psfClump.Y) / sigma[i]);
-
-	    psMetadataAddS32 (recipe, PS_LIST_TAIL, "PSF.CLUMP.NREGIONS",  PS_META_REPLACE, "psf clump regions", 1);
-	    psMetadata *regionMD = psMetadataLookupPtr (&status, recipe, "PSF.CLUMP.REGION.000");
-	    if (!regionMD) {
-		regionMD = psMetadataAlloc();
-		psMetadataAddMetadata (recipe, PS_LIST_TAIL, "PSF.CLUMP.REGION.000", PS_META_REPLACE, "psf clump region", regionMD);
-		psFree (regionMD);
-	    }
-	    psMetadataAddF32 (regionMD, PS_LIST_TAIL, "PSF.CLUMP.X",  PS_META_REPLACE, "psf clump center", psfClump.X);
-	    psMetadataAddF32 (regionMD, PS_LIST_TAIL, "PSF.CLUMP.Y",  PS_META_REPLACE, "psf clump center", psfClump.Y);
-	    psMetadataAddF32 (regionMD, PS_LIST_TAIL, "PSF.CLUMP.DX", PS_META_REPLACE, "psf clump center", psfClump.dX);
-	    psMetadataAddF32 (regionMD, PS_LIST_TAIL, "PSF.CLUMP.DY", PS_META_REPLACE, "psf clump center", psfClump.dY);
-	    
-	    // psphotVisualPlotMoments (recipe, sources);
-
-	    Sout[i] = sqrt(0.5*(psfClump.X + psfClump.Y)) / sigma[i];
+    if (setWindow) {
+	if (!psphotSetMomentsWindow(recipe, sources)) {
+	    psError(PS_ERR_UNEXPECTED_NULL, false, "Failed to determine Moments Window!");
+	    return NULL;
 	}
-
-	// we are looking for sigma for which Sout = 0.65 (or some other value)
-
-	float Sigma = NAN;
-	float minS = Sout[0];
-	float maxS = Sout[0];
-	for (int i = 0; i < 4; i++) {
-	    minS = PS_MIN(Sout[i], minS);
-	    maxS = PS_MAX(Sout[i], maxS);
-	}
-	if (minS > 0.65) Sigma = sigma[3];
-	if (maxS < 0.65) Sigma = sigma[0];
-
-	for (int i = 0; i < 3; i++) {
-	    if ((Sout[i] > 0.65) && (Sout[i+1] > 0.65)) continue;
-	    if ((Sout[i] < 0.65) && (Sout[i+1] < 0.65)) continue;
-	    Sigma = sigma[i] + (0.65 - Sout[i])*(sigma[i+1] - sigma[i])/(Sout[i+1] - Sout[i]);
-	}
-	psAssert (isfinite(Sigma), "did we miss a case?");
-	
-	// choose a grid scale that is a fixed fraction of the psf sigma^2
-	psMetadataAddF32(recipe, PS_LIST_TAIL, "PSF_CLUMP_GRID_SCALE", PS_META_REPLACE, "clump grid", 0.1*PS_SQR(Sigma));
-	psMetadataAddF32(recipe, PS_LIST_TAIL, "MOMENTS_SX_MAX", PS_META_REPLACE, "moments limit", 2.0*PS_SQR(Sigma));
-	psMetadataAddF32(recipe, PS_LIST_TAIL, "MOMENTS_SY_MAX", PS_META_REPLACE, "moments limit", 2.0*PS_SQR(Sigma));
-	psMetadataAddF32(recipe, PS_LIST_TAIL, "MOMENTS_GAUSS_SIGMA", PS_META_REPLACE, "moments limit", Sigma);
-	psMetadataAddF32(recipe, PS_LIST_TAIL, "PSF_MOMENTS_RADIUS", PS_META_REPLACE, "moments limit", 4.0*Sigma);
-
-	psLogMsg ("psphot", 3, "using window function with sigma = %f\n", Sigma);
     }
 
@@ -318,4 +243,94 @@
 }
 
+// this function attempts to iteratively determine the best value for sigma of the moments weighting Gaussian
+bool psphotSetMomentsWindow (psMetadata *recipe, psArray *sources) {
+
+    bool status;
+
+    float MIN_SN       = psMetadataLookupF32 (&status, recipe, "MOMENTS_SN_MIN");
+    if (!status) return false;
+
+    // XXX move this to a config file?
+    float sigma[4] = {1.0, 2.0, 3.0, 4.5};
+    float Sout[4];
+
+    // this sorts by peak->SN
+    sources = psArraySort (sources, pmSourceSortBySN);
+
+    // loop over radii:
+    for (int i = 0; i < 4; i++) {
+
+	// XXX move max source number to config
+	for (int j = 0; (j < sources->n) && (j < 400); j++) {
+ 
+	    pmSource *source = sources->data[j];
+	    psAssert (source->moments, "force moments to exist");
+	    source->moments->nPixels = 0;
+
+	    // skip faint sources for moments measurement
+	    if (source->peak->SN < MIN_SN) {
+		continue;
+	    }
+
+	    // measure basic source moments
+	    status = pmSourceMoments (source, 4*sigma[i], sigma[i], 0.0);
+	}
+
+	// choose a grid scale that is a fixed fraction of the psf sigma^2
+	psMetadataAddF32(recipe, PS_LIST_TAIL, "PSF_CLUMP_GRID_SCALE", PS_META_REPLACE, "clump grid", 0.1*PS_SQR(sigma[i]));
+	psMetadataAddF32(recipe, PS_LIST_TAIL, "MOMENTS_SX_MAX", PS_META_REPLACE, "moments limit", 2.0*PS_SQR(sigma[i]));
+	psMetadataAddF32(recipe, PS_LIST_TAIL, "MOMENTS_SY_MAX", PS_META_REPLACE, "moments limit", 2.0*PS_SQR(sigma[i]));
+
+	// determine the PSF parameters from the source moment values
+	pmPSFClump psfClump = pmSourcePSFClump (NULL, sources, recipe);
+	psLogMsg ("psphot", 3, "radius %.1f, nStars: %d, nSigma: %5.2f, X,  Y: %f, %f (%f, %f)\n", sigma[i], psfClump.nStars, psfClump.nSigma, psfClump.X, psfClump.Y, sqrt(psfClump.X) / sigma[i], sqrt(psfClump.Y) / sigma[i]);
+
+	psMetadataAddS32 (recipe, PS_LIST_TAIL, "PSF.CLUMP.NREGIONS",  PS_META_REPLACE, "psf clump regions", 1);
+	psMetadata *regionMD = psMetadataLookupPtr (&status, recipe, "PSF.CLUMP.REGION.000");
+	if (!regionMD) {
+	    regionMD = psMetadataAlloc();
+	    psMetadataAddMetadata (recipe, PS_LIST_TAIL, "PSF.CLUMP.REGION.000", PS_META_REPLACE, "psf clump region", regionMD);
+	    psFree (regionMD);
+	}
+	psMetadataAddF32 (regionMD, PS_LIST_TAIL, "PSF.CLUMP.X",  PS_META_REPLACE, "psf clump center", psfClump.X);
+	psMetadataAddF32 (regionMD, PS_LIST_TAIL, "PSF.CLUMP.Y",  PS_META_REPLACE, "psf clump center", psfClump.Y);
+	psMetadataAddF32 (regionMD, PS_LIST_TAIL, "PSF.CLUMP.DX", PS_META_REPLACE, "psf clump center", psfClump.dX);
+	psMetadataAddF32 (regionMD, PS_LIST_TAIL, "PSF.CLUMP.DY", PS_META_REPLACE, "psf clump center", psfClump.dY);
+	    
+	// psphotVisualPlotMoments (recipe, sources);
+
+	Sout[i] = sqrt(0.5*(psfClump.X + psfClump.Y)) / sigma[i];
+    }
+
+    // we are looking for sigma for which Sout = 0.65 (or some other value)
+
+    float Sigma = NAN;
+    float minS = Sout[0];
+    float maxS = Sout[0];
+    for (int i = 0; i < 4; i++) {
+	minS = PS_MIN(Sout[i], minS);
+	maxS = PS_MAX(Sout[i], maxS);
+    }
+    if (minS > 0.65) Sigma = sigma[3];
+    if (maxS < 0.65) Sigma = sigma[0];
+
+    for (int i = 0; i < 3; i++) {
+	if ((Sout[i] > 0.65) && (Sout[i+1] > 0.65)) continue;
+	if ((Sout[i] < 0.65) && (Sout[i+1] < 0.65)) continue;
+	Sigma = sigma[i] + (0.65 - Sout[i])*(sigma[i+1] - sigma[i])/(Sout[i+1] - Sout[i]);
+    }
+    psAssert (isfinite(Sigma), "did we miss a case?");
+	
+    // choose a grid scale that is a fixed fraction of the psf sigma^2
+    psMetadataAddF32(recipe, PS_LIST_TAIL, "PSF_CLUMP_GRID_SCALE", PS_META_REPLACE, "clump grid", 0.1*PS_SQR(Sigma));
+    psMetadataAddF32(recipe, PS_LIST_TAIL, "MOMENTS_SX_MAX", PS_META_REPLACE, "moments limit", 2.0*PS_SQR(Sigma));
+    psMetadataAddF32(recipe, PS_LIST_TAIL, "MOMENTS_SY_MAX", PS_META_REPLACE, "moments limit", 2.0*PS_SQR(Sigma));
+    psMetadataAddF32(recipe, PS_LIST_TAIL, "MOMENTS_GAUSS_SIGMA", PS_META_REPLACE, "moments limit", Sigma);
+    psMetadataAddF32(recipe, PS_LIST_TAIL, "PSF_MOMENTS_RADIUS", PS_META_REPLACE, "moments limit", 4.0*Sigma);
+
+    psLogMsg ("psphot", 3, "using window function with sigma = %f\n", Sigma);
+    return true;
+}
+
 # if (0)
 bool psphotSourceStats_Unthreaded (int *nfail, int *nmoments, psArray *sources, psMetadata *recipe) {