Index: branches/pap/psphot/src/Makefile.am
===================================================================
--- branches/pap/psphot/src/Makefile.am	(revision 25025)
+++ branches/pap/psphot/src/Makefile.am	(revision 25027)
@@ -25,5 +25,5 @@
 libpsphot_la_LDFLAGS = $(PSPHOT_LIBS) $(PSMODULE_LIBS) $(PSLIB_LIBS)
 
-bin_PROGRAMS = psphot psphotTest
+bin_PROGRAMS = psphot psphotTest psphotMomentsStudy
 
 psphot_CFLAGS = $(PSPHOT_CFLAGS) $(PSMODULE_CFLAGS) $(PSLIB_CFLAGS)
@@ -34,4 +34,8 @@
 psphotTest_LDFLAGS = $(PSPHOT_LIBS) $(PSMODULE_LIBS) $(PSLIB_LIBS)
 psphotTest_LDADD = libpsphot.la
+
+psphotMomentsStudy_CFLAGS = $(PSPHOT_CFLAGS) $(PSMODULE_CFLAGS) $(PSLIB_CFLAGS)
+psphotMomentsStudy_LDFLAGS = $(PSPHOT_LIBS) $(PSMODULE_LIBS) $(PSLIB_LIBS)
+psphotMomentsStudy_LDADD = libpsphot.la
 
 psphot_SOURCES = \
@@ -54,4 +58,7 @@
 psphotTest_SOURCES = \
         psphotTest.c
+
+psphotMomentsStudy_SOURCES = \
+        psphotMomentsStudy.c
 
 libpsphot_la_SOURCES = \
Index: branches/pap/psphot/src/psphot.c
===================================================================
--- branches/pap/psphot/src/psphot.c	(revision 25025)
+++ branches/pap/psphot/src/psphot.c	(revision 25027)
@@ -1,8 +1,3 @@
 # include "psphotStandAlone.h"
-
-static void usage (void) {
-    fprintf (stderr, "USAGE: psphot [-file image(s)] [-list imagelist] (output)\n");
-    exit (PS_EXIT_CONFIG_ERROR);
-}
 
 int main (int argc, char **argv) {
@@ -14,8 +9,5 @@
     // load command-line arguments, options, and system config data
     pmConfig *config = psphotArguments (argc, argv);
-    if (!config) {
-        psErrorStackPrint(stderr, "Error reading arguments\n");
-        usage ();
-    }
+    assert(config);
 
     psphotVersionPrint();
Index: branches/pap/psphot/src/psphot.h
===================================================================
--- branches/pap/psphot/src/psphot.h	(revision 25025)
+++ branches/pap/psphot/src/psphot.h	(revision 25027)
@@ -157,5 +157,5 @@
 bool            psphotSubWithTest (pmSource *source, bool useState, psImageMaskType maskVal);
 bool            psphotSetState (pmSource *source, bool curState, psImageMaskType maskVal);
-bool            psphotDeblendSatstars (psArray *sources, psMetadata *recipe);
+bool            psphotDeblendSatstars (pmReadout *readout, psArray *sources, psMetadata *recipe);
 bool            psphotSourceSize (pmConfig *config, pmReadout *readout, psArray *sources, psMetadata *recipe, long first);
 
Index: branches/pap/psphot/src/psphotApResid.c
===================================================================
--- branches/pap/psphot/src/psphotApResid.c	(revision 25025)
+++ branches/pap/psphot/src/psphotApResid.c	(revision 25027)
@@ -168,4 +168,7 @@
         if (source->mode &  PM_SOURCE_MODE_POOR) SKIPSTAR ("POOR STAR");
 
+	if (source->mode &  PM_SOURCE_MODE_EXT_LIMIT) SKIPSTAR ("EXTENDED");
+	if (source->mode &  PM_SOURCE_MODE_CR_LIMIT) SKIPSTAR ("COSMIC RAY");
+	    
         if (!isfinite(source->apMag) || !isfinite(source->psfMag)) {
             continue;
@@ -317,4 +320,5 @@
 */
 
+// XXX this still sucks...  need a better way to estimate the error floor...
 bool psphotMagErrorScale (float *errorScale, float *errorFloor, psVector *dMag, psVector *dap, psVector *mask, int nGroup) {
 
@@ -355,12 +359,15 @@
 
         if (j > 2) {
-            bool status = true;
-            status &= psVectorStats (statsS, dASubset, NULL, mkSubset, 0xff);
-            status &= psVectorStats (statsM, dMSubset, NULL, mkSubset, 0xff);
-            if (!status) { psErrorClear (); }
+            if (!psVectorStats (statsS, dASubset, NULL, mkSubset, 0xff)) {
+		psError(PS_ERR_UNKNOWN, false, "failure to measure stats");
+		return false;
+	    }
+            if (!psVectorStats (statsM, dMSubset, NULL, mkSubset, 0xff)) {
+		psError(PS_ERR_UNKNOWN, false, "failure to measure stats");
+		return false;
+	    }
             dSo->data.F32[i] = statsS->robustStdev;
             dMo->data.F32[i] = statsM->sampleMean;
             dRo->data.F32[i] = statsS->robustStdev / statsM->sampleMean;
-            //      fprintf (stderr, "%d (%d) : sys: %f, phot: %f, rat: %f\n", i, j, dSo->data.F32[i], dMo->data.F32[i], dRo->data.F32[i]);
         } else {
             dSo->data.F32[i] = NAN;
@@ -375,6 +382,6 @@
     psStats *stats = psStatsAlloc (PS_STAT_SAMPLE_MEDIAN);
     if (!psVectorStats (stats, dRo, NULL, NULL, 0)) {
-        // XXX better testing of raised error
-        psErrorClear();
+	psError(PS_ERR_UNKNOWN, false, "failure to measure stats");
+	return false;
     }
 
@@ -382,4 +389,5 @@
     for (int i = 0; i < dSo->n; i++) {
         *errorFloor = dSo->data.F32[i];
+        if (fabs(*errorFloor) <= FLT_EPSILON) continue;
         if (isfinite(*errorFloor)) break;
     }
Index: branches/pap/psphot/src/psphotArguments.c
===================================================================
--- branches/pap/psphot/src/psphotArguments.c	(revision 25025)
+++ branches/pap/psphot/src/psphotArguments.c	(revision 25027)
@@ -1,3 +1,99 @@
 # include "psphotStandAlone.h"
+
+static void writeHelpInfo(const char* program, pmConfig* config, FILE* ofile)
+{
+  fprintf(ofile,
+	  "Usage: one of the following\n"
+	  "%s -file fname1[,fname2,...] -mask maskfile1[,maskfile2,...]\n"
+	  "     -variance varfile1[,varfile2,...] OutFileBaseName\n"
+	  "\n"
+	  "%s -list FileNameList [-masklist MaskFileNameList] \n"
+	  "     -variancelist VarFileNameList OutFileBaseName\n"
+	  "\n"
+	  "%s -help\n"
+	  "\n"
+	  "%s -version\n"
+	  "\n"
+	  "where:\n"
+	  "  FileNameList is a text file containing filenames, one per line\n"
+	  "  MaskFileNameList is a text file of mask filenames, one per line\n"
+	  "  VarFileNameList is a text file of variance filenames, one per line\n"
+	  "  OutFileBaseName is the 'root name' for output files\n"
+	  "\n"
+	  "additional options:\n"
+	  "  -modeltest xObj yObj [-model DEFAULT|ModelName] \n"
+	  "        [-fitmode DEFAULT|PSF|CONV] [-fitset FitFileName]\n"
+	  "     Test fit for object at the given coordinates.  ModelName\n"
+	  "     is one of PS_MODEL_GAUSS, PS_MODEL_PGAUSS, PS_MODEL_QGAUSS,\n"
+	  "     PS_MODEL_RGAUSS, PS_MODEL_PS1_V1, or PS_MODEL_SERSIC.\n"
+	  "     FitFileName is a file of x,y,Io triples.\n"
+	  "  -psf PsfFile1[,PsfFile2,...] or -psflist PsfFileNameList\n"
+	  "     specify PSF rather than letting %s estimate it\n"
+	  "  -src SrcFile1[,SrcFile2,...] or -srclist SrcFileNameList\n"
+	  "     specify additional sources for PSF generation\n"
+	  "  -chip nn[,nn,...]\n"
+	  "     select detector chips to process; default is all.\n"
+	  "     Indices correspond to zero-based offset in the FPA metadata table.\n"
+	  "  -photcode PhotoCodeName\n"
+	  "     specify photocode\n"
+	  "  -region RegionString\n"
+	  "     specify analysis region.  String is of form '[x0:x1,y0:y1]'\n"
+	  "     To use this option you must define a default in psphot.config\n"
+	  "  -visual\n"
+	  "     turns on interactive display mode\n"
+	  "  -dumpconfig CfgFileName\n"
+          "     causes config info to be dumped to the named file.\n"
+	  "  -break NOTHING|BACKMDL|PEAKS|MOMENTS|PSFMODEL|ENSEMBLE|PASS1\n"
+	  "     choose a point at which to exit processing early\n"
+	  "  -nthreads n\n"
+	  "     set number of parallel threads of execution\n"
+	  "  -F OldFileRule ReplacementFileRule\n"
+	  "     change file naming rule; e.g. '-F PSPHOT.OUTPUT PSPHOT.OUT.CMF.MEF'\n"
+	  "  -D name stringval\n"
+	  "     set a string-valued config parameter\n"
+	  "  -Di name intval\n"
+	  "     set an integer-valued config parameter\n"
+	  "  -Df name fval\n"
+	  "     set a float-valued config parameter\n"
+	  "  -Db name boolval\n"
+	  "     set a boolean-valued config parameter\n"
+	  "  -v, -vv, -vvv\n"
+	  "     set increasing levels of verbosity\n"
+	  "  -logfmt FormatString\n"
+	  "     set format string used for log messages\n"
+	  "  -trace Fac Lvl\n"
+	  "     set tracing for facility Fac to integer Lvl, e.g. '-trace err 10'\n"
+	  "  -trace-levels\n"
+	  "     print current trace levels\n",
+	  program,program,program,program,program);
+    psFree(config);
+    pmConfigDone();
+    psLibFinalize();
+    exit(PS_EXIT_SUCCESS);
+}
+
+static void usage(const char *program,  // Name of the program
+                  psMetadata *arguments, // Command-line arguments
+                  pmConfig *config,      // Configuration
+		  int exitCode
+		  ) 
+{
+  fprintf(stderr,
+	  "Usage: one of the following\n"
+	  "%s -file fname1[,fname2,...] -mask maskfile1[,maskfile2,...]\n"
+	  "     -variance varfile1[,varfile2,...] OutFileBaseName\n"
+	  "\n"
+	  "%s -list FileNameList [-masklist MaskFileNameList] \n"
+	  "     -variancelist VarFileNameList OutFileBaseName\n"
+	  "\n"
+	  "Try '%s -help' for more options and explanation\n",
+	  program,program,program);
+    if (exitCode != PS_EXIT_SUCCESS)
+      psErrorStackPrint(stderr, "Error reading arguments\n");
+    psFree(config);
+    pmConfigDone();
+    psLibFinalize();
+    exit(exitCode);
+}
 
 pmConfig *psphotArguments(int argc, char **argv) {
@@ -5,24 +101,13 @@
     int N;
     bool status;
-
-    if (argc == 1) {
-        psError(PSPHOT_ERR_ARGUMENTS, true, "Too few arguments: %d", argc);
-        return NULL;
-    }
-
-    if ((N = psArgumentGet (argc, argv, "-version"))) {
-        psString version;
-        version = psphotVersionLong();    fprintf (stdout, "%s\n", version); psFree (version);
-        version = psModulesVersionLong(); fprintf (stdout, "%s\n", version); psFree (version);
-        version = psLibVersionLong();     fprintf (stdout, "%s\n", version); psFree (version);
-        exit (0);
-    }
 
     // load config data from default locations
     pmConfig *config = pmConfigRead(&argc, argv, PSPHOT_RECIPE);
     if (config == NULL) {
-        psError(PSPHOT_ERR_CONFIG, false, "Can't read site configuration");
-        return NULL;
-    }
+      psErrorStackPrint(stderr, "Can't read site configuration");
+	exit(PS_EXIT_CONFIG_ERROR);
+    }
+
+    PSARGUMENTS_INSTANTIATE_GENERICS( psphot, config, argc, argv );
 
     // save the following additional recipe values based on command-line options
@@ -30,18 +115,14 @@
     psMetadata *options = pmConfigRecipeOptions (config, PSPHOT_RECIPE);
 
-    // Number of threads
-    if ((N = psArgumentGet(argc, argv, "-threads"))) {
-        psArgumentRemove(N, &argc, argv);
-        int nThreads = atoi(argv[N]);
-        psMetadataAddS32(config->arguments, PS_LIST_TAIL, "NTHREADS", 0, "number of psphot threads", nThreads);
-        psArgumentRemove(N, &argc, argv);
-
-        // create the thread pool with number of desired threads, supplying our thread launcher function
-        psThreadPoolInit (nThreads);
-    }
-    psphotSetThreads();
+    // Number of threads is handled
+    PSARGUMENTS_INSTANTIATE_THREADSARG( psphot, config, argc, argv )
 
     // run the test model (requires X,Y coordinate)
     if ((N = psArgumentGet (argc, argv, "-modeltest"))) {
+        if (argc<=N+2) {
+          psErrorStackPrint(stderr, "Expected to see 2 more arguments; saw %d", argc - 1);
+	  exit(PS_EXIT_CONFIG_ERROR);
+	}
+
         psMetadataAddBool (options, PS_LIST_TAIL, "TEST_FIT",   0, "", true);
         psMetadataAddF32  (options, PS_LIST_TAIL, "TEST_FIT_X", 0, "", atof(argv[N+1]));
@@ -54,4 +135,8 @@
         // specify the modeltest model
         if ((N = psArgumentGet (argc, argv, "-model"))) {
+	    if (argc<=N+1) {
+	      psErrorStackPrint(stderr, "Expected to see 1 more argument; saw %d", argc - 1);
+	      exit(PS_EXIT_CONFIG_ERROR);
+	    }
             psArgumentRemove (N, &argc, argv);
             psMetadataAddStr (options, PS_LIST_TAIL, "TEST_FIT_MODEL", 0, "", argv[N]);
@@ -61,4 +146,8 @@
         // specify the test fit mode
         if ((N = psArgumentGet (argc, argv, "-fitmode"))) {
+	    if (argc<=N+1) {
+	      psErrorStackPrint(stderr, "Expected to see 1 more argument; saw %d", argc - 1);
+	      exit(PS_EXIT_CONFIG_ERROR);
+	    }
             psArgumentRemove (N, &argc, argv);
             psMetadataAddStr (options, PS_LIST_TAIL, "TEST_FIT_MODE", 0, "", argv[N]);
@@ -66,4 +155,8 @@
         }
         if ((N = psArgumentGet (argc, argv, "-fitset"))) {
+	    if (argc<=N+1) {
+	      psErrorStackPrint(stderr, "Expected to see 1 more argument; saw %d", argc - 1);
+	      exit(PS_EXIT_CONFIG_ERROR);
+	    }
             psArgumentRemove (N, &argc, argv);
             psMetadataAddStr (options, PS_LIST_TAIL, "TEST_FIT_SET", 0, "", argv[N]);
@@ -74,4 +167,8 @@
     // photcode : used in output to supplement header data (argument or recipe?)
     if ((N = psArgumentGet (argc, argv, "-photcode"))) {
+        if (argc<=N+1) {
+	  psErrorStackPrint(stderr, "Expected to see 1 more argument; saw %d", argc - 1);
+	  exit(PS_EXIT_CONFIG_ERROR);
+	}
         psArgumentRemove (N, &argc, argv);
         psMetadataAddStr (options, PS_LIST_TAIL, "PHOTCODE", PS_META_REPLACE, "", argv[N]);
@@ -87,13 +184,10 @@
     // break : used from recipe throughout psphotReadout
     if ((N = psArgumentGet (argc, argv, "-break"))) {
+	if (argc<=N+1) {
+	  psErrorStackPrint(stderr, "Expected to see 1 more argument; saw %d", argc - 1);
+	  exit(PS_EXIT_CONFIG_ERROR);
+	}
         psArgumentRemove (N, &argc, argv);
         psMetadataAddStr (options, PS_LIST_TAIL, "BREAK_POINT", PS_META_REPLACE, "", argv[N]);
-        psArgumentRemove (N, &argc, argv);
-    }
-
-    // fitmode : used from recipe throughout psphotReadout
-    if ((N = psArgumentGet (argc, argv, "-fitmode"))) {
-        psArgumentRemove (N, &argc, argv);
-        psMetadataAddStr (options, PS_LIST_TAIL, "FITMODE", PS_META_REPLACE, "", argv[N]);
         psArgumentRemove (N, &argc, argv);
     }
@@ -101,4 +195,8 @@
     // analysis region : overrides recipe value, used in psphotReadout/psphotEnsemblePSF
     if ((N = psArgumentGet (argc, argv, "-region"))) {
+	if (argc<=N+1) {
+	  psErrorStackPrint(stderr, "Expected to see 1 more argument; saw %d", argc - 1);
+	  exit(PS_EXIT_CONFIG_ERROR);
+	}
         psArgumentRemove (N, &argc, argv);
         psMetadataAddStr (options, PS_LIST_TAIL, "ANALYSIS_REGION", 0, "", argv[N]);
@@ -108,4 +206,8 @@
     // chip selection is used to limit chips to be processed
     if ((N = psArgumentGet (argc, argv, "-chip"))) {
+	if (argc<=N+1) {
+	  psErrorStackPrint(stderr, "Expected to see 1 more argument; saw %d", argc - 1);
+	  exit(PS_EXIT_CONFIG_ERROR);
+	}
         psArgumentRemove (N, &argc, argv);
         psMetadataAddStr (config->arguments, PS_LIST_TAIL, "CHIP_SELECTIONS", PS_DATA_STRING, "", argv[N]);
@@ -121,14 +223,23 @@
     pmConfigFileSetsMD (config->arguments, &argc, argv, "SRC",    "-src",    "-srclist");
 
+    if (argc == 1) {
+        psError(PSPHOT_ERR_ARGUMENTS, true, "Too few arguments: %d", argc);
+	usage(argv[0], config->arguments, config, PS_EXIT_CONFIG_ERROR);
+    }
+
+    if (psArgumentGet(argc, argv, "-help")
+	|| psArgumentGet(argc, argv, "-h"))
+      writeHelpInfo(argv[0], config, stdout);
+      
     // the input file is a required argument; if not found, we will exit
     status = pmConfigFileSetsMD (config->arguments, &argc, argv, "INPUT", "-file", "-list");
     if (!status) {
         psError(PSPHOT_ERR_ARGUMENTS, false, "pmConfigFileSetsMD failed to parse arguments");
-        return NULL;
+	usage(argv[0], config->arguments, config, PS_EXIT_CONFIG_ERROR);
     }
 
     if (argc != 2) {
         psError(PSPHOT_ERR_ARGUMENTS, true, "Expected to see one more argument; saw %d", argc - 1);
-        return NULL;
+	usage(argv[0], config->arguments, config, PS_EXIT_CONFIG_ERROR);
     }
 
Index: branches/pap/psphot/src/psphotBlendFit.c
===================================================================
--- branches/pap/psphot/src/psphotBlendFit.c	(revision 25025)
+++ branches/pap/psphot/src/psphotBlendFit.c	(revision 25027)
@@ -197,5 +197,5 @@
         if (source->type == PM_SOURCE_TYPE_SATURATED) continue;
 
-        // skip DBL second sources (ie, added by psphotFitBlob)
+        // skip DBL second sources (ie, added by psphotFitBlob
         if (source->mode &  PM_SOURCE_MODE_PAIR) continue;
 
Index: branches/pap/psphot/src/psphotChoosePSF.c
===================================================================
--- branches/pap/psphot/src/psphotChoosePSF.c	(revision 25025)
+++ branches/pap/psphot/src/psphotChoosePSF.c	(revision 25027)
@@ -361,5 +361,9 @@
 
     psStats *stats = psStatsAlloc (PS_STAT_SAMPLE_MEAN | PS_STAT_SAMPLE_STDEV | PS_STAT_SAMPLE_QUARTILE);
-    psVectorStats (stats, fwhmMajor, NULL, NULL, 0);
+    if (!psVectorStats (stats, fwhmMajor, NULL, NULL, 0)) {
+        psError(PS_ERR_UNKNOWN, false, "failure to measure stats for FWHM MAJOR");
+        return false;
+    }
+
     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);
@@ -367,5 +371,8 @@
     psMetadataAddF32 (recipe, PS_LIST_TAIL, "FW_MJ_UQ",   PS_META_REPLACE, "PSF FWHM Major axis (upper quartile)", stats->sampleUQ);
 
-    psVectorStats (stats, fwhmMinor, NULL, NULL, 0);
+    if (!psVectorStats (stats, fwhmMinor, NULL, NULL, 0)) {
+        psError(PS_ERR_UNKNOWN, false, "failure to measure stats for FWHM MINOR");
+        return false;
+    }
     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);
Index: branches/pap/psphot/src/psphotCleanup.c
===================================================================
--- branches/pap/psphot/src/psphotCleanup.c	(revision 25025)
+++ branches/pap/psphot/src/psphotCleanup.c	(revision 25027)
@@ -2,4 +2,13 @@
 
 void psphotCleanup (pmConfig *config) {
+
+    // Dump configuration if requested
+    bool status;
+    psString dump_file = psMetadataLookupStr(&status, config->arguments, "DUMP_CONFIG");
+    if (dump_file) {
+        (void)pmConfigCamerasCull(config,NULL);
+	(void)pmConfigRecipesCull(config,NULL);
+        pmConfigDump(config, dump_file);
+    }
 
     psFree (config);
Index: branches/pap/psphot/src/psphotCullPeaks.c
===================================================================
--- branches/pap/psphot/src/psphotCullPeaks.c	(revision 25025)
+++ branches/pap/psphot/src/psphotCullPeaks.c	(revision 25027)
@@ -18,25 +18,14 @@
         nsigma_min = 0;
     }
+    float fPadding = psMetadataLookupF32(&status, recipe, "FOOTPRINT_CULL_NSIGMA_PAD");
+    if (!status) {
+        fPadding = 0;
+    }
     const float skyStdev = psMetadataLookupF32(NULL, recipe, "SKY_STDEV");
-
-    return pmFootprintArrayCullPeaks(image, weight, footprints,
-                                     nsigma_delta, nsigma_min*skyStdev);
-}
-
-
-/*
- * Cull an entire psArray of pmFootprints
- * XXX drop this intermediate level function?
- */
-psErrorCode
-pmFootprintArrayCullPeaks(const psImage *img, // the image wherein lives the footprint
-			  const psImage *weight,	// corresponding variance image
-			  psArray *footprints, // array of pmFootprints
-			  const float nsigma_delta, // how many sigma above local background a peak
-    					// needs to be to survive
-			  const float min_threshold) { // minimum permitted coll height
+    const float min_threshold = nsigma_min*skyStdev;
+    
     for (int i = 0; i < footprints->n; i++) {
 	pmFootprint *fp = footprints->data[i];
-	if (pmFootprintCullPeaks(img, weight, fp, nsigma_delta, min_threshold) != PS_ERR_NONE) {
+	if (pmFootprintCullPeaks(image, weight, fp, nsigma_delta, fPadding, min_threshold) != PS_ERR_NONE) {
 	    return psError(PS_ERR_UNKNOWN, false, "Culling pmFootprint %d", fp->id);
 	}
@@ -45,220 +34,2 @@
     return PS_ERR_NONE;
 }
-
- /*
-  * Examine the peaks in a pmFootprint, and throw away the ones that are not sufficiently
-  * isolated.  More precisely, for each peak find the highest coll that you'd have to traverse
-  * to reach a still higher peak --- and if that coll's more than nsigma DN below your
-  * starting point, discard the peak.
-  */
-psErrorCode pmFootprintCullPeaks_OLD(const psImage *img, // the image wherein lives the footprint
-				 const psImage *weight,	// corresponding variance image
-				 pmFootprint *fp, // Footprint containing mortal peaks
-				 const float nsigma_delta, // how many sigma above local background a peak
-				 	// needs to be to survive
-				 const float min_threshold) { // minimum permitted coll height
-    assert (img != NULL); assert (img->type.type == PS_TYPE_F32);
-    assert (weight != NULL); assert (weight->type.type == PS_TYPE_F32);
-    assert (img->row0 == weight->row0 && img->col0 == weight->col0);
-    assert (fp != NULL);
-
-    if (fp->peaks == NULL || fp->peaks->n == 0) { // nothing to do
-	return PS_ERR_NONE;
-    }
-
-    psRegion subRegion;			// desired subregion; 1 larger than bounding box (grr)
-    subRegion.x0 = fp->bbox.x0; subRegion.x1 = fp->bbox.x1 + 1;
-    subRegion.y0 = fp->bbox.y0; subRegion.y1 = fp->bbox.y1 + 1;
-    const psImage *subImg = psImageSubset((psImage *)img, subRegion);
-    const psImage *subWt = psImageSubset((psImage *)weight, subRegion);
-    assert (subImg != NULL && subWt != NULL);
-    //
-    // We need a psArray of peaks brighter than the current peak.  We'll fake this
-    // by reusing the fp->peaks but lying about n.
-    //
-    // We do this for efficiency (otherwise I'd need two peaks lists), and we are
-    // rather too chummy with psArray in consequence.  But it works.
-    //
-    psArray *brightPeaks = psArrayAlloc(0);
-    psFree(brightPeaks->data);
-    brightPeaks->data = psMemIncrRefCounter(fp->peaks->data);// use the data from fp->peaks
-    //
-    // The brightest peak is always safe; go through other peaks trying to cull them
-    //
-    for (int i = 1; i < fp->peaks->n; i++) { // n.b. fp->peaks->n can change within the loop
-	const pmPeak *peak = fp->peaks->data[i];
-	int x = peak->x - subImg->col0;
-	int y = peak->y - subImg->row0;
-	//
-	// Find the level nsigma below the peak that must separate the peak
-	// from any of its friends
-	//
-	assert (x >= 0 && x < subImg->numCols && y >= 0 && y < subImg->numRows);
-	const float stdev = sqrt(subWt->data.F32[y][x]);
-	float threshold = subImg->data.F32[y][x] - nsigma_delta*stdev;
-	if (isnan(threshold) || threshold < min_threshold) {
-#if 1	  // min_threshold is assumed to be below the detection threshold,
-	  // so all the peaks are pmFootprint, and this isn't the brightest
-	    // XXX mark peak to be dropped
-	    (void)psArrayRemoveIndex(fp->peaks, i);
-	    i--;			// we moved everything down one
-	    continue;
-#else
-#error n.b. We will be running LOTS of checks at this threshold, so only find the footprint once
-	    threshold = min_threshold;
-#endif
-	}
-
-	// XXX EAM : if stdev >= 0, i'm not sure how this can ever be true?
-	if (threshold > subImg->data.F32[y][x]) {
-	    threshold = subImg->data.F32[y][x] - 10*FLT_EPSILON;
-	}
-
-	// XXX this is a bit expensive: psImageAlloc for every peak contained in this footprint
-	// perhaps this should alloc a single ID image above and pass it in to be set.
-
-	const int peak_id = 1;		// the ID for the peak of interest
-	brightPeaks->n = i;		// only stop at a peak brighter than we are
-
-	// XXX optionally use the faster pmFootprintsFind if the subimage size is large (eg, M31)
-
-	pmFootprint *peakFootprint = pmFootprintsFindAtPoint(subImg, threshold, brightPeaks, peak->y, peak->x);
-	brightPeaks->n = 0;		// don't double free
-	psImage *idImg = pmSetFootprintID(NULL, peakFootprint, peak_id);
-	psFree(peakFootprint);
-
-	// Check if any of the previous (brighter) peaks are within the footprint of this peak
-	// If so, the current peak is bogus; drop it.
-	int j;
-	for (j = 0; j < i; j++) {
-	    const pmPeak *peak2 = fp->peaks->data[j];
-	    int x2 = peak2->x - subImg->col0;
-	    int y2 = peak2->y - subImg->row0;
-	    const int peak2_id = idImg->data.S32[y2][x2]; // the ID for some other peak
-
-	    if (peak2_id == peak_id) {	// There's a brighter peak within the footprint above
-		;			// threshold; so cull our initial peak
-		(void)psArrayRemoveIndex(fp->peaks, i);
-		i--;			// we moved everything down one
-		break;
-	    }
-	}
-	if (j == i) {
-	    j++;
-	}
-
-	psFree(idImg);
-    }
-
-    brightPeaks->n = 0; psFree(brightPeaks);
-    psFree((psImage *)subImg);
-    psFree((psImage *)subWt);
-
-    return PS_ERR_NONE;
-}
-
- /*
-  * Examine the peaks in a pmFootprint, and throw away the ones that are not sufficiently
-  * isolated.  More precisely, for each peak find the highest coll that you'd have to traverse
-  * to reach a still higher peak --- and if that coll's more than nsigma DN below your
-  * starting point, discard the peak.
-  */
-
-# define IN_PEAK 1 
-psErrorCode pmFootprintCullPeaks(const psImage *img, // the image wherein lives the footprint
-				 const psImage *weight,	// corresponding variance image
-				 pmFootprint *fp, // Footprint containing mortal peaks
-				 const float nsigma_delta, // how many sigma above local background a peak
-				 // needs to be to survive
-				 const float min_threshold) { // minimum permitted coll height
-    assert (img != NULL); assert (img->type.type == PS_TYPE_F32);
-    assert (weight != NULL); assert (weight->type.type == PS_TYPE_F32);
-    assert (img->row0 == weight->row0 && img->col0 == weight->col0);
-    assert (fp != NULL);
-
-    if (fp->peaks == NULL || fp->peaks->n == 0) { // nothing to do
-	return PS_ERR_NONE;
-    }
-
-    psRegion subRegion;			// desired subregion; 1 larger than bounding box (grr)
-    subRegion.x0 = fp->bbox.x0; subRegion.x1 = fp->bbox.x1 + 1;
-    subRegion.y0 = fp->bbox.y0; subRegion.y1 = fp->bbox.y1 + 1;
-
-    psImage *subImg = psImageSubset((psImage *)img, subRegion);
-    psImage *subWt = psImageSubset((psImage *)weight, subRegion);
-    assert (subImg != NULL && subWt != NULL);
-
-    psImage *idImg = psImageAlloc(subImg->numCols, subImg->numRows, PS_TYPE_S32);
-
-    // We need a psArray of peaks brighter than the current peak.  
-    // We reject peaks which either:
-    // 1) are below the local threshold
-    // 2) have a brighter peak within their threshold
-
-    // allocate the full-sized array.  if the final array is much smaller, we can realloc
-    // at that point.
-    psArray *brightPeaks = psArrayAllocEmpty(fp->peaks->n);
-    psArrayAdd (brightPeaks, 128, fp->peaks->data[0]);
-
-    // The brightest peak is always safe; go through other peaks trying to cull them
-    for (int i = 1; i < fp->peaks->n; i++) { // n.b. fp->peaks->n can change within the loop
-	const pmPeak *peak = fp->peaks->data[i];
-	int x = peak->x - subImg->col0;
-	int y = peak->y - subImg->row0;
-	//
-	// Find the level nsigma below the peak that must separate the peak
-	// from any of its friends
-	//
-	assert (x >= 0 && x < subImg->numCols && y >= 0 && y < subImg->numRows);
-	const float stdev = sqrt(subWt->data.F32[y][x]);
-	float threshold = subImg->data.F32[y][x] - nsigma_delta*stdev;
-	if (isnan(threshold) || threshold < min_threshold) {
-	    // min_threshold is assumed to be below the detection threshold,
-	    // so all the peaks are pmFootprint, and this isn't the brightest
-	    continue;
-	}
-
-	// XXX EAM : if stdev >= 0, i'm not sure how this can ever be true?
-	if (threshold > subImg->data.F32[y][x]) {
-	    threshold = subImg->data.F32[y][x] - 10*FLT_EPSILON;
-	}
-
-	// XXX this is a bit expensive: psImageAlloc for every peak contained in this footprint
-	// perhaps this should alloc a single ID image above and pass it in to be set.
-
-	// XXX optionally use the faster pmFootprintsFind if the subimage size is large (eg, M31)
-
-	// at this point brightPeaks only has the peaks brighter than the current
-	pmFootprint *peakFootprint = pmFootprintsFindAtPoint(subImg, threshold, brightPeaks, peak->y, peak->x);
-
-	// XXX need to supply the image here
-	// we set the IDs to either 1 (in peak) or 0 (not in peak)
-	pmSetFootprintID (idImg, peakFootprint, IN_PEAK);
-	psFree(peakFootprint);
-
-	// Check if any of the previous (brighter) peaks are within the footprint of this peak
-	// If so, the current peak is bogus; drop it.
-	bool keep = true;
-	for (int j = 0; keep && (j < brightPeaks->n); j++) {
-	    const pmPeak *peak2 = fp->peaks->data[j];
-	    int x2 = peak2->x - subImg->col0;
-	    int y2 = peak2->y - subImg->row0;
-	    if (idImg->data.S32[y2][x2] == IN_PEAK) 
-		// There's a brighter peak within the footprint above threshold; so cull our initial peak
-		keep = false;
-	}
-	if (!keep) continue;
-
-	psArrayAdd (brightPeaks, 128, fp->peaks->data[i]);
-    }
-
-    psFree (fp->peaks);
-    fp->peaks = brightPeaks;
-
-    psFree(idImg);
-    psFree(subImg);
-    psFree(subWt);
-
-    return PS_ERR_NONE;
-}
-
Index: branches/pap/psphot/src/psphotDeblendSatstars.c
===================================================================
--- branches/pap/psphot/src/psphotDeblendSatstars.c	(revision 25025)
+++ branches/pap/psphot/src/psphotDeblendSatstars.c	(revision 25027)
@@ -1,5 +1,5 @@
 # include "psphotInternal.h"
 
-bool psphotDeblendSatstars (psArray *sources, psMetadata *recipe) {
+bool psphotDeblendSatstars (pmReadout *readout, psArray *sources, psMetadata *recipe) {
 
     int N;
@@ -9,6 +9,10 @@
 
     int Nblend = 0;
-    float SAT_TEST_LEVEL = 50000;
     float SAT_MIN_RADIUS = 5.0;
+
+    bool status;
+    pmCell *cell = readout->parent;
+    float SATURATION = 0.75*psMetadataLookupF32 (&status, cell->concepts, "CELL.SATURATION");
+    float SAT_TEST_LEVEL = 0.5*SATURATION;
 
     // we need sources spatially-sorted to find overlaps
@@ -32,5 +36,5 @@
         // XXX filter? if (source->mode & PM_SOURCE_MODE_SATSTAR) continue;
         if (source->mode & PM_SOURCE_MODE_BLEND) continue;
-        if (source->peak->flux < SAT_TEST_LEVEL) continue;
+        if (source->peak->flux < SATURATION) continue;
 
 	// save these for reference below
@@ -45,6 +49,5 @@
 	psVector *xVec = contour->data[0];
 	psVector *yVec = contour->data[1];
-
-	// XXX should we filter based on the number of pixels in the contour?
+	if (xVec->n < 5) continue;
 
 	// find the center of the contour (let's just use mid[x,y])
@@ -62,4 +65,9 @@
 	int yCenter = 0.5*(yMin + yMax);
 	psFree (contour);
+
+	psAssert (xCenter >= source->pixels->col0, "invalid shift in object center");
+	psAssert (xCenter <  source->pixels->col0 + source->pixels->numCols, "invalid shift in object center");
+	psAssert (yCenter >= source->pixels->row0, "invalid shift in object center");
+	psAssert (yCenter <  source->pixels->row0 + source->pixels->numRows, "invalid shift in object center");
 
 	// reset the peak for this source to the value of the center pixel
Index: branches/pap/psphot/src/psphotDiagnosticPlots.c
===================================================================
--- branches/pap/psphot/src/psphotDiagnosticPlots.c	(revision 25025)
+++ branches/pap/psphot/src/psphotDiagnosticPlots.c	(revision 25027)
@@ -40,5 +40,8 @@
 
     psStats *stats = psStatsAlloc (PS_STAT_MAX | PS_STAT_MIN);
-    psVectorStats (stats, values, NULL, NULL, 0);
+    if (!psVectorStats (stats, values, NULL, NULL, 0)) {
+	psError(PS_ERR_UNKNOWN, false, "failure to measure stats for histogram");
+        return false;
+    }
 
     psHistogram *histogram = psHistogramAlloc (stats->min, stats->max, 1000);
@@ -63,5 +66,8 @@
     psStatsInit (stats);
     stats->options = PS_STAT_MAX | PS_STAT_MIN;
-    psVectorStats (stats, histogram->nums, NULL, NULL, 0);
+    if (!psVectorStats (stats, histogram->nums, NULL, NULL, 0)) {
+	psError(PS_ERR_UNKNOWN, false, "failure to measure stats for histogram");
+        return false;
+    }
 
     // scale the plot to hold the histogram
@@ -100,5 +106,8 @@
     psStatsInit (stats);
     stats->options = PS_STAT_MAX | PS_STAT_MIN;
-    psVectorStats (stats, histogram->nums, NULL, NULL, 0);
+    if (!psVectorStats (stats, histogram->nums, NULL, NULL, 0)) {
+	psError(PS_ERR_UNKNOWN, false, "failure to measure stats for histogram");
+        return false;
+    }
 
     // scale the plot to hold the histogram
Index: branches/pap/psphot/src/psphotFindFootprints.c
===================================================================
--- branches/pap/psphot/src/psphotFindFootprints.c	(revision 25025)
+++ branches/pap/psphot/src/psphotFindFootprints.c	(revision 25027)
@@ -59,5 +59,5 @@
 
         psArray *mergedFootprints = pmFootprintArraysMerge(detections->footprints, footprints, includePeaks);
-        psLogMsg ("psphot", PS_LOG_MINUTIA, "merged %ld new footprints with %ld existing ones: %f sec\n", footprints->n, detections->footprints->n, psTimerMark ("psphot.footprints"));
+        psLogMsg ("psphot", PS_LOG_MINUTIA, "merged %ld new footprints with %ld existing ones: %f sec\n", footprints->n, (detections->footprints ? detections->footprints->n : 0), psTimerMark ("psphot.footprints"));
 
         psFree(footprints);
Index: branches/pap/psphot/src/psphotFitSourcesLinearStack.c
===================================================================
--- branches/pap/psphot/src/psphotFitSourcesLinearStack.c	(revision 25027)
+++ branches/pap/psphot/src/psphotFitSourcesLinearStack.c	(revision 25027)
@@ -0,0 +1,252 @@
+# include "psphotInternal.h"
+
+// fit flux (and optionally sky model) to all reasonable sources
+// with the linear fitting process.  sources must have an associated
+// model with selected pixels, and the fit radius must be defined
+
+// given the set of sources, each of which points to the pixels in the
+// science image, we construct a set of simulated sources with their own pixels.
+// these are used to determine the simultaneous linear fit of fluxes.
+// the analysis is performed wrt the simulated pixel values
+
+static bool SetBorderMatrixElements (psSparseBorder *border, pmReadout *readout, psArray *sources, bool constant_weights, psImageMaskType markVal);
+
+bool psphotFitSourcesLinear (pmReadout *readout, psArray *sources, psMetadata *recipe, pmPSF *psf, bool final) {
+
+    bool status;
+    float x;
+    float y;
+    float f;
+    // float r;
+
+    psTimerStart ("psphot.linear");
+
+    // bit-masks to test for good/bad pixels
+    psImageMaskType maskVal = psMetadataLookupImageMask(&status, recipe, "MASK.PSPHOT");
+    assert (maskVal);
+
+    // bit-mask to mark pixels not used in analysis
+    psImageMaskType markVal = psMetadataLookupImageMask(&status, recipe, "MARK.PSPHOT");
+    assert (markVal);
+
+    // maskVal is used to test for rejected pixels, and must include markVal
+    maskVal |= markVal;
+
+    // source analysis is done in spatial order
+    sources = psArraySort (sources, pmSourceSortByY);
+
+    // storage array for fitSources
+    psArray *fitSources = psArrayAllocEmpty (sources->n);
+
+    bool CONSTANT_PHOTOMETRIC_WEIGHTS =
+        psMetadataLookupBool(&status, recipe, "CONSTANT_PHOTOMETRIC_WEIGHTS");
+    if (!status) {
+        psAbort("You must provide a value for the BOOL recipe CONSTANT_PHOTOMETRIC_WEIGHTS");
+    }
+
+    // select the sources which will be used for the fitting analysis
+    for (int i = 0; i < sources->n; i++) {
+        pmSource *source = sources->data[i];
+
+        // turn this bit off and turn it on again if we pass this test
+        source->mode &= ~PM_SOURCE_MODE_LINEAR_FIT;
+
+        // skip non-astronomical objects (very likely defects)
+        if (source->type == PM_SOURCE_TYPE_DEFECT) continue;
+        if (source->type == PM_SOURCE_TYPE_SATURATED) continue;
+
+        // do not include CRs in the full ensemble fit
+        if (source->mode & PM_SOURCE_MODE_CR_LIMIT) continue;
+
+        if (final) {
+            if (source->tmpFlags & PM_SOURCE_TMPF_SUBTRACTED) continue;
+        } else {
+            if (source->mode & PM_SOURCE_MODE_BLEND) continue;
+        }
+
+        // generate model for sources without, or skip if we can't
+        if (!source->modelFlux) {
+            if (!pmSourceCacheModel (source, maskVal)) continue;
+        }
+
+        source->mode |= PM_SOURCE_MODE_LINEAR_FIT;
+        psArrayAdd (fitSources, 100, source);
+    }
+    psLogMsg ("psphot.ensemble", PS_LOG_MINUTIA, "built fitSources: %f sec (%ld objects)\n", psTimerMark ("psphot.linear"), sources->n);
+
+    if (fitSources->n == 0) {
+        return true;
+    }
+
+    // vectors to store stats for each object
+    // psVector *variance = psVectorAlloc (fitSources->n, PS_TYPE_F32);
+    psVector *errors = psVectorAlloc (fitSources->n, PS_TYPE_F32);
+
+    // create the border matrix (includes the sparse matrix)
+    // for just sky: 1 row; for x,y terms: 3 rows
+    psSparse *sparse = psSparseAlloc (fitSources->n, 100);
+    psSparseBorder *border = psSparseBorderAlloc (sparse, 1);
+
+    // fill out the sparse matrix elements and border elements (B)
+    // SRCi is the current source of interest
+    // SRCj is a possibly overlapping source
+    for (int i = 0; i < fitSources->n; i++) {
+        pmSource *SRCi = fitSources->data[i];
+
+        // diagonal elements of the sparse matrix (auto-cross-product)
+        f = pmSourceModelDotModel (SRCi, SRCi, CONSTANT_PHOTOMETRIC_WEIGHTS);
+        psSparseMatrixElement (sparse, i, i, f);
+
+        // the formal error depends on the weighting scheme
+        if (CONSTANT_PHOTOMETRIC_WEIGHTS) {
+            float var = pmSourceModelDotModel (SRCi, SRCi, false);
+            errors->data.F32[i] = 1.0 / sqrt(var);
+        } else {
+            errors->data.F32[i] = 1.0 / sqrt(f);
+        }
+
+
+        // find the image x model value
+        f = pmSourceDataDotModel (SRCi, SRCi, CONSTANT_PHOTOMETRIC_WEIGHTS);
+        psSparseVectorElement (sparse, i, f);
+
+	f = pmSourceModelWeight (SRCi, 0, CONSTANT_PHOTOMETRIC_WEIGHTS);
+	psSparseBorderElementB (border, i, 0, f);
+
+        // loop over all other stars following this one
+        for (int j = i + 1; j < fitSources->n; j++) {
+            pmSource *SRCj = fitSources->data[j];
+
+            // skip over disjoint source images, break after last possible overlap
+            if (SRCi->pixels->row0 + SRCi->pixels->numRows < SRCj->pixels->row0) break;
+            if (SRCj->pixels->row0 + SRCj->pixels->numRows < SRCi->pixels->row0) continue;
+            if (SRCi->pixels->col0 + SRCi->pixels->numCols < SRCj->pixels->col0) continue;
+            if (SRCj->pixels->col0 + SRCj->pixels->numCols < SRCi->pixels->col0) continue;
+
+            // got an overlap; calculate cross-product and add to output array
+            f = pmSourceModelDotModel (SRCi, SRCj, CONSTANT_PHOTOMETRIC_WEIGHTS);
+            psSparseMatrixElement (sparse, j, i, f);
+        }
+    }
+
+    psSparseResort (sparse);
+    psLogMsg ("psphot.ensemble", PS_LOG_MINUTIA, "built matrix: %f sec (%d elements)\n", psTimerMark ("psphot.linear"), sparse->Nelem);
+
+    // set the sky, sky_x, sky_y components of border matrix
+    SetBorderMatrixElements (border, readout, fitSources, CONSTANT_PHOTOMETRIC_WEIGHTS, SKY_FIT_ORDER, markVal);
+
+    psSparseConstraint constraint;
+    constraint.paramMin   = 0.0;
+    constraint.paramMax   = 1e8;
+    constraint.paramDelta = 1e8;
+
+    // solve for normalization terms (need include local sky?)
+    psVector *norm = psSparseSolve (NULL, constraint, sparse, 5);
+
+    psLogMsg ("psphot.ensemble", PS_LOG_MINUTIA, "solve matrix: %f sec (%d elements)\n", psTimerMark ("psphot.linear"), sparse->Nelem);
+
+    // XXXX **** philosophical question:
+    // we measure bright objects in three passes: 1) linear fit; 2) non-linear fit; 3) linear fit:
+    // should retain the chisq and errors from the intermediate non-linear fit?
+    // the non-linear fit provides better values for the position errors, and for
+    // extended sources, the shape errors
+
+    // adjust I0 for fitSources and subtract
+    for (int i = 0; i < fitSources->n; i++) {
+        pmSource *source = fitSources->data[i];
+        pmModel *model = pmSourceGetModel (NULL, source);
+
+        // assign linearly-fitted normalization
+        if (isnan(norm->data.F32[i])) {
+            psAbort("linear fitted source is nan");
+        }
+
+        model->params->data.F32[PM_PAR_I0] = norm->data.F32[i];
+        model->dparams->data.F32[PM_PAR_I0] = errors->data.F32[i];
+        // XXX is the value of 'errors' modified by the sky fit?
+
+        // subtract object
+        pmSourceSub (source, PM_MODEL_OP_FULL, maskVal);
+        source->tmpFlags |= PM_SOURCE_TMPF_SUBTRACTED;
+    }
+    psLogMsg ("psphot.ensemble", PS_LOG_MINUTIA, "sub models: %f sec (%d elements)\n", psTimerMark ("psphot.linear"), sparse->Nelem);
+
+    // measure chisq for each source
+    for (int i = 0; final && (i < fitSources->n); i++) {
+        pmSource *source = fitSources->data[i];
+        if (source->mode & PM_SOURCE_MODE_NONLINEAR_FIT) continue;
+        pmModel *model = pmSourceGetModel (NULL, source);
+        pmSourceChisq (model, source->pixels, source->maskObj, source->variance, maskVal);
+    }
+    psLogMsg ("psphot.ensemble", PS_LOG_MINUTIA, "get chisqs: %f sec (%d elements)\n", psTimerMark ("psphot.linear"), sparse->Nelem);
+
+    // psFree (index);
+    psFree (sparse);
+    psFree (fitSources);
+    psFree (norm);
+    psFree (errors);
+    psFree (border);
+
+    psLogMsg ("psphot.ensemble", PS_LOG_INFO, "measure ensemble of PSFs: %f sec\n", psTimerMark ("psphot.linear"));
+
+    psphotVisualShowResidualImage (readout);
+    psphotVisualShowFlags (sources);
+
+    return true;
+}
+
+// Calculate the weight terms for the sky fit component of the matrix.  This function operates
+// on the pixels which correspond to all of the sources of interest.  These elements fill in
+// the border matrix components in the sparse matrix equation.
+static bool SetBorderMatrixElements (psSparseBorder *border, pmReadout *readout, psArray *sources, bool constant_weights, psImageMaskType markVal) {
+
+    // generate the image-wide weight terms
+    // turn on MARK for all image pixels
+    psRegion fullArray = psRegionSet (0, 0, 0, 0);
+    fullArray = psRegionForImage (readout->mask, fullArray);
+    psImageMaskRegion (readout->mask, fullArray, "OR", markVal);
+
+    // turn off MARK for all object pixels
+    for (int i = 0; i < sources->n; i++) {
+        pmSource *source = sources->data[i];
+        pmModel *model = pmSourceGetModel (NULL, source);
+        if (model == NULL) continue;
+        float x = model->params->data.F32[PM_PAR_XPOS];
+        float y = model->params->data.F32[PM_PAR_YPOS];
+        psImageMaskCircle (source->maskView, x, y, model->radiusFit, "AND", PS_NOT_IMAGE_MASK(markVal));
+    }
+
+    // accumulate the image statistics from the masked regions
+    psF32 **image  = readout->image->data.F32;
+    psF32 **variance = readout->variance->data.F32;
+    psImageMaskType  **mask   = readout->mask->data.PS_TYPE_IMAGE_MASK_DATA;
+
+    double w, x, y, x2, xy, y2, xc, yc, wt, f, fo, fx, fy;
+    w = x = y = x2 = xy = y2 = fo = fx = fy = 0;
+
+    int col0 = readout->image->col0;
+    int row0 = readout->image->row0;
+
+    for (int j = 0; j < readout->image->numRows; j++) {
+        for (int i = 0; i < readout->image->numCols; i++) {
+            if (mask[j][i]) continue;
+            if (constant_weights) {
+                wt = 1.0;
+            } else {
+                wt = variance[j][i];
+            }
+            f = image[j][i];
+            w   += 1/wt;
+            fo  += f/wt;
+        }
+    }
+
+    // turn off MARK for all image pixels
+    psImageMaskRegion (readout->mask, fullArray, "AND", PS_NOT_IMAGE_MASK(markVal));
+
+    // set the Border T elements
+    psSparseBorderElementG (border, 0, fo);
+    psSparseBorderElementT (border, 0, 0, w);
+
+    return true;
+}
Index: branches/pap/psphot/src/psphotGuessModels.c
===================================================================
--- branches/pap/psphot/src/psphotGuessModels.c	(revision 25025)
+++ branches/pap/psphot/src/psphotGuessModels.c	(revision 25027)
@@ -6,19 +6,4 @@
 // 2) loop over the sources once and associate them with their cell
 // 3) define the threaded function to work with sources for a given cell
-
-// A guess for when the moments aren't available
-static pmModel *wildGuess(pmSource *source, // Source for which to guess
-                          pmPSF *psf    // The point-spread function
-    )
-{
-    pmModel *model = pmModelAlloc(psf->type);
-    psF32 *PAR = model->params->data.F32;
-    PAR[PM_PAR_SKY]  = 0;
-    // XXX get this from the image pixels
-    PAR[PM_PAR_I0]   = source->peak->flux;
-    PAR[PM_PAR_XPOS] = source->peak->xf;
-    PAR[PM_PAR_YPOS] = source->peak->yf;
-    return model;
-}
 
 // construct an initial PSF model for each object
@@ -81,11 +66,4 @@
 	    }
 	    psFree(job);
-
-# if (0)		
-		if (!psphotGuessModel_Unthreaded (readout, cells->data[j], psf, maskVal, markVal)) {
-		    psError(PS_ERR_UNKNOWN, false, "Unable to guess model.");
-		    return false;
-		}
-# endif
 	}
 
@@ -152,47 +130,46 @@
 	nSrc ++;
 	
-	// XXX if a source is faint, it will not have moments measured.
-	// it must be modelled as a PSF.  In this case, we need to use
-	// the peak centroid to get the coordinates and get the peak flux
-	// from the image?
-	pmModel *modelEXT;
-	if (!source->moments) {
-	    modelEXT = wildGuess(source, psf);
+	// the guess central intensity comes from the peak:
+	float Io = source->peak->flux;
+
+	// We have two options to get a guess for the object position: the position from the
+	// peak and the position from the moments.  Use the peak position if (a) there are no
+	// moments and (b) the sources is not saturated
+
+	bool useMoments = false;
+	useMoments = (source->mode & PM_SOURCE_MODE_SATSTAR);  // we only want to try if SATSTAR is set, but..
+	useMoments = (useMoments && source->moments);	       // can't if there are no moments
+	useMoments = (useMoments && source->moments->nPixels); // can't if the moments were not measured
+	useMoments = (useMoments && !(source->mode && PM_SOURCE_MODE_MOMENTS_FAILURE)); // can't if the moments failed...
+
+	float Xo, Yo;
+	if (useMoments) {
+	    Xo = source->moments->Mx;
+	    Yo = source->moments->My;
 	} else {
-	    // use the source moments, etc to guess basic model parameters
-	    modelEXT = pmSourceModelGuess (source, psf->type); // ALLOC X5
-	    if (!modelEXT) {
-		modelEXT = wildGuess(source, psf);
-	    }
-	    // these valuse are set in pmSourceModelGuess, should this rule be in there as well?
-	    if (source->mode &  PM_SOURCE_MODE_SATSTAR) {
-		modelEXT->params->data.F32[PM_PAR_XPOS] = source->moments->Mx;
-		modelEXT->params->data.F32[PM_PAR_YPOS] = source->moments->My;
-	    } else {
-		modelEXT->params->data.F32[PM_PAR_XPOS] = source->peak->xf;
-		modelEXT->params->data.F32[PM_PAR_YPOS] = source->peak->yf;
-	    }
+	    Xo = source->peak->xf;
+	    Yo = source->peak->yf;
 	}
 
-	// set PSF parameters for this model (apply 2D shape model)
-	pmModel *modelPSF = pmModelFromPSF (modelEXT, psf); // ALLOC X5
+	// set PSF parameters for this model (apply 2D shape model to coordinates Xo, Yo)
+	pmModel *modelPSF = pmModelFromPSFforXY(psf, Xo, Yo, Io);
+
 	if (modelPSF == NULL) {
-	    psWarning ("Failed to determine PSF model at r,c = (%d,%d); trying centre of image",
-		    source->peak->y, source->peak->x);
+	    psWarning ("Failed to determine PSF model at (%f,%f); trying image center", Xo, Yo);
 
-	    // Try the center of the image
-	    modelEXT->params->data.F32[PM_PAR_XPOS] = 0.5*readout->image->numCols;
-	    modelEXT->params->data.F32[PM_PAR_YPOS] = 0.5*readout->image->numRows;
-	    modelPSF = pmModelFromPSF (modelEXT, psf);
+	    float Xc = 0.5*readout->image->numCols;
+	    float Yc = 0.5*readout->image->numRows;
+	    pmModel *modelPSF = pmModelFromPSFforXY(psf, Xc, Yc, Io);
 	    if (modelPSF == NULL) {
 		psError(PSPHOT_ERR_PSF, false, "Failed to determine PSF model at center of image");
-		psFree(modelEXT);
 		return false;
 	    }
+
+	    // Now set the object position at the expected location:
+	    modelPSF->params->data.F32[PM_PAR_XPOS] = Xo;
+	    modelPSF->params->data.F32[PM_PAR_YPOS] = Yo;
 	    source->mode |= PM_SOURCE_MODE_BADPSF;
 	}
-	psFree (modelEXT); // FREE (x3)
 
-	// XXX need to define the guess flux?
 	// set the fit radius based on the object flux limit and the model
 	// this function affects the mask pixels
@@ -209,85 +186,2 @@
     return true;
 }
-
-# if (0)
-// construct models only for sources in the specified region
-bool psphotGuessModel_Unthreaded (pmReadout *readout, psArray *sources, pmPSF *psf, psImageMaskType maskVal, psImageMaskType markVal) {
-
-    int nSrc = 0;
-
-    for (int i = 0; i < sources->n; i++) {
-	pmSource *source = sources->data[i];
-
-	// XXXX this is just for a test: use this to mark sources for which the model is measured
-	// check later that all are used.
-	source->mode |= PM_SOURCE_MODE_EXT_LIMIT;
-
-	// skip non-astronomical objects (very likely defects)
-	if (source->type == PM_SOURCE_TYPE_DEFECT) continue;
-	if (source->type == PM_SOURCE_TYPE_SATURATED) continue;
-	if (!source->peak) continue;
-
-	nSrc ++;
-	
-	// XXX if a source is faint, it will not have moments measured.
-	// it must be modelled as a PSF.  In this case, we need to use
-	// the peak centroid to get the coordinates and get the peak flux
-	// from the image?
-	pmModel *modelEXT;
-	if (!source->moments) {
-	    modelEXT = wildGuess(source, psf);
-	} else {
-	    // use the source moments, etc to guess basic model parameters
-	    modelEXT = pmSourceModelGuess (source, psf->type); // ALLOC
-	    if (!modelEXT) {
-		modelEXT = wildGuess(source, psf);
-	    }
-	    // these valuse are set in pmSourceModelGuess, should this rule be in there as well?
-	    if (source->mode &  PM_SOURCE_MODE_SATSTAR) {
-		modelEXT->params->data.F32[PM_PAR_XPOS] = source->moments->Mx;
-		modelEXT->params->data.F32[PM_PAR_YPOS] = source->moments->My;
-	    } else {
-		modelEXT->params->data.F32[PM_PAR_XPOS] = source->peak->xf;
-		modelEXT->params->data.F32[PM_PAR_YPOS] = source->peak->yf;
-	    }
-	}
-
-	// set PSF parameters for this model (apply 2D shape model)
-	pmModel *modelPSF = pmModelFromPSF (modelEXT, psf); // ALLOC
-	if (modelPSF == NULL) {
-	    psError(PSPHOT_ERR_PSF, false,
-		    "Failed to determine PSF model at r,c = (%d,%d); trying centre of image",
-		    source->peak->y, source->peak->x);
-	    //
-	    // Try the centre of the image
-	    //
-	    modelEXT->params->data.F32[PM_PAR_XPOS] = 0.5*readout->image->numCols;
-	    modelEXT->params->data.F32[PM_PAR_YPOS] = 0.5*readout->image->numRows;
-	    modelPSF = pmModelFromPSF (modelEXT, psf);
-	    if (modelPSF == NULL) {
-		psError(PSPHOT_ERR_PSF, false,
-			"Failed to determine PSF model at centre of image");
-		psFree(modelEXT);
-		return false;
-	    }
-
-	    source->mode |= PM_SOURCE_MODE_BADPSF;
-	}
-	psFree (modelEXT);
-
-	// XXX need to define the guess flux?
-	// set the fit radius based on the object flux limit and the model
-	// this function affects the mask pixels
-	psphotCheckRadiusPSF (readout, source, modelPSF, markVal);
-
-	// set the source PSF model
-	source->modelPSF = modelPSF;
-	source->modelPSF->residuals = psf->residuals;
-
-	pmSourceCacheModel (source, maskVal);
-
-    }
-
-    return true;
-}
-# endif
Index: branches/pap/psphot/src/psphotImageLoop.c
===================================================================
--- branches/pap/psphot/src/psphotImageLoop.c	(revision 25025)
+++ branches/pap/psphot/src/psphotImageLoop.c	(revision 25027)
@@ -75,4 +75,11 @@
                 }
 
+                psImageMaskType maskSat = pmConfigMaskGet("SAT", config); // Mask value for saturated pixels
+                if (!pmReadoutMaskNonfinite(readout, maskSat)) {
+                    psError(psErrorCodeLast(), false, "Unable to mask non-finite pixels.");
+                    psFree(view);
+                    return false;
+                }
+
                 // run the actual photometry analysis on this chip/cell/readout
                 if (!psphotReadout (config, view)) {
Index: branches/pap/psphot/src/psphotImageQuality.c
===================================================================
--- branches/pap/psphot/src/psphotImageQuality.c	(revision 25025)
+++ branches/pap/psphot/src/psphotImageQuality.c	(revision 25027)
@@ -88,4 +88,5 @@
                      "Number of stars used for IQ measurements", M2->n);
 
+// XXX make this a recipe option
 #if (USE_SAMPLE)
     psStats *stats = psStatsAlloc (PS_STAT_SAMPLE_MEAN | PS_STAT_SAMPLE_STDEV | PS_STAT_SAMPLE_QUARTILE);
Index: branches/pap/psphot/src/psphotMakeResiduals.c
===================================================================
--- branches/pap/psphot/src/psphotMakeResiduals.c	(revision 25025)
+++ branches/pap/psphot/src/psphotMakeResiduals.c	(revision 25027)
@@ -200,6 +200,15 @@
             // measure the robust median to determine a baseline reference value
             *fluxClip = *fluxClipDef;
-            psVectorStats (fluxClip, fluxes, NULL, fmasks, fmaskVal);
-            psErrorClear();             // clear (ignore) any outstanding errors
+            if (!psVectorStats (fluxClip, fluxes, NULL, fmasks, fmaskVal)) {
+		psError(PSPHOT_ERR_CONFIG, false, "Error calculating residual stats");
+		return false;
+	    }
+	    if (isnan(fluxClip->robustMedian)) {
+                resid->Ro->data.F32[oy][ox] = 0.0;
+                resid->Rx->data.F32[oy][ox] = 0.0;
+                resid->Ry->data.F32[oy][ox] = 0.0;
+                resid->mask->data.PM_TYPE_RESID_MASK_DATA[oy][ox] = badMask;
+                continue;
+	    }
 
             // mark input pixels which are more than N sigma from the median
@@ -220,10 +229,19 @@
                 // measure the desired statistic on the unclipped pixels
                 *fluxStats = *fluxStatsDef;
-                psVectorStats (fluxStats, fluxes, NULL, fmasks, fmaskVal);
-                psErrorClear();         // clear (ignore) any outstanding errors
+                if (!psVectorStats (fluxStats, fluxes, NULL, fmasks, fmaskVal)) {
+		    psError(PSPHOT_ERR_CONFIG, false, "Error calculating residual stats");
+		    return false;
+		}
 
                 resid->Ro->data.F32[oy][ox] = psStatsGetValue(fluxStats, statOption);
                 resid->Rx->data.F32[oy][ox] = resid->Ry->data.F32[oy][ox] = 0.0;
-                //resid->variance->data.F32[oy][ox] = fluxStats->sampleStdev;
+
+		if (isnan(resid->Ro->data.F32[oy][ox])) {
+		    resid->Ro->data.F32[oy][ox] = 0.0;
+		    resid->Rx->data.F32[oy][ox] = 0.0;
+		    resid->Ry->data.F32[oy][ox] = 0.0;
+		    resid->mask->data.PM_TYPE_RESID_MASK_DATA[oy][ox] = badMask;
+		    continue;
+		}
 
                 if (fabs(resid->Ro->data.F32[oy][ox]) < pixelSN*fluxStats->sampleStdev/sqrt(nKeep)) {
Index: branches/pap/psphot/src/psphotMaskReadout.c
===================================================================
--- branches/pap/psphot/src/psphotMaskReadout.c	(revision 25025)
+++ branches/pap/psphot/src/psphotMaskReadout.c	(revision 25027)
@@ -27,5 +27,5 @@
     }
     if (!readout->variance) {
-        if (!pmReadoutGenerateVariance(readout, true)) {
+        if (!pmReadoutGenerateVariance(readout, NULL, true)) {
             psError (PSPHOT_ERR_CONFIG, false, "trouble creating variance");
             return false;
Index: branches/pap/psphot/src/psphotModelBackground.c
===================================================================
--- branches/pap/psphot/src/psphotModelBackground.c	(revision 25025)
+++ branches/pap/psphot/src/psphotModelBackground.c	(revision 25027)
@@ -189,9 +189,9 @@
 
             if (gotX && gotY) {
-                psTraceSetLevel ("psLib.math.vectorFittedStats_v4", 6);
-                psTraceSetLevel ("psLib.math.vectorRobustStats", 6);
+                (void) psTraceSetLevel ("psLib.math.vectorFittedStats_v4", 6);
+                (void) psTraceSetLevel ("psLib.math.vectorRobustStats", 6);
             } else {
-                psTraceSetLevel ("psLib.math.vectorFittedStats_v4", 0);
-                psTraceSetLevel ("psLib.math.vectorRobustStats", 0);
+                (void) psTraceSetLevel ("psLib.math.vectorFittedStats_v4", 0);
+                (void) psTraceSetLevel ("psLib.math.vectorRobustStats", 0);
             }
             # endif
Index: branches/pap/psphot/src/psphotModelTest.c
===================================================================
--- branches/pap/psphot/src/psphotModelTest.c	(revision 25025)
+++ branches/pap/psphot/src/psphotModelTest.c	(revision 25027)
@@ -139,5 +139,5 @@
 
     // get the source moments
-    status = pmSourceMoments (source, mRADIUS);
+    status = pmSourceMoments (source, mRADIUS, 0.0, 1.0);
     if (!status) psAbort("psSourceMoments error");
     source->peak->value = source->moments->Peak;
Index: branches/pap/psphot/src/psphotMomentsStudy.c
===================================================================
--- branches/pap/psphot/src/psphotMomentsStudy.c	(revision 25027)
+++ branches/pap/psphot/src/psphotMomentsStudy.c	(revision 25027)
@@ -0,0 +1,171 @@
+# include "psphotInternal.h"
+
+int SetUpKapa (Graphdata *graphdata);
+int PlotVectors (int kapa, Graphdata *graphdata, psVector *inVector, psVector *outVector);
+int PlotModelSet (int kapa, Graphdata *graphdata, pmReadout *readout, char *name, float par7);
+
+FILE *output = NULL;
+
+int main (int argc, char **argv) {
+
+  Graphdata graphdata;
+  int kapa;
+
+  pmErrorRegister();                  // register psModule's error codes/messages
+  pmModelClassInit();
+
+  // loop over sigma_window
+  // loop over sigma_input 
+  // generate a fake source using the desired model
+  // add noise?
+  // measure the moments for the object
+  // accumulate Mxx, Myy
+
+  // plot sigma_input vs sigma_output == sqrt(Mxx)
+
+  // pmSourceMoments requires a source with a peak and pixels
+
+  if (argc != 2) {
+    fprintf (stderr, "USAGE: psphotMomentsStudy (output)\n");
+    exit (2);
+  }
+
+  output = fopen (argv[1], "w");
+  if (!output) exit (1);
+  fprintf (output, "# %4s %5s  %6s %6s %6s\n", "type", "par7", "sigWin", "sigIn", "sigOut");
+
+  kapa = SetUpKapa(&graphdata);
+
+  // create a containing image & associated readout
+  pmReadout *readout = pmReadoutAlloc(NULL);
+  readout->image = psImageAlloc(128, 128, PS_TYPE_F32);
+
+  // create a dummy variance, but don't populate -- it is not used by pmSourceMoments if the sigma parameters is set to 0.0
+  readout->variance = psImageAlloc(128, 128, PS_TYPE_F32);
+
+  PlotModelSet (kapa, &graphdata, readout, "PS_MODEL_GAUSS",  0.0);
+  PlotModelSet (kapa, &graphdata, readout, "PS_MODEL_PGAUSS", 0.0);
+
+  PlotModelSet (kapa, &graphdata, readout, "PS_MODEL_PS1_V1", 0.3);
+  PlotModelSet (kapa, &graphdata, readout, "PS_MODEL_PS1_V1", 1.0);
+  PlotModelSet (kapa, &graphdata, readout, "PS_MODEL_PS1_V1", 3.0);
+
+  PlotModelSet (kapa, &graphdata, readout, "PS_MODEL_SERSIC", 0.125);
+  PlotModelSet (kapa, &graphdata, readout, "PS_MODEL_SERSIC", 0.250);
+  PlotModelSet (kapa, &graphdata, readout, "PS_MODEL_SERSIC", 0.500);
+
+  // pause and wait for user input:
+  // continue, save (provide name), ??
+  char key[10];
+  fprintf (stdout, "[c]ontinue? ");
+  if (!fgets(key, 8, stdin)) {
+    psWarning("Unable to read option");
+  }
+  fclose (output);
+
+  // psphotSaveImage(NULL, readout->image, "test.fits");
+  exit (0);
+}
+
+int SetUpKapa (Graphdata *graphdata) {
+
+  int kapa = KapaOpenNamedSocket ("kapa", "psphot");
+
+  KapaClearPlots (kapa);
+  KapaInitGraph (graphdata);
+  KapaSetFont (kapa, "courier", 14);
+
+  graphdata->color = KapaColorByName ("black");
+  graphdata->xmin = -0.1;
+  graphdata->xmax =  1.6;
+  graphdata->ymin = -0.1;
+  graphdata->ymax =  1.6;
+  // graphdata->xmax = 10.1;
+  // graphdata->ymax = 10.1;
+  KapaSetLimits (kapa, graphdata);
+
+  KapaBox (kapa, graphdata);
+  // KapaSendLabel (kapa, "&ss&h_in| (pixels)", KAPA_LABEL_XM);
+  // KapaSendLabel (kapa, "&ss&h_out| (pixels)", KAPA_LABEL_YM);
+  KapaSendLabel (kapa, "&ss&h_out| / &ss&h_window|", KAPA_LABEL_XM);
+  KapaSendLabel (kapa, "&ss&h_out| / &ss&h_in|", KAPA_LABEL_YM);
+
+  return (kapa);
+}
+
+int PlotVectors (int kapa, Graphdata *graphdata, psVector *inVector, psVector *outVector) {
+
+  graphdata->color = KapaColorByName ("black");
+  graphdata->ptype = 2;
+  graphdata->size = 1.0;
+  graphdata->style = 2;
+  KapaPrepPlot (kapa, inVector->n, graphdata);
+  KapaPlotVector (kapa, inVector->n, inVector->data.F32, "x");
+  KapaPlotVector (kapa, inVector->n, outVector->data.F32, "y");
+
+  return 0;
+}
+
+
+int PlotModelSet (int kapa, Graphdata *graphdata, pmReadout *readout, char *name, float par7) {
+
+  // create a model & associated source
+  // pmModelType type = pmModelClassGetType("PS_MODEL_GAUSS");
+
+  pmModelType type = pmModelClassGetType(name);
+  pmModel *model = pmModelAlloc(type);
+
+  // set the model parameters
+  model->params->data.F32[PM_PAR_SKY]  = 0.0;
+  model->params->data.F32[PM_PAR_I0]   = 100.0;
+  model->params->data.F32[PM_PAR_XPOS] = 64.0;
+  model->params->data.F32[PM_PAR_YPOS] = 64.0;
+
+  model->params->data.F32[PM_PAR_SXX]  = 1.0 * M_SQRT2;
+  model->params->data.F32[PM_PAR_SYY]  = 1.0 * M_SQRT2;
+  model->params->data.F32[PM_PAR_SXY]  = 0.0;
+
+  if (model->params->n > 7) {
+    model->params->data.F32[PM_PAR_7]  = par7;
+  }
+
+  pmSource *source = pmSourceFromModel(model, readout, 32.0, PM_SOURCE_TYPE_STAR);
+
+  psVector *inVector  = psVectorAllocEmpty(100, PS_TYPE_F32);
+  psVector *outVector = psVectorAllocEmpty(100, PS_TYPE_F32);
+
+  for (float sigWindow = 1.0; sigWindow < 10.0; sigWindow += 1.0) {
+
+    for (float sigIn = 0.5; sigIn < 10.0; sigIn += 0.1) {
+
+      model->params->data.F32[PM_PAR_SXX]  = sigIn * M_SQRT2;
+      model->params->data.F32[PM_PAR_SYY]  = sigIn * M_SQRT2;
+
+      // generate the modelFlux 
+      pmSourceCacheModel(source, 0);
+
+      // instantiate the source
+      pmSourceAdd(source, PM_MODEL_OP_FUNC, 0); 
+
+      pmSourceMoments (source, 32.0, sigWindow, 0.0);
+      // fprintf (stderr, "sigOut : %f\n", sqrt(source->moments->Mxx));
+
+      psVectorAppend(inVector,  sqrt(source->moments->Mxx) / sigWindow);
+      psVectorAppend(outVector, sqrt(source->moments->Mxx) / sigIn);
+
+      psImageInit (readout->image, 0.0);
+
+      fprintf (output, "  %4d %5.2f  %6.3f %6.3f %6.3f\n", type, par7, sigWindow, sigIn, sqrt(source->moments->Mxx));
+    }    
+  }
+
+  PlotVectors (kapa, graphdata, inVector, outVector);
+  inVector->n = 0;
+  outVector->n = 0;
+
+  psFree (source);
+  psFree (inVector);
+  psFree (outVector);
+
+  return 1;
+}
Index: branches/pap/psphot/src/psphotReadout.c
===================================================================
--- branches/pap/psphot/src/psphotReadout.c	(revision 25025)
+++ branches/pap/psphot/src/psphotReadout.c	(revision 25027)
@@ -61,5 +61,5 @@
     // display the backsub and backgnd images
     psphotVisualShowBackground (config, view, readout);
-
+    
     // run a single-model test if desired (exits from here if test is run)
     psphotModelTest (config, view, recipe);
@@ -87,5 +87,5 @@
     // find blended neighbors of very saturated stars
     // XXX merge this with Basic Deblend?
-    psphotDeblendSatstars (sources, recipe);
+    psphotDeblendSatstars (readout, sources, recipe);
 
     // mark blended peaks PS_SOURCE_BLEND
@@ -236,2 +236,3 @@
     return psphotReadoutCleanup(config, readout, recipe, detections, psf, sources);
 }
+
Index: branches/pap/psphot/src/psphotReadoutCleanup.c
===================================================================
--- branches/pap/psphot/src/psphotReadoutCleanup.c	(revision 25025)
+++ branches/pap/psphot/src/psphotReadoutCleanup.c	(revision 25027)
@@ -9,11 +9,11 @@
     if (psErrorCodeLast() == PSPHOT_ERR_DATA) {
         psErrorStackPrint(stderr, "Error in the psphot readout analysis");
-	psErrorClear();
-    } 
-    if (psErrorCodeLast() != PS_ERR_NONE) { 
-	psFree (psf);
-	psFree (sources);
-	psFree (detections);
-	return false;
+        psErrorClear();
+    }
+    if (psErrorCodeLast() != PS_ERR_NONE) {
+        psFree (psf);
+        psFree (sources);
+        psFree (detections);
+        return false;
     }
 
@@ -22,7 +22,7 @@
         if (!psphotPSFstats (readout, recipe, psf)) {
             psError(PSPHOT_ERR_PROG, false, "Failed to measure PSF shape parameters");
-	    psFree (psf);
-	    psFree (sources);
-	    psFree (detections);
+            psFree (psf);
+            psFree (sources);
+            psFree (detections);
             return false;
         }
@@ -32,10 +32,24 @@
         if (!psphotMomentsStats (readout, recipe, sources)) {
             psError(PSPHOT_ERR_PROG, false, "Failed to measure Moment shape parameters");
-	    psFree (psf);
-	    psFree (sources);
-	    psFree (detections);
+            psFree (psf);
+            psFree (sources);
+            psFree (detections);
             return false;
         }
     }
+
+    // Check to see if the image quality was measured
+    if (!psf) {
+        bool mdok;                      // Status of MD lookup
+        int nIQ = psMetadataLookupS32(&mdok, recipe, "IQ_NSTAR"); // Number of stars for IQ measurement
+        if (!mdok || nIQ <= 0) {
+            psError(PSPHOT_ERR_DATA, false, "Unable to measure image quality");
+            psFree (psf);
+            psFree (sources);
+            psFree (detections);
+            return false;
+        }
+    }
+
 
     // write NSTARS to the image header
@@ -48,19 +62,19 @@
     psMetadataAdd (readout->analysis, PS_LIST_TAIL, "PSPHOT.HEADER",  PS_DATA_METADATA, "header stats", header);
     if (sources) {
-	psMetadataAdd (readout->analysis, PS_LIST_TAIL, "PSPHOT.SOURCES", PS_DATA_ARRAY, "psphot sources", sources);
+        psMetadataAdd (readout->analysis, PS_LIST_TAIL, "PSPHOT.SOURCES", PS_DATA_ARRAY, "psphot sources", sources);
     }
     if (psf) {
-	// save the psf for possible output.  if there was already an entry, it was loaded from external sources
-	// the new one may have been updated or modified, so replace the existing entry.  We
-	// are required to save it on the chip, but this will cause problems if we ever want to
-	// run psphot on an unmosaiced image
-	pmCell *cell = readout->parent;
-	pmChip *chip = cell->parent;
-	psMetadataAdd (chip->analysis, PS_LIST_TAIL, "PSPHOT.PSF", PS_DATA_UNKNOWN | PS_META_REPLACE,  "psphot psf", psf);
+        // save the psf for possible output.  if there was already an entry, it was loaded from external sources
+        // the new one may have been updated or modified, so replace the existing entry.  We
+        // are required to save it on the chip, but this will cause problems if we ever want to
+        // run psphot on an unmosaiced image
+        pmCell *cell = readout->parent;
+        pmChip *chip = cell->parent;
+        psMetadataAdd (chip->analysis, PS_LIST_TAIL, "PSPHOT.PSF", PS_DATA_UNKNOWN | PS_META_REPLACE,  "psphot psf", psf);
     }
 
     if (psErrorCodeLast() != PS_ERR_NONE) {
         psErrorStackPrint(stderr, "unexpected remaining errors");
-	abort();
+        abort();
     }
 
Index: branches/pap/psphot/src/psphotReadoutMinimal.c
===================================================================
--- branches/pap/psphot/src/psphotReadoutMinimal.c	(revision 25025)
+++ branches/pap/psphot/src/psphotReadoutMinimal.c	(revision 25027)
@@ -46,5 +46,6 @@
 
     // find the detections (by peak and/or footprint) in the image.
-    pmDetections *detections = psphotFindDetections (NULL, readout, recipe);
+    pmDetections *detections = pmDetectionsAlloc(); // New detections; allocated to ensure pass=2
+    detections = psphotFindDetections(detections, readout, recipe);
     if (!detections) {
         psLogMsg ("psphot", 3, "unable to find detections in this image");
@@ -58,5 +59,5 @@
     // find blended neighbors of very saturated stars
     // XXX merge this with Basic Deblend?
-    psphotDeblendSatstars (sources, recipe);
+    psphotDeblendSatstars (readout, sources, recipe);
 
     // mark blended peaks PS_SOURCE_BLEND
Index: branches/pap/psphot/src/psphotRoughClass.c
===================================================================
--- branches/pap/psphot/src/psphotRoughClass.c	(revision 25025)
+++ branches/pap/psphot/src/psphotRoughClass.c	(revision 25027)
@@ -91,5 +91,5 @@
 	return false;
     }
-    if (!psfClump.X || !psfClump.Y) {
+    if (!psfClump.X || !psfClump.Y || isnan(psfClump.X) || isnan(psfClump.Y)) {
 	psLogMsg ("psphot", 4, "Failed to find a valid PSF clump for region %f,%f - %f,%f\n", region->x0, region->y0, region->x1, region->y1);
 	return false;
Index: branches/pap/psphot/src/psphotSourceFits.c
===================================================================
--- branches/pap/psphot/src/psphotSourceFits.c	(revision 25025)
+++ branches/pap/psphot/src/psphotSourceFits.c	(revision 25027)
@@ -224,5 +224,6 @@
 
     // recalculate the source moments using the larger extended-source moments radius
-    if (!pmSourceMoments (source, EXT_MOMENTS_RAD)) return false;
+    // at this stage, skip Gaussian windowing, and do not clip pixels by S/N
+    if (!pmSourceMoments (source, EXT_MOMENTS_RAD, 0.0, 0.0)) return false;
 
     psTrace ("psphot", 5, "trying blob...\n");
@@ -234,8 +235,10 @@
     pmSource *tmpSrc = pmSourceAlloc ();
 
+    // XXX need to handle failures better here
     pmModel *EXT = psphotFitEXT (readout, source, modelTypeEXT, maskVal, markVal);
     okEXT = psphotEvalEXT (tmpSrc, EXT);
-    chiEXT = EXT->chisq / EXT->nDOF;
-
+    chiEXT = EXT ? EXT->chisq / EXT->nDOF : NAN;
+
+    // DBL will always be defined, but DBL->data[n] might not
     psArray *DBL = psphotFitDBL (readout, source, maskVal, markVal);
     okDBL  = psphotEvalDBL (tmpSrc, DBL->data[0]);
@@ -244,19 +247,22 @@
 
     // correct first model chisqs for flux trend
+    chiDBL = NAN;
     ONE = DBL->data[0];
-    chiTrend = psPolynomial1DEval (psf->ChiTrend, ONE->params->data.F32[1]);
-    ONE->chisqNorm = ONE->chisq / chiTrend;
-
-    // save chisq for double-star/galaxy comparison
-    chiDBL = ONE->chisq / ONE->nDOF;
+    if (ONE) {
+      chiTrend = psPolynomial1DEval (psf->ChiTrend, ONE->params->data.F32[1]);
+      ONE->chisqNorm = ONE->chisq / chiTrend;
+      chiDBL = ONE->chisq / ONE->nDOF; // save chisq for double-star/galaxy comparison
+    }
 
     // correct second model chisqs for flux trend
     ONE = DBL->data[1];
-    chiTrend = psPolynomial1DEval (psf->ChiTrend, ONE->params->data.F32[1]);
-    ONE->chisqNorm = ONE->chisq / chiTrend;
+    if (ONE) {
+      chiTrend = psPolynomial1DEval (psf->ChiTrend, ONE->params->data.F32[1]);
+      ONE->chisqNorm = ONE->chisq / chiTrend;
+    }
 
     psFree (tmpSrc);
 
-    // psTraceSetLevel("psModules.objects.pmSourceFitSet", 0);
+    // (void) psTraceSetLevel("psModules.objects.pmSourceFitSet", 0);
 
     if (okEXT && okDBL) {
Index: branches/pap/psphot/src/psphotSourceStats.c
===================================================================
--- branches/pap/psphot/src/psphotSourceStats.c	(revision 25025)
+++ branches/pap/psphot/src/psphotSourceStats.c	(revision 25027)
@@ -89,25 +89,14 @@
 
             if (!psThreadJobAddPending(job)) {
-                psError(PS_ERR_UNKNOWN, false, "Unable to guess model.");
+                psError(PS_ERR_UNKNOWN, false, "Unable to launch thread job PSPHOT_SOURCE_STATS");
                 psFree (job);
                 return NULL;
             }
             psFree(job);
-
-# if (0)
-                int nfail = 0;
-                int nmoments = 0;
-                if (!psphotSourceStats_Unthreaded (&nfail, &nmoments, cells->data[j], recipe)) {
-                    psError(PS_ERR_UNKNOWN, false, "Unable to guess model.");
-                    return NULL;
-                }
-                Nfail += nfail;
-                Nmoments += nmoments;
-# endif
         }
 
         // wait for the threads to finish and manage results
         if (!psThreadPoolWait (false)) {
-            psError(PS_ERR_UNKNOWN, false, "Unable to guess model.");
+            psError(PS_ERR_UNKNOWN, false, "Failure in thread job PSPHOT_SOURCE_STATS");
             return NULL;
         }
@@ -149,9 +138,13 @@
     psMetadata *recipe              = job->args->data[1];
 
-    float INNER    = psMetadataLookupF32 (&status, recipe, "SKY_INNER_RADIUS");
-    if (!status) return false;
-    float MIN_SN   = psMetadataLookupF32 (&status, recipe, "MOMENTS_SN_MIN");
-    if (!status) return false;
-    float RADIUS   = psMetadataLookupF32 (&status, recipe, "PSF_MOMENTS_RADIUS");
+    float INNER        = psMetadataLookupF32 (&status, recipe, "SKY_INNER_RADIUS");
+    if (!status) return false;
+    float MIN_SN       = psMetadataLookupF32 (&status, recipe, "MOMENTS_SN_MIN");
+    if (!status) return false;
+    float RADIUS       = psMetadataLookupF32 (&status, recipe, "PSF_MOMENTS_RADIUS");
+    if (!status) return false;
+    float MIN_PIXEL_SN = psMetadataLookupF32 (&status, recipe, "MOMENTS_MIN_PIXEL_SN");
+    if (!status) return false;
+    float SIGMA        = psMetadataLookupF32 (&status, recipe, "MOMENTS_GAUSS_SIGMA");
     if (!status) return false;
 
@@ -202,5 +195,5 @@
 
         // measure basic source moments
-        status = pmSourceMoments (source, RADIUS);
+        status = pmSourceMoments (source, RADIUS, SIGMA, MIN_PIXEL_SN);
         if (status) {
             Nmoments ++;
@@ -212,5 +205,5 @@
         BIG_RADIUS = PS_MIN (INNER, 3*RADIUS);
         psTrace ("psphot", 4, "retrying moments for %d, %d\n", source->peak->x, source->peak->y);
-        status = pmSourceMoments (source, BIG_RADIUS);
+        status = pmSourceMoments (source, BIG_RADIUS, 3.0*SIGMA, MIN_PIXEL_SN);
         if (status) {
             source->mode |= PM_SOURCE_MODE_BIG_RADIUS;
@@ -292,5 +285,5 @@
 
         // measure basic source moments
-        status = pmSourceMoments (source, RADIUS);
+        status = pmSourceMoments (source, RADIUS, SIGMA, MIN_PIXEL_SN);
         if (status) {
             Nmoments ++;
@@ -302,5 +295,5 @@
         BIG_RADIUS = PS_MIN (INNER, 3*RADIUS);
         psTrace ("psphot", 4, "retrying moments for %d, %d\n", source->peak->x, source->peak->y);
-        status = pmSourceMoments (source, BIG_RADIUS);
+        status = pmSourceMoments (source, BIG_RADIUS, 3.0*SIGMA, MIN_PIXEL_SN);
         if (status) {
             Nmoments ++;
Index: branches/pap/psphot/src/psphotTest.c
===================================================================
--- branches/pap/psphot/src/psphotTest.c	(revision 25025)
+++ branches/pap/psphot/src/psphotTest.c	(revision 25027)
@@ -50,5 +50,5 @@
     }
 
-    psTraceSetLevel ("psLib.sys.mutex", 3);
+    (void) psTraceSetLevel ("psLib.sys.mutex", 3);
 
     int nThreads = atoi (argv[2]);
Index: branches/pap/psphot/src/psphotVisual.c
===================================================================
--- branches/pap/psphot/src/psphotVisual.c	(revision 25025)
+++ branches/pap/psphot/src/psphotVisual.c	(revision 25027)
@@ -471,4 +471,5 @@
     bool status;
     Graphdata graphdata;
+    KapaSection section;
 
     if (!pmVisualIsVisual()) return true;
@@ -485,27 +486,43 @@
     KapaClearPlots (kapa3);
     KapaInitGraph (&graphdata);
-
-    // there are N regions: use the first (guaranteed to exist) to get the overall limits
-    psMetadata *regionMD = psMetadataLookupPtr (&status, recipe, "PSF.CLUMP.REGION.000");
-
-    float psfX  = psMetadataLookupF32 (&status, regionMD, "PSF.CLUMP.X");
-    float psfY  = psMetadataLookupF32 (&status, regionMD, "PSF.CLUMP.Y");
-
-    // examine sources to set data range
-    graphdata.xmin = -0.05;
-    graphdata.ymin = -0.05;
-    graphdata.xmax = 2.0*psfX;
-    graphdata.ymax = 2.0*psfY;
-    KapaSetLimits (kapa3, &graphdata);
-
-    KapaSetFont (kapa3, "helvetica", 14);
-    KapaBox (kapa3, &graphdata);
-    KapaSendLabel (kapa3, "&ss&h_x| (pixels)", KAPA_LABEL_XM);
-    KapaSendLabel (kapa3, "&ss&h_y| (pixels)", KAPA_LABEL_YM);
-
+    KapaSetFont (kapa3, "courier", 14);
+
+    float SN_LIM = psMetadataLookupF32(&status, recipe, "PSF_SN_LIM");
+
+    // select the max psfX,Y values for the plot limits
+    float Xmin = 0.0, Xmax = 0.0;
+    float Ymin = 0.0, Ymax = 0.0;
+    {
+        int nRegions = psMetadataLookupS32 (&status, recipe, "PSF.CLUMP.NREGIONS");
+        for (int n = 0; n < nRegions; n++) {
+
+            char regionName[64];
+            snprintf (regionName, 64, "PSF.CLUMP.REGION.%03d", n);
+            psMetadata *regionMD = psMetadataLookupPtr (&status, recipe, regionName);
+
+	    float psfX = psMetadataLookupF32 (&status, regionMD, "PSF.CLUMP.X");
+            float psfY = psMetadataLookupF32 (&status, regionMD, "PSF.CLUMP.Y");
+            float psfdX = psMetadataLookupF32 (&status, regionMD, "PSF.CLUMP.DX");
+            float psfdY = psMetadataLookupF32 (&status, regionMD, "PSF.CLUMP.DY");
+
+	    float X0 = psfX - 4.0*psfdX;
+	    float X1 = psfX + 4.0*psfdX;
+	    float Y0 = psfY - 4.0*psfdY;
+	    float Y1 = psfY + 4.0*psfdY;
+
+	    if (isfinite(X0)) { Xmin = PS_MAX(Xmin, X0); }
+	    if (isfinite(X1)) { Xmax = PS_MAX(Xmax, X1); }
+	    if (isfinite(Y0)) { Ymin = PS_MAX(Ymin, Y0); }
+	    if (isfinite(Y1)) { Ymax = PS_MAX(Ymax, Y1); }
+        }
+    }
+
+    // storage vectors for data to be plotted
     psVector *xBright = psVectorAllocEmpty (sources->n, PS_TYPE_F32);
     psVector *yBright = psVectorAllocEmpty (sources->n, PS_TYPE_F32);
+    psVector *mBright = psVectorAllocEmpty (sources->n, PS_TYPE_F32);
     psVector *xFaint  = psVectorAllocEmpty (sources->n, PS_TYPE_F32);
     psVector *yFaint  = psVectorAllocEmpty (sources->n, PS_TYPE_F32);
+    psVector *mFaint  = psVectorAllocEmpty (sources->n, PS_TYPE_F32);
 
     // construct the vectors
@@ -519,19 +536,45 @@
         xFaint->data.F32[nF] = source->moments->Mxx;
         yFaint->data.F32[nF] = source->moments->Myy;
+        mFaint->data.F32[nF] = -2.5*log10(source->moments->Sum);
         nF++;
 
         // XXX make this a user-defined cutoff
-        if (source->moments->SN < 50)
+        if (source->moments->SN < SN_LIM)
             continue;
 
         xBright->data.F32[nB] = source->moments->Mxx;
         yBright->data.F32[nB] = source->moments->Myy;
+        mBright->data.F32[nB] = -2.5*log10(source->moments->Sum);
         nB++;
     }
     xFaint->n = nF;
     yFaint->n = nF;
+    mFaint->n = nF;
 
     xBright->n = nB;
     yBright->n = nB;
+    mBright->n = nB;
+
+    // three sections: MxxMyy, MagMxx, MagMyy
+
+    // first section: MxxMyy
+    section.dx = 0.75;
+    section.dy = 0.75;
+    section.x  = 0.00;
+    section.y  = 0.00;
+    section.name = psStringCopy ("MxxMyy");
+    KapaSetSection (kapa3, &section);
+    psFree (section.name);
+
+    graphdata.color = KapaColorByName ("black");
+    graphdata.xmin = Xmin;
+    graphdata.ymin = Ymin;
+    graphdata.xmax = Xmax;
+    graphdata.ymax = Ymax;
+    KapaSetLimits (kapa3, &graphdata);
+
+    KapaBox (kapa3, &graphdata);
+    KapaSendLabel (kapa3, "M_xx| (pixels)", KAPA_LABEL_XM);
+    KapaSendLabel (kapa3, "M_yy| (pixels)", KAPA_LABEL_YM);
 
     graphdata.color = KapaColorByName ("black");
@@ -551,6 +594,82 @@
     KapaPlotVector (kapa3, nB, yBright->data.F32, "y");
 
-    // XXX draw N circles to outline the clumps
+    // second section: MagMyy
+    section.dx = 0.75;
+    section.dy = 0.25;
+    section.x  = 0.00;
+    section.y  = 0.80;
+    section.name = psStringCopy ("MagMyy");
+    KapaSetSection (kapa3, &section);
+    psFree (section.name);
+
+    graphdata.color = KapaColorByName ("black");
+    graphdata.xmin = -17.1;
+    graphdata.xmax =  -7.9;
+    graphdata.ymin = Ymin;
+    graphdata.ymax = Ymax;
+    KapaSetLimits (kapa3, &graphdata);
+
+    strcpy (graphdata.labels, "0210");
+    KapaBox (kapa3, &graphdata);
+    KapaSendLabel (kapa3, "inst mag", KAPA_LABEL_XP);
+    KapaSendLabel (kapa3, "M_yy| (pixels)", KAPA_LABEL_YM);
+
+    graphdata.color = KapaColorByName ("black");
+    graphdata.ptype = 0;
+    graphdata.size = 0.3;
+    graphdata.style = 2;
+    KapaPrepPlot (kapa3, nF, &graphdata);
+    KapaPlotVector (kapa3, nF, mFaint->data.F32, "x");
+    KapaPlotVector (kapa3, nF, yFaint->data.F32, "y");
+
+    graphdata.color = KapaColorByName ("red");
+    graphdata.ptype = 2;
+    graphdata.size = 0.5;
+    graphdata.style = 2;
+    KapaPrepPlot (kapa3, nB, &graphdata);
+    KapaPlotVector (kapa3, nB, mBright->data.F32, "x");
+    KapaPlotVector (kapa3, nB, yBright->data.F32, "y");
+
+    // third section: MagMxx
+    section.dx = 0.25;
+    section.dy = 0.75;
+    section.x  = 0.80;
+    section.y  = 0.00;
+    section.name = psStringCopy ("MagMxx");
+    KapaSetSection (kapa3, &section);
+    psFree (section.name);
+
+    graphdata.color = KapaColorByName ("black");
+    graphdata.xmin = Xmin;
+    graphdata.xmax = Xmax;
+    graphdata.ymin =  -7.9;
+    graphdata.ymax = -17.1;
+    KapaSetLimits (kapa3, &graphdata);
+
+    strcpy (graphdata.labels, "2001");
+    KapaBox (kapa3, &graphdata);
+    KapaSendLabel (kapa3, "M_xx| (pixels)", KAPA_LABEL_XM);
+    KapaSendLabel (kapa3, "inst mag", KAPA_LABEL_YP);
+
+    graphdata.color = KapaColorByName ("black");
+    graphdata.ptype = 0;
+    graphdata.size = 0.3;
+    graphdata.style = 2;
+    KapaPrepPlot (kapa3, nF, &graphdata);
+    KapaPlotVector (kapa3, nF, xFaint->data.F32, "x");
+    KapaPlotVector (kapa3, nF, mFaint->data.F32, "y");
+
+    graphdata.color = KapaColorByName ("red");
+    graphdata.ptype = 2;
+    graphdata.size = 0.5;
+    graphdata.style = 2;
+    KapaPrepPlot (kapa3, nB, &graphdata);
+    KapaPlotVector (kapa3, nB, xBright->data.F32, "x");
+    KapaPlotVector (kapa3, nB, mBright->data.F32, "y");
+
+    // draw N circles to outline the clumps
     {
+        KapaSelectSection (kapa3, "MxxMyy");
+
         // draw a circle centered on psfX,Y with size of the psf limit
         psVector *xLimit  = psVectorAlloc (120, PS_TYPE_F32);
@@ -562,4 +681,10 @@
         graphdata.color = KapaColorByName ("blue");
         graphdata.style = 0;
+
+	graphdata.xmin = Xmin;
+	graphdata.ymin = Ymin;
+	graphdata.xmax = Xmax;
+	graphdata.ymax = Ymax;
+	KapaSetLimits (kapa3, &graphdata);
 
         for (int n = 0; n < nRegions; n++) {
@@ -626,6 +751,8 @@
     psFree (xBright);
     psFree (yBright);
+    psFree (mBright);
     psFree (xFaint);
     psFree (yFaint);
+    psFree (mFaint);
 
     // pause and wait for user input: