Index: /trunk/psphot/Makefile
===================================================================
--- /trunk/psphot/Makefile	(revision 5801)
+++ /trunk/psphot/Makefile	(revision 5802)
@@ -51,5 +51,6 @@
 $(SRC)/pmSourceFitFixed.$(ARCH).o  \
 $(SRC)/psSparse.$(ARCH).o            \
-$(SRC)/psImageData.$(ARCH).o
+$(SRC)/psImageData.$(ARCH).o        \
+$(SRC)/psphotApResid.$(ARCH).o
 
 PSMODULES = \
@@ -72,4 +73,5 @@
 $(SRC)/modelTestFitSource.$(ARCH).o \
 $(SRC)/modelTestArguments.$(ARCH).o \
+$(SRC)/psphotModelGroupInit.$(ARCH).o  \
 $(SRC)/psModulesUtils.$(ARCH).o	    \
 $(SRC)/psImageData.$(ARCH).o        \
@@ -111,12 +113,12 @@
 
 $(BIN)/%.$(ARCH) : $(SRC)/%.$(ARCH).o
-	@if [ ! -d $(BIN) ]; then mkdir -p $(BIN); fi
-	$(CC) $^ -o $@ $(LFLAGS)
+	@if [ ! -d $(BIN) ]; then mkdir -p $(BIN) || exit; fi
+	$(CC) $^ -o $@ $(LFLAGS) || exit
 	@echo "done with $@"
 
 $(DESTBIN)/%: $(BIN)/%.$(ARCH)
-	@if [ ! -d $(DESTBIN) ]; then mkdir -p $(DESTBIN); fi
-	rm -f $(DESTBIN)/$*
-	cp $(BIN)/$*.$(ARCH) $(DESTBIN)/$*
+	@if [ ! -d $(DESTBIN) ]; then mkdir -p $(DESTBIN) || exit; fi
+	rm -f $(DESTBIN)/$* || exit
+	cp $(BIN)/$*.$(ARCH) $(DESTBIN)/$* || exit
 
 $(INSTALL): % : $(BIN)/%.$(ARCH)
@@ -127,10 +129,10 @@
 
 %.install:
-	make $(DESTBIN)/$*
+	make $(DESTBIN)/$* || exit
 
 # utilities #################################################
 
 install:
-	for i in $(INSTALL); do make $$i.install; done
+	for i in $(INSTALL); do make $$i.install || exit; done
 
 clean:	
Index: /trunk/psphot/src/modelTestFitSource.c
===================================================================
--- /trunk/psphot/src/modelTestFitSource.c	(revision 5801)
+++ /trunk/psphot/src/modelTestFitSource.c	(revision 5802)
@@ -67,4 +67,6 @@
     source->peak->counts = source->moments->Peak;
 
+    fprintf (stderr, "sum: %f\n", source->moments->Sum);
+
     // get the initial model parameter guess
     pmModel *model = pmSourceModelGuess (source, modelType);
@@ -73,10 +75,24 @@
     psF32 *params = model->params->data.F32;
     for (int i = 0; i < nParams; i++) {
+	if (i == 2) {
+	    params[i] = xObj;
+	    continue;
+	}
+	if (i == 3) {
+	    params[i] = yObj;
+	    continue;
+	}
 	sprintf (name, "TEST_FIT_PAR%d", i);
 	value = psMetadataLookupF32 (&status, config, name);
 	if (status) {
+	    fprintf (stderr, "using supplied value %f for PAR %d\n", value, i);
 	    params[i] = value;
+	} else {
+	    fprintf (stderr, "using guessed  value %f for PAR %d\n", params[i], i);
 	}
     }
+
+    float area = params[4]*params[5];
+    fprintf (stderr, "peak: %f\n", source->moments->Sum*area);
 
     // what fitting mode to use?
@@ -98,4 +114,8 @@
     pmSourceSubModel (source->pixels, source->mask, model, false, false);
     
+    for (int i = 0; i < nParams; i++) {
+	fprintf (stderr, "result value %f for PAR %d\n", params[i], i);
+    }
+
     // write out 
     DumpImage (source->pixels, "resid.fits");
Index: /trunk/psphot/src/psModulesUtils.c
===================================================================
--- /trunk/psphot/src/psModulesUtils.c	(revision 5801)
+++ /trunk/psphot/src/psModulesUtils.c	(revision 5802)
@@ -65,6 +65,12 @@
     float sky = model->params->data.F32[0];
 
+    *fitMag = 99.0;
+    *obsMag = 99.0;
+
     pmModelFlux modelFluxFunc = pmModelFlux_GetFunction (model->type);
     fitSum = modelFluxFunc (model->params);
+    if (fitSum <= 0) return false;
+    if (!isfinite(fitSum)) return false;
+    *fitMag = -2.5*log10(fitSum);
 
     for (int ix = 0; ix < image->numCols; ix++) {
@@ -75,8 +81,6 @@
     }
     if (obsSum <= 0) return false;
-    if (fitSum <= 0) return false;
+    *obsMag = -2.5*log10(obsSum);
 
-    *fitMag = -2.5*log10(fitSum);
-    *obsMag = -2.5*log10(obsSum);
     return (true);
 } 
Index: /trunk/psphot/src/psphot.c
===================================================================
--- /trunk/psphot/src/psphot.c	(revision 5801)
+++ /trunk/psphot/src/psphot.c	(revision 5802)
@@ -46,13 +46,4 @@
     pmSourceRoughClass (sources, config, psfClump);
 
-    // optional dump of all rough source data
-    char *output = psMetadataLookupPtr (&status, config, "MOMENTS_OUTPUT_FILE");
-    if (status && (output != NULL) && (output[0]))
-    {
-      pmMomentsWriteText (sources, output);
-      psFree (output);
-    }
-    if (!strcasecmp (breakPt, "CLASS")) exit (0);
-
     // mark blended peaks PS_SOURCE_BLEND
     psphotBasicDeblend (sources, config, sky);
@@ -61,4 +52,7 @@
     // source analysis is done in S/N order (brightest first)
     sources = psArraySort (sources, psphotSortBySN);
+
+    psphotDumpMoments (config, sources);
+    if (!strcasecmp (breakPt, "CLASS")) exit (0);
 
     // use bright stellar objects to measure PSF
@@ -75,9 +69,14 @@
       case 0:
 	psphotEnsemblePSF (imdata, config, sources, psf, sky);
-	psphotFullFit (imdata, config, sources, psf, sky);
-//	psphotReplaceUnfit (sources);
 	break;
 
       case 1:
+	psphotEnsemblePSF (imdata, config, sources, psf, sky);
+	psphotFullFit (imdata, config, sources, psf, sky);
+	psphotReplaceUnfit (sources);
+	psphotApResid (sources, config, psf);
+	break;
+
+      case 2:
 	// fit extended objects with galaxy models
 	psphotApplyPSF (imdata, config, sources, psf, sky);
@@ -85,10 +84,10 @@
 	break;
 
-      case 2:
+      case 3:
 	// fit extended objects with galaxy models
 	psphotApplyPSF (imdata, config, sources, psf, sky);
 	break;
 
-      case 3:
+      case 4:
 	// fit extended objects with galaxy models
 	psphotFixedPSF (imdata, config, sources, psf, sky);
Index: /trunk/psphot/src/psphot.h
===================================================================
--- /trunk/psphot/src/psphot.h	(revision 5801)
+++ /trunk/psphot/src/psphot.h	(revision 5802)
@@ -85,2 +85,4 @@
 bool psphotSamplePSFs (pmPSF *psf, psImage *image);
 bool psphotReplaceUnfit (psArray *sources);
+bool psphotDumpMoments (psMetadata *config, psArray *sources);
+bool psphotApResid (psArray *sources, psMetadata *config, pmPSF *psf);
Index: /trunk/psphot/src/psphotApResid.c
===================================================================
--- /trunk/psphot/src/psphotApResid.c	(revision 5802)
+++ /trunk/psphot/src/psphotApResid.c	(revision 5802)
@@ -0,0 +1,85 @@
+# include "psphot.h"
+
+bool psphotApResid (psArray *sources, psMetadata *config, pmPSF *psf) { 
+
+    int Npsf;
+    bool status;
+    float x, y;
+    float fitMag, obsMag;
+
+    // XXX EAM : check that PSF_FIT_RADIUS < SKY_OUTER_RADIUS
+    float RADIUS = psMetadataLookupF32 (&status, config, "PSF_FIT_RADIUS");
+
+    psTimerStart ("psphot");
+
+    // XXX EAM : drop this if we know the list is sorted
+    sources = psArraySort (sources, psphotSortBySN);
+
+    psVector *mask    = psVectorAlloc (300, PS_TYPE_U8);
+    psVector *xPos    = psVectorAlloc (300, PS_TYPE_F64);
+    psVector *yPos    = psVectorAlloc (300, PS_TYPE_F64);
+    psVector *rflux   = psVectorAlloc (300, PS_TYPE_F64);
+    psVector *apResid = psVectorAlloc (300, PS_TYPE_F64);
+    mask->n = rflux->n = apResid->n = 0;
+    Npsf = 0;
+
+    // select the NNN brightest, non-saturated sources, or just select PSFSTARs?
+    for (int i = 0; (i < sources->n) && (Npsf < 300); i++) {
+	pmSource *source = sources->data[i];
+	pmModel  *model  = source->modelPSF;
+
+	if (model == NULL) continue;
+	if (source->type != PM_SOURCE_STAR) continue; 
+	if (source->mode &  PM_SOURCE_SATSTAR) continue;
+	if (source->mode &  PM_SOURCE_BLEND) continue;
+	if (source->mode &  PM_SOURCE_FAIL) continue;
+	if (source->mode &  PM_SOURCE_POOR) continue;
+
+	// set shape for this model based on PSF
+	x = model->params->data.F32[2];
+	y = model->params->data.F32[3];
+
+	pmSourceAddModel (source->pixels, source->mask, model, false, false);
+
+	psImageKeepCircle (source->mask, x, y, RADIUS, "OR", PSPHOT_MASK_MARKED);
+	// XXX EAM : add in source flux
+	status = pmSourcePhotometry (&fitMag, &obsMag, model, source->pixels, source->mask);
+	psImageKeepCircle (source->mask, x, y, RADIUS, "AND", ~PSPHOT_MASK_MARKED);
+	if (!status) continue;
+
+	mask->data.U8[Npsf] = 0;
+	xPos->data.F64[Npsf] = x;
+	yPos->data.F64[Npsf] = y;
+	rflux->data.F64[Npsf] = pow(10.0, 0.4*fitMag);
+	apResid->data.F64[Npsf] = obsMag - fitMag;
+
+	psVectorExtend (mask, 100, 1);
+	psVectorExtend (rflux, 100, 1);
+	psVectorExtend (apResid, 100, 1);
+	Npsf ++;
+
+	pmSourceSubModel (source->pixels, source->mask, model, false, false);
+	psMemCheckCorruption (true);
+    }
+
+    // 3hi/1lo sigma clipping on the rflux vs metric fit
+    psStats *stats = psStatsAlloc (PS_STAT_SAMPLE_MEDIAN | PS_STAT_SAMPLE_STDEV);
+    stats->min = 1.0;
+    stats->max = 3.0;
+    stats->clipIter = 3;
+
+    // linear clipped fit of apResid to rflux, xPos, yPos
+    psf->ApTrend  = psVectorClipFitPolynomial3D (psf->ApTrend, stats, mask, PSFTRY_MASK_ALL, apResid, NULL, rflux, xPos, yPos);
+    psf->skyBias  = psf->ApTrend->coeff[0][0][1] / (M_PI * PS_SQR(RADIUS));
+    psf->ApResid  = psf->ApTrend->coeff[0][0][0];
+    psf->dApResid = stats->sampleStdev;
+
+    psLogMsg ("ApResid", 4, "aperture residual: %f +/- %f : %f bias\n", psf->ApResid, psf->dApResid, psf->skyBias);
+
+    psFree (stats);
+    psFree (mask);
+    psFree (rflux);
+    psFree (apResid);
+
+    return true;
+}
Index: /trunk/psphot/src/psphotChoosePSF.c
===================================================================
--- /trunk/psphot/src/psphotChoosePSF.c	(revision 5801)
+++ /trunk/psphot/src/psphotChoosePSF.c	(revision 5802)
@@ -66,11 +66,19 @@
 	}
     }
-    // XXX EAM : need to unflag the PSF stars which are not used to build the PSF
-    // XXX EAM : each pmPSFtry needs to have its own mask array
+    
+    // use the best model:
+    try = models->data[bestN];
 
-    // keep only the selected model:
-    try = models->data[bestN];
+    // unset the PSFSTAR flag for stars not used for PSF model
+    for (int i = 0; i < try->sources->n; i++) {
+	pmSource *source = try->sources->data[i];
+	if (try->mask->data.U8[i] & PSFTRY_MASK_FLT_FAIL) {
+	    source->mode &= ~PM_SOURCE_PSFSTAR;
+	}
+    }
+
+    // save only the best model;
     pmPSF *psf = psMemCopy(try->psf);
-    psFree (models);				 // keep only the pmPSF resulting from this analysis
+    psFree (models);
 
     modelName = pmModelGetType (psf->type);
Index: /trunk/psphot/src/psphotEnsemblePSF.c
===================================================================
--- /trunk/psphot/src/psphotEnsemblePSF.c	(revision 5801)
+++ /trunk/psphot/src/psphotEnsemblePSF.c	(revision 5802)
@@ -6,5 +6,4 @@
     float x;
     float y;
-    float Sky;
 
     psTimerStart ("psphot");
@@ -14,4 +13,5 @@
 
     float OUTER_RADIUS     = psMetadataLookupF32 (&status, config, "SKY_OUTER_RADIUS");
+    float INNER_RADIUS     = psMetadataLookupF32 (&status, config, "SKY_INNER_RADIUS");
     float PSF_FIT_NSIGMA   = psMetadataLookupF32 (&status, config, "PSF_FIT_NSIGMA");
     float PSF_FIT_PADDING  = psMetadataLookupF32 (&status, config, "PSF_FIT_PADDING");
@@ -38,7 +38,19 @@
 	pmSource *otSource = pmSourceAlloc ();
 
+	// really saturated stars should be re-measured for a better centroid
+	// XXX EAM : move this to a 'clear satstar function'
+	// XXX EAM : extend size of fit box around SATSTAR
+	if (inSource->mode &  PM_SOURCE_SATSTAR) {
+	    status = pmSourceMoments (inSource, INNER_RADIUS);
+	}
+
 	// XXX EAM : add option to use FLT or PSF form
 	// use the source moments, etc to guess basic model parameters
 	pmModel *modelFLT = pmSourceModelGuess (inSource, psf->type); 
+	if (inSource->mode &  PM_SOURCE_SATSTAR) {
+	    modelFLT->params->data.F32[2] = inSource->moments->x;
+	    modelFLT->params->data.F32[3] = inSource->moments->y;
+	}
+	// XXX EAM : add option to peak-up on peak (for non-sat objects)
 
 	// set PSF parameters for this model
@@ -69,9 +81,4 @@
 	psImage *flux = otSource->pixels;
 	psImage *mask = otSource->mask;
-
-	// XXX EAM : use these lines to fit to the peak
-	// model->params->data.F32[2] = inSource->peak->x;
-	// model->params->data.F32[3] = inSource->peak->y;
-	// XXX EAM : better option: improve the peak with 2D poly fit 3x3
 
 	// set model to unit peak, zero sky (we assume sky is constant)
@@ -136,11 +143,10 @@
 
 	Fi->modelPSF = Mi->modelPSF;
+
+	// assign linearly-fitted normalization
 	Fi->modelPSF->params->data.F32[1] = norm->data.F32[i];
 
 	// subtract object
 	pmSourceSubModel (Fi->pixels, Fi->mask, Fi->modelPSF, false, false);
-
-	Fi->modelPSF->params->data.F32[0] = Sky; 
-	// need to set this!
     }
 
Index: /trunk/psphot/src/psphotEvalFLT.c
===================================================================
--- /trunk/psphot/src/psphotEvalFLT.c	(revision 5801)
+++ /trunk/psphot/src/psphotEvalFLT.c	(revision 5802)
@@ -37,5 +37,5 @@
 
     // if the object has a fitted peak below 0, the fit did not converge cleanly
-    if (model->params->data.F32[1] < 0) {
+    if (model->params->data.F32[1] <= 0) {
 	source->mode |= PM_SOURCE_FAIL;
 	return false;
Index: /trunk/psphot/src/psphotEvalPSF.c
===================================================================
--- /trunk/psphot/src/psphotEvalPSF.c	(revision 5801)
+++ /trunk/psphot/src/psphotEvalPSF.c	(revision 5802)
@@ -94,5 +94,5 @@
 
     // if the object has a fitted peak below 0, the fit did not converge cleanly
-    if (model->params->data.F32[1] < 0) {
+    if (model->params->data.F32[1] <= 0) {
 	source->mode |= PM_SOURCE_FAIL;
 	return false;
@@ -125,10 +125,11 @@
     if (keep) return true;
 
-    // this source is not a star, unflag as PSFSTAR
-    // XXX : if this object was used to build the PSF, this flag should
-    //       be set even if the object is not a star...
+    // this source is not a star, warn if it was a PSFSTAR
     if (source->mode & PM_SOURCE_PSFSTAR) {
-	source->mode &= ~PM_SOURCE_PSFSTAR;
-	psLogMsg ("psphot", 5, "PSFSTAR demoted based on fit quality\n");
+	psphotSaveImage (NULL, source->pixels, "failpx.fits");
+	psphotSaveImage (NULL, source->mask, "failmk.fits");
+	psphotSaveImage (NULL, source->weight, "failwt.fits");
+	psLogMsg ("psphot", 5, "PSFSTAR demoted based on fit quality   (%f, %f  :  %f %f %f %f)\n", 
+		  model->params->data.F32[2], model->params->data.F32[3], nSx, nSy, SN, Chi);
     }
 
Index: /trunk/psphot/src/psphotFullFit.c
===================================================================
--- /trunk/psphot/src/psphotFullFit.c	(revision 5801)
+++ /trunk/psphot/src/psphotFullFit.c	(revision 5802)
@@ -1,7 +1,13 @@
 # include "psphot.h"
+
+# define CLEAR_TRACE \
+  if (TEST_SRC_ON) { \
+    psTraceSetLevel (".pmObjects.pmSourceFitModel", 0); \
+    psTraceSetLevel ("psMinimizeLMChi2", 0); \
+    TEST_SRC_ON = false; } \
 
 bool psphotFullFit (eamReadout *imdata, psMetadata *config, psArray *sources, pmPSF *psf, psStats *sky) 
 { 
-    bool  skip, status;
+    bool  status;
     float x, y;
     int   Nfit = 0;
@@ -26,4 +32,10 @@
     psphotInitRadiusFLT (config, sky, modelTypeFLT);
 
+    float TEST_SRC_X = psMetadataLookupF32 (&status, config, "TEST_SRC_X");
+    float TEST_SRC_Y = psMetadataLookupF32 (&status, config, "TEST_SRC_Y");
+    bool TEST_SRC_ON = false;
+    bool TEST_SRC    = status;
+    fprintf (stderr, "test src: %f %f\n", TEST_SRC_X, TEST_SRC_Y);
+
     for (int i = 0; i < sources->n; i++) {
 
@@ -37,5 +49,5 @@
 
 	// save the PSF model from the Ensemble fit
-	pmModel *save  = pmModelCopy (source->modelPSF);
+	pmModel *LIN  = pmModelCopy (source->modelPSF);
 
 	// attempt to fit the PSF model
@@ -49,4 +61,10 @@
 	x = modelPSF->params->data.F32[2];
 	y = modelPSF->params->data.F32[3];
+
+	if (TEST_SRC && (fabs(TEST_SRC_X - x) < 5) && (fabs(TEST_SRC_Y - y) < 5)) {
+	    psTraceSetLevel (".pmObjects.pmSourceFitModel", 6);
+	    psTraceSetLevel ("psMinimizeLMChi2", 6);
+	    TEST_SRC_ON = true;
+	}
 
 	// fit PSF model (set/unset the pixel mask)
@@ -63,20 +81,13 @@
 	    pmSourceSubModel (source->pixels, source->mask, source->modelPSF, false, false);
 	    source->mode |= PM_SOURCE_SUBTRACTED;
-	    psFree (save);
+	    psFree (LIN);
 	    Nsub ++;
+	    CLEAR_TRACE;
 	    continue;
 	} 
 
 	// skip the source if we don't think it is extended
-	skip  = (source->type != PM_SOURCE_GALAXY);
-	skip |= (source->moments->SN < GAL_MIN_SN);
-	if (skip) {
-	    source->modelPSF = save;
-	    pmSourceSubModel (source->pixels, source->mask, source->modelPSF, false, false);
-	    source->mode |= PM_SOURCE_SUBTRACTED;
-	    psFree (modelPSF);
-	    Nsub ++;
-	    continue;
-	}
+	if (source->type != PM_SOURCE_GALAXY) goto subLINEAR;
+	if (source->moments->SN < GAL_MIN_SN) goto subLINEAR;
 
 	// add the double-PSF fit mode
@@ -84,8 +95,5 @@
 	
 	// recalculate the source moments using the larger galaxy moments radius
-	if (!pmSourceMoments (source, GAL_MOMENTS_RAD)) {
-	    source->mode |= PM_SOURCE_FAIL;  // better choice?
-	    continue;
-	}
+	if (!pmSourceMoments (source, GAL_MOMENTS_RAD)) goto subLINEAR;
 
 	// use the source moments, etc to guess basic model parameters
@@ -111,16 +119,19 @@
 	    pmSourceSubModel (source->pixels, source->mask, source->modelFLT, false, false);
 	    source->mode |= PM_SOURCE_SUBTRACTED;
-	    psFree (save);
+	    psFree (LIN);
 	    Nsub ++;
+	    CLEAR_TRACE;
 	    continue;
 	}
 
 	// subtract PSF for object, leave local sky
-	// XXX : we should keep the non-linear PSF fit if object was marked extended but
-	//       PSF fit is OK otherwise.
+    subLINEAR:
 	psFree (source->modelPSF);
-	source->modelPSF = save;
+	source->modelPSF = LIN;
 	pmSourceSubModel (source->pixels, source->mask, source->modelPSF, false, false);
 	source->mode |= PM_SOURCE_SUBTRACTED;
+	source->mode |= PM_SOURCE_LINEAR;
+	Nsub ++;
+	CLEAR_TRACE;
     }
 
Index: /trunk/psphot/src/psphotMagnitudes.c
===================================================================
--- /trunk/psphot/src/psphotMagnitudes.c	(revision 5801)
+++ /trunk/psphot/src/psphotMagnitudes.c	(revision 5802)
@@ -8,5 +8,5 @@
     int status;
     float x, y;
-    float rflux, apMag, fitMag;
+    float rflux;
     pmModel *model;
 
Index: /trunk/psphot/src/psphotModelTest.c
===================================================================
--- /trunk/psphot/src/psphotModelTest.c	(revision 5801)
+++ /trunk/psphot/src/psphotModelTest.c	(revision 5802)
@@ -2,4 +2,6 @@
 
 int main (int argc, char **argv) {
+
+    psphotModelGroupInit ();
 
     psMetadata *config = modelTestArguments (&argc, argv);
Index: /trunk/psphot/src/psphotOutput.c
===================================================================
--- /trunk/psphot/src/psphotOutput.c	(revision 5801)
+++ /trunk/psphot/src/psphotOutput.c	(revision 5802)
@@ -370,5 +370,5 @@
 bool pmModelWritePSFs (psArray *sources, psMetadata *config, char *filename, pmPSF *psf) {
 
-    double dP, dmag;
+    double dPos, dMag;
     int i, j;
     FILE *f;
@@ -393,5 +393,5 @@
 	// if (source->mode & PM_SOURCE_POOR) continue;
 
-	model = pmSourceMagnitudes (source, psf, RADIUS, true);
+	model = pmSourceMagnitudes (source, psf, RADIUS);
 	if (model == NULL) continue;
 
@@ -399,14 +399,10 @@
 	dPAR = model->dparams->data.F32;
 	
-	// dP is positional error
-	dP = 0;
-	dP += PS_SQR(dPAR[2]);
-	dP += PS_SQR(dPAR[3]);
-	dP = sqrt (dP);
-	
-	dmag = dPAR[1] / PAR[1];
-
-	fprintf (f, "%7.1f %7.1f  %5.1f %7.4f  %7.4f %7.4f  ", 
-		 PAR[2], PAR[3], PAR[0], source->fitMag, dmag, dP);
+	// dPos is positional error, dMag is mag error
+	dPos = hypot (dPAR[2], dPAR[3]);
+	dMag = dPAR[1] / PAR[1];
+
+	fprintf (f, "%7.1f %7.1f  %5.1f %8.4f  %7.4f %7.4f  ", 
+		 PAR[2], PAR[3], PAR[0], source->fitMag, dMag, dPos);
 
 	for (j = 4; j < model->params->n; j++) {
@@ -417,5 +413,5 @@
 	    fprintf (f, "%9.6f ", dPAR[j]);
 	}
-	fprintf (f, ": %2d %2d %7.3f %7.3f %7.2f %4d %2d\n", 
+	fprintf (f, ": %2d %#5x %7.3f %7.1f %7.2f %4d %2d\n", 
 		 source[0].type, source[0].mode, 
 		 log10(model[0].chisq/model[0].nDOF), 
@@ -432,9 +428,8 @@
 bool pmModelWriteFLTs (psArray *sources, char *filename) {
 
-    double dP;
+    double dPos, dMag;
     int i, j;
     FILE *f;
-    psVector *params;
-    psVector *dparams;
+    psF32 *PAR, *dPAR;
     pmModel  *model;
 
@@ -450,25 +445,22 @@
 
 	if (source->type != PM_SOURCE_GALAXY) continue;
-	// if (source->mode & PM_SOURCE_FAIL) continue;
-	// if (source->mode & PM_SOURCE_POOR) continue;
 	
-	model = pmSourceMagnitudes (source->modelFLT, NULL, model->radius);
-	if (model == NULL) continue;
+	if (source->modelFLT == NULL) continue;
+	model = pmSourceMagnitudes (source, NULL, source->modelFLT->radius);
 
 	PAR  = model->params->data.F32;
 	dPAR = model->dparams->data.F32;
 	
-	// dP is shape error
+	// dPos is shape error
 	// XXX these are hardwired for SGAUSS
-	dP = 0;
-	dP += PS_SQR(dPAR[4] / PAR[4]);
-	dP += PS_SQR(dPAR[5] / PAR[5]);
-	dP += PS_SQR(dPAR[7] / PAR[7]);
-	dP = sqrt (dP);
-
-	dmag = dPAR[1] / PAR[1];
-
-	fprintf (f, "%7.1f %7.1f  %5.1f %7.3f  %7.4f %7.4f  ", 
-		 PAR[2], PAR[3], PAR[0], source->fitMag, dmag, dP);
+	dPos = 0;
+	dPos += PS_SQR(dPAR[4] / PAR[4]);
+	dPos += PS_SQR(dPAR[5] / PAR[5]);
+	// dPos += PS_SQR(dPAR[7] / PAR[7]);
+	dPos = sqrt (dPos);
+	dMag = dPAR[1] / PAR[1];
+
+	fprintf (f, "%7.1f %7.1f  %7.1f %7.3f  %7.4f %7.4f  ", 
+		 PAR[2], PAR[3], PAR[0], source->fitMag, dMag, dPos);
 
 	for (j = 4; j < model->params->n; j++) {
@@ -479,5 +471,5 @@
 	    fprintf (f, "%9.6f ", dPAR[j]);
 	}
-	fprintf (f, ": %2d %2d %7.3f %7.3f %7.2f %4d %2d\n", 
+	fprintf (f, ": %2d %#5x %7.3f %7.1f %7.2f %4d %2d\n", 
 		 source[0].type, source[0].mode,
 		 log10(model[0].chisq/model[0].nDOF), 
@@ -548,9 +540,10 @@
     }
 
+    fprintf (stderr, "writing out moments\n");
     for (i = 0; i < sources->n; i++) {
         source = sources->data[i];
         if (source->moments == NULL)
             continue;
-        fprintf (f, "%5d %5d  %7.1f  %7.1f %7.1f  %6.3f %6.3f  %8.1f %7.1f %7.1f %7.1f  %4d %2d\n",
+        fprintf (f, "%5d %5d  %7.1f  %7.1f %7.1f  %6.3f %6.3f  %10.1f %7.1f %7.1f %7.1f  %4d %2d %#5x\n",
                  source->peak->x, source->peak->y, source->peak->counts,
                  source->moments->x, source->moments->y,
@@ -558,5 +551,5 @@
                  source->moments->Sum, source->moments->Peak,
                  source->moments->Sky, source->moments->SN,
-                 source->moments->nPixels, source->type);
+                 source->moments->nPixels, source->type, source->mode);
     }
     fclose (f);
@@ -599,5 +592,6 @@
 
     psRegion region = {x - dx, x + dx, y - dy, y + dy};
-    psImage *sample = psImageSubset (image, region);
+    psImage *view = psImageSubset (image, region);
+    psImage *sample = psImageCopy (NULL, view, PS_TYPE_F32);
     psImageInit (sample, 0);
     modelFLT->params->data.F32[2] = x;
@@ -635,2 +629,17 @@
     return (TRUE);
 }
+
+bool psphotDumpMoments (psMetadata *config, psArray *sources) {
+
+    bool status;
+
+    // optional dump of all rough source data
+    char *output = psMetadataLookupPtr (&status, config, "MOMENTS_OUTPUT_FILE");
+    if (!status) return false;
+    if (output == NULL) return false;
+    if (output[0] == 0) return false;
+
+    pmMomentsWriteText (sources, output);
+    psFree (output);
+    return true;
+}
Index: /trunk/psphot/src/psphotRadiusChecks.c
===================================================================
--- /trunk/psphot/src/psphotRadiusChecks.c	(revision 5801)
+++ /trunk/psphot/src/psphotRadiusChecks.c	(revision 5802)
@@ -29,4 +29,8 @@
     if (isnan(model->radius)) psAbort ("apply_psf_model", "error in radius");
 	
+    if (source->mode &  PM_SOURCE_SATSTAR) {
+	model->radius *= 2;
+    }
+
     // check if we need to redefine the pixels
     // XXX EAM : a better test would examine the source pixels
Index: /trunk/psphot/src/psphotReplaceUnfit.c
===================================================================
--- /trunk/psphot/src/psphotReplaceUnfit.c	(revision 5801)
+++ /trunk/psphot/src/psphotReplaceUnfit.c	(revision 5802)
@@ -3,53 +3,41 @@
 bool psphotReplaceUnfit (psArray *sources) { 
 
-    bool status;
-    float x;
-    float y;
+    pmModel *model;
+    pmSource *source;
 
     psTimerStart ("psphot");
 
     for (int i = 0; i < sources->n; i++) {
-	pmSource *source = sources->data[i];
+      source = sources->data[i];
 
 	// skip non-astronomical objects (very likely defects)
 	// these were never fitted and subtracted
 	if (source->mode &  PM_SOURCE_BLEND) continue;
-	if (source->type == PM_SOURCE_DEFECT) continue; 
 	if (source->type == PM_SOURCE_SATURATED) continue;
+
+	if (source->type == PM_SOURCE_DEFECT) {
+	    if (source->mode &  PM_SOURCE_SUBTRACTED) goto addPSF;
+	    continue; 
+	}
+
+	if (source->type == PM_SOURCE_GALAXY) {
+	    if (source->mode & PM_SOURCE_FAIL) goto addPSF;
+	    if (source->mode & PM_SOURCE_POOR) goto addPSF;
+	    continue;
+	}
 
 	// need to skip the successful fits
 	if (source->type == PM_SOURCE_STAR) {
-	    if (source->mode & PM_SOURCE_FAIL) goto subPSF;
-	    if (source->mode & PM_SOURCE_POOR) goto subPSF;
-	    continue;
-
-	subPSF:
-	    pmModel  *model = source->modelPSF;
-
-	    x = model->params->data.F32[2];
-	    y = model->params->data.F32[3];
-
-	    psImageKeepCircle (source->mask, x, y, model->radius, "OR", PSPHOT_MASK_MARKED);
-	    pmSourceAddModel (source->pixels, source->mask, model, false, false);
-	    psImageKeepCircle (source->mask, x, y, model->radius, "AND", ~PSPHOT_MASK_MARKED);
+	    if (source->mode & PM_SOURCE_FAIL) goto addPSF;
+	    if (source->mode & PM_SOURCE_POOR) goto addPSF;
 	    continue;
 	}
 
-	if (source->type == PM_SOURCE_GALAXY) {
-	    if (source->mode & PM_SOURCE_FAIL) goto subFLT;
-	    if (source->mode & PM_SOURCE_POOR) goto subFLT;
-	    continue;
+    addPSF:
+	model = source->modelPSF;
+	if (model == NULL) continue;
 
-	subPSF:
-	    pmModel  *model = source->modelFLT;
-
-	    x = model->params->data.F32[2];
-	    y = model->params->data.F32[3];
-
-	    psImageKeepCircle (source->mask, x, y, model->radius, "OR", PSPHOT_MASK_MARKED);
-	    pmSourceAddModel (source->pixels, source->mask, model, false, false);
-	    psImageKeepCircle (source->mask, x, y, model->radius, "AND", ~PSPHOT_MASK_MARKED);
-	    continue;
-	}
+	pmSourceAddModel (source->pixels, source->mask, model, false, false);
+	source->mode &= ~PM_SOURCE_SUBTRACTED;
     }
     psLogMsg ("psphot.replace", 4, "replace models: %f (%d objects)\n", psTimerMark ("psphot"), sources->n);
Index: /trunk/psphot/src/psphotSetup.c
===================================================================
--- /trunk/psphot/src/psphotSetup.c	(revision 5801)
+++ /trunk/psphot/src/psphotSetup.c	(revision 5802)
@@ -55,5 +55,5 @@
 	for (int iy = 0; iy < image->numRows; iy++) {
 	    for (int ix = 0; ix < image->numCols; ix++) {
-		weight->data.F32[iy][ix] = image->data.F32[iy][ix] / GAIN + PS_SQR(RDNOISE/GAIN);
+		weight->data.F32[iy][ix] = PS_MAX (image->data.F32[iy][ix] / GAIN + PS_SQR(RDNOISE/GAIN), 0.0);
 	    }
 	}
@@ -104,4 +104,7 @@
     psLogMsg ("psphot", 4, "load data: %f sec\n", psTimerMark ("psphot"));
 
+    psphotSaveImage (NULL, weight, "weight.fits");
+    psphotSaveImage (NULL, mask, "mask.fits");
+
     // save the image data & return it
     eamReadout *imdata = eamReadoutAlloc(image, weight, mask, header);
Index: /trunk/psphot/src/psphotSourceStats.c
===================================================================
--- /trunk/psphot/src/psphotSourceStats.c	(revision 5801)
+++ /trunk/psphot/src/psphotSourceStats.c	(revision 5802)
@@ -5,4 +5,5 @@
     bool     status  = false;
     psArray *sources = NULL;
+    float BIG_RADIUS;
 
     psTimerStart ("psphot");
@@ -35,14 +36,26 @@
 
 	// measure basic source moments
-	// XXX EAM : choose between these two versions
 	status = pmSourceMoments (source, RADIUS);
-	if (!status) {
-	  psFree (source);
-	  continue;
+	if (status) {
+	    // add to the source array
+	    psArrayAdd (sources, 100, source);
+	    psFree (source);
+	    continue;
 	}
-	
-	// add to the source array
-	psArrayAdd (sources, 100, source);
+
+	// if no valid pixels, or massive swing, likely saturated source,
+	// try a much larger box
+	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);
+	if (status) {
+	    // add to the source array
+	    psArrayAdd (sources, 100, source);
+	    psFree (source);
+	    continue;
+	}
+
 	psFree (source);
+	continue;
     }