Index: /trunk/psphot/src/Makefile.am
===================================================================
--- /trunk/psphot/src/Makefile.am	(revision 6480)
+++ /trunk/psphot/src/Makefile.am	(revision 6481)
@@ -37,4 +37,5 @@
 	psphotRadiusChecks.c	\
 	psphotImageMedian.c	\
+	psphotSkyReplace.c	\
 	psLine.c		\
 	psSparse.c		\
Index: /trunk/psphot/src/pmModelFitSet.c
===================================================================
--- /trunk/psphot/src/pmModelFitSet.c	(revision 6480)
+++ /trunk/psphot/src/pmModelFitSet.c	(revision 6481)
@@ -42,5 +42,5 @@
     model = params->data.F32[0];
 
-    PAR[0] = 0;
+    PAR[0] = params->data.F32[0];
     for (int i = 0; i < nSrc; i++) {
 	for (int n = 1; n < nPar; n++) {
@@ -53,4 +53,5 @@
 	}
     }
+    deriv->data.F32[0] = dPAR[0]*2.0;
     return (model);
 }
Index: /trunk/psphot/src/pmSourceFitSet.c
===================================================================
--- /trunk/psphot/src/pmSourceFitSet.c	(revision 6480)
+++ /trunk/psphot/src/pmSourceFitSet.c	(revision 6481)
@@ -74,10 +74,30 @@
     int nPar = model->params->n - 1;  // number of object parameters (excluding sky)
 
+    // define parameter vectors for source set
     psVector *params = psVectorAlloc (nSrc*nPar + 1, PS_TYPE_F32);
     psVector *dparams = psVectorAlloc (nSrc*nPar + 1, PS_TYPE_F32);
-    psVector *paramMask = PSF ? psVectorAlloc (nSrc*nPar + 1, PS_TYPE_U8) : NULL;
+
+    pmModelLimits modelLimits = pmModelLimits_GetFunction (model->type);
+
+    // define limits for single-source model
+    psVector *oneDelta;
+    psVector *oneParMin;
+    psVector *oneParMax;
+    modelLimits (&oneDelta, &oneParMin, &oneParMax);
+
+    psMinConstrain *constrain = psMinConstrainAlloc();
+    constrain->paramDelta = psVectorAlloc (nSrc*nPar + 1, PS_TYPE_F32);
+    constrain->paramMin = psVectorAlloc (nSrc*nPar + 1, PS_TYPE_F32);
+    constrain->paramMax = psVectorAlloc (nSrc*nPar + 1, PS_TYPE_F32);
+    constrain->paramMask = PSF ? psVectorAlloc (nSrc*nPar + 1, PS_TYPE_U8) : NULL;
 
     // all but the sky are allowed to vary independently (subject to PSF)
+    // set the parameters from the multiple models
+    // also, set limits based on single-source limits
     params->data.F32[0] = model->params->data.F32[0];
+    constrain->paramDelta->data.F32[0] = oneDelta->data.F32[0];
+    constrain->paramMin->data.F32[0]   = oneParMin->data.F32[0];
+    constrain->paramMax->data.F32[0]   = oneParMax->data.F32[0];
+    constrain->paramMask->data.U8[0]   = 0;
     for (int i = 0; i < nSrc; i++) {
       model = modelSet->data[i];
@@ -85,8 +105,20 @@
 	params->data.F32[i*nPar + n] = model->params->data.F32[n];
 	dparams->data.F32[i*nPar + n] = model->dparams->data.F32[n];
+	constrain->paramDelta->data.F32[i*nPar + n] = oneDelta->data.F32[n];
+	constrain->paramMin->data.F32[i*nPar + n]   = oneParMin->data.F32[n];
+	constrain->paramMax->data.F32[i*nPar + n]   = oneParMax->data.F32[n];
 	if (PSF) {
-	    paramMask->data.U8[i*nPar + n] = (n < 4) ? 0 : 1;
+	    constrain->paramMask->data.U8[i*nPar + n] = (n < 4) ? 0 : 1;
 	}
       }
+    }
+    psFree (oneDelta);
+    psFree (oneParMin);
+    psFree (oneParMax);
+
+    if (psTraceGetLevel(__func__) >= 5) {
+	for (int i = 0; i < params->n; i++) {
+	    fprintf (stderr, "%d %f %d\n", i, params->data.F32[i], constrain->paramMask->data.U8[i]);
+	}
     }
 
@@ -126,6 +158,6 @@
     }
     if (nPix <  nParams + 1) {
-        psTrace (".pmObjects.pmSourceFitModel", 4, "insufficient valid pixels\n");
-        psTrace(__func__, 3, "---- %s(false) end ----\n", __func__);
+        psTrace (__func__, 4, "insufficient valid pixels\n");
+        psTrace(__func__, 3, "---- %s() end : fail pixels ----\n", __func__);
         model->status = PM_MODEL_BADARGS;
 	psFree (x);
@@ -134,5 +166,9 @@
 	psFree (params);
 	psFree (dparams);
-	psFree (paramMask);
+	psFree(constrain->paramMask);
+	psFree(constrain->paramMin);
+	psFree(constrain->paramMax);
+	psFree(constrain->paramDelta);
+	psFree(constrain);
         return(false);
     }
@@ -143,19 +179,13 @@
     psMinimization *myMin = psMinimizationAlloc(PM_SOURCE_FIT_MODEL_NUM_ITERATIONS,
                             PM_SOURCE_FIT_MODEL_TOLERANCE);
-    psMinConstrain *constrain = psMinConstrainAlloc();
-    constrain->paramMask = paramMask;
-
-    // Set the parameter range checks
-    pmModelLimits modelLimits = pmModelLimits_GetFunction (model->type);
-    modelLimits (&constrain->paramDelta, &constrain->paramMin, &constrain->paramMax);
 
     psImage *covar = psImageAlloc (params->n, params->n, PS_TYPE_F64);
 
-    psTrace (".pmObjects.pmSourceFitSet", 5, "fitting function\n");
+    psTrace (__func__, 5, "fitting function\n");
     fitStatus = psMinimizeLMChi2(myMin, covar, params, constrain, x, y, yErr, pmModelFitSet);
 
     // parameter errors from the covariance matrix
     for (int i = 0; i < dparams->n; i++) {
-        if ((paramMask != NULL) && paramMask->data.U8[i])
+        if ((constrain->paramMask != NULL) && constrain->paramMask->data.U8[i])
             continue;
         dparams->data.F32[i] = sqrt(covar->data.F64[i][i]);
@@ -163,9 +193,9 @@
 
     // get the Gauss-Newton distance for fixed model parameters
-    if (paramMask != NULL) {
+    if (constrain->paramMask != NULL) {
         psVector *delta = psVectorAlloc (params->n, PS_TYPE_F64);
         psMinimizeGaussNewtonDelta(delta, params, NULL, x, y, yErr, pmModelFitSet);
         for (int i = 0; i < dparams->n; i++) {
-            if (!paramMask->data.U8[i])
+            if (!constrain->paramMask->data.U8[i])
                 continue;
             dparams->data.F32[i] = delta->data.F64[i];
@@ -220,5 +250,6 @@
 
     rc = (onPic && fitStatus);
-    psTrace(__func__, 3, "---- %s(%d) end ----\n", __func__, rc);
+    psTrace (__func__, 5, "onPic: %d, fitStatus: %d, nIter: %d, chisq: %f, nDof: %d\n", onPic, fitStatus, model->nIter, model->chisq, model->nDOF);
+    psTrace(__func__, 3, "---- %s end : status %d ----\n", __func__, rc);
     return(rc);
 }
Index: /trunk/psphot/src/psModulesUtils.c
===================================================================
--- /trunk/psphot/src/psModulesUtils.c	(revision 6480)
+++ /trunk/psphot/src/psModulesUtils.c	(revision 6481)
@@ -144,8 +144,10 @@
     int xIs, xJs, yIs, yJs;
     int xIe, yIe;
-    float flux;
+    float flux, wt;
 
     psImage *Pi = Mi->pixels;
     psImage *Pj = Mj->pixels;
+
+    psImage *Wi = Mi->weight;
 
     psImage *Ti = Mi->mask;
@@ -166,4 +168,6 @@
     yIe = Ye - Pi->row0;
 
+    // note that this is addressing the same image pixels,
+    // though only if both are source not model images
     flux = 0;
     for (yi = yIs, yj = yJs; yi < yIe; yi++, yj++) {
@@ -171,5 +175,54 @@
 	    if (Ti->data.U8[yi][xi]) continue;
 	    if (Tj->data.U8[yj][xj]) continue;
-	    flux += Pi->data.F32[yi][xi] * Pj->data.F32[yj][xj];
+	    wt = Wi->data.F32[yi][xi];
+	    if (wt > 0) {
+		flux += (Pi->data.F32[yi][xi] * Pj->data.F32[yj][xj]) / wt;
+	    } 
+	}
+    }
+    return (flux);
+}
+
+float pmSourceCrossWeight (pmSource *Mi, pmSource *Mj) {
+
+    int Xs, Xe, Ys, Ye;
+    int xi, xj, yi, yj;
+    int xIs, xJs, yIs, yJs;
+    int xIe, yIe;
+    float flux, wt;
+
+    psImage *Pi = Mi->pixels;
+    psImage *Pj = Mj->pixels;
+
+    psImage *Wi = Mi->weight;
+
+    psImage *Ti = Mi->mask;
+    psImage *Tj = Mj->mask;
+
+    Xs = PS_MAX (Pi->col0, Pj->col0);
+    Xe = PS_MIN (Pi->col0 + Pi->numCols, Pj->col0 + Pj->numCols);
+    
+    Ys = PS_MAX (Pi->row0, Pj->row0);
+    Ye = PS_MIN (Pi->row0 + Pi->numRows, Pj->row0 + Pj->numRows);
+    
+    xIs = Xs - Pi->col0;
+    xJs = Xs - Pj->col0;
+    yIs = Ys - Pi->row0;
+    yJs = Ys - Pj->row0;
+
+    xIe = Xe - Pi->col0;
+    yIe = Ye - Pi->row0;
+
+    // note that this is addressing the same image pixels,
+    // though only if both are source not model images
+    flux = 0;
+    for (yi = yIs, yj = yJs; yi < yIe; yi++, yj++) {
+	for (xi = xIs, xj = xJs; xi < xIe; xi++, xj++) {
+	    if (Ti->data.U8[yi][xi]) continue;
+	    if (Tj->data.U8[yj][xj]) continue;
+	    wt = Wi->data.F32[yi][xi];
+	    if (wt > 0) {
+		flux += 1.0 / wt;
+	    } 
 	}
     }
@@ -250,2 +303,11 @@
     return NULL;
 }
+
+bool psRegionIsNaN (psRegion region) {
+
+    if (!isfinite(region.x0)) return true;
+    if (!isfinite(region.x1)) return true;
+    if (!isfinite(region.y0)) return true;
+    if (!isfinite(region.y1)) return true;
+    return false;
+}
Index: /trunk/psphot/src/psModulesUtils.h
===================================================================
--- /trunk/psphot/src/psModulesUtils.h	(revision 6480)
+++ /trunk/psphot/src/psModulesUtils.h	(revision 6481)
@@ -44,9 +44,11 @@
 // unify with paul's image/header/metadata functions
 psF32        pmConfigLookupF32 (bool *status, psMetadata *config, psMetadata *header, char *name);
-pmModel *pmModelSelect (pmSource *source);
+pmModel     *pmModelSelect (pmSource *source);
 
-ppConfig       *ppConfigAlloc (void);
-ppFile         *ppFileAlloc (void);
-psMetadata    *pmReadoutGetHeader (pmReadout *readout);
+ppConfig    *ppConfigAlloc (void);
+ppFile      *ppFileAlloc (void);
+psMetadata  *pmReadoutGetHeader (pmReadout *readout);
+
+bool psRegionIsNaN (psRegion region);
 
 # endif
Index: /trunk/psphot/src/psSparse.c
===================================================================
--- /trunk/psphot/src/psSparse.c	(revision 6480)
+++ /trunk/psphot/src/psSparse.c	(revision 6481)
@@ -34,17 +34,22 @@
     sparse->Bfj->data.F32[2] = B->data.F32[2];
 
-    x = psSparseSolve (x, sparse, 0);
-    fprintf (stderr, "x: %f %f %f\n", x->data.F32[0], x->data.F32[1], x->data.F32[2]);
-
-    x = psSparseSolve (x, sparse, 1);
-    fprintf (stderr, "x: %f %f %f\n", x->data.F32[0], x->data.F32[1], x->data.F32[2]);
-
-    x = psSparseSolve (x, sparse, 2);
-    fprintf (stderr, "x: %f %f %f\n", x->data.F32[0], x->data.F32[1], x->data.F32[2]);
-
-    x = psSparseSolve (x, sparse, 3);
-    fprintf (stderr, "x: %f %f %f\n", x->data.F32[0], x->data.F32[1], x->data.F32[2]);
-
-    x = psSparseSolve (x, sparse, 4);
+    psSparseConstraint constraint;
+    constraint.paramMin   = -1e8;
+    constraint.paramMax   = +1e8;
+    constraint.paramDelta = +1e8;
+
+    x = psSparseSolve (x, constraint, sparse, 0);
+    fprintf (stderr, "x: %f %f %f\n", x->data.F32[0], x->data.F32[1], x->data.F32[2]);
+
+    x = psSparseSolve (x, constraint, sparse, 1);
+    fprintf (stderr, "x: %f %f %f\n", x->data.F32[0], x->data.F32[1], x->data.F32[2]);
+
+    x = psSparseSolve (x, constraint, sparse, 2);
+    fprintf (stderr, "x: %f %f %f\n", x->data.F32[0], x->data.F32[1], x->data.F32[2]);
+
+    x = psSparseSolve (x, constraint, sparse, 3);
+    fprintf (stderr, "x: %f %f %f\n", x->data.F32[0], x->data.F32[1], x->data.F32[2]);
+
+    x = psSparseSolve (x, constraint, sparse, 4);
     fprintf (stderr, "x: %f %f %f\n", x->data.F32[0], x->data.F32[1], x->data.F32[2]);
     return;
@@ -166,5 +171,5 @@
 }
 
-psVector *psSparseSolve (psVector *guess, psSparse *sparse, int Niter) {
+psVector *psSparseSolve (psVector *guess, psSparseConstraint constraint, psSparse *sparse, int Niter) {
 
     psF32 dG;
@@ -182,5 +187,14 @@
 	for (int i = 0; i < dQ->n; i++) {
 	    dG = (dQ->data.F32[i] - Bfj->data.F32[i]) / Qii->data.F32[i];
+	    if (fabs (dG) > constraint.paramDelta) {
+		if (dG > 0) {
+		    dG = +constraint.paramDelta;
+		} else {
+		    dG = -constraint.paramDelta;
+		}
+	    }
 	    guess->data.F32[i] -= dG;
+	    guess->data.F32[i] = PS_MAX (guess->data.F32[i], constraint.paramMin);
+	    guess->data.F32[i] = PS_MIN (guess->data.F32[i], constraint.paramMax);
 	}
     }
@@ -220,3 +234,2 @@
 }
 
-
Index: /trunk/psphot/src/psSparse.h
===================================================================
--- /trunk/psphot/src/psSparse.h	(revision 6480)
+++ /trunk/psphot/src/psSparse.h	(revision 6481)
@@ -1,2 +1,8 @@
+
+typedef struct {
+    double paramDelta;
+    double paramMin;
+    double paramMax;
+} psSparseConstraint;
 
 typedef struct {
@@ -15,4 +21,4 @@
 void psSparseVectorElement (psSparse *sparse, int i, float value);
 psVector *psSparseMatrixTimesVector (psVector *output, psSparse *matrix, psVector *vector);
-psVector *psSparseSolve (psVector *guess, psSparse *sparse, int Niter);
+psVector *psSparseSolve (psVector *guess, psSparseConstraint constraint, psSparse *sparse, int Niter);
 psSparse *psSparseAlloc (int Nrows, int Nelem);
Index: /trunk/psphot/src/psphot.c
===================================================================
--- /trunk/psphot/src/psphot.c	(revision 6480)
+++ /trunk/psphot/src/psphot.c	(revision 6481)
@@ -26,5 +26,5 @@
     psFree (config);
     psphotCleanup ();
-
+ 
     exit (0);
 }
Index: /trunk/psphot/src/psphot.h
===================================================================
--- /trunk/psphot/src/psphot.h	(revision 6480)
+++ /trunk/psphot/src/psphot.h	(revision 6481)
@@ -42,5 +42,5 @@
 
 // optional object analysis steps
-bool 	        psphotEnsemblePSF (pmReadout *readout, psMetadata *config, psArray *sources, pmPSF *psf);
+bool 	        psphotEnsemblePSF (pmReadout *readout, psMetadata *config, psArray *sources, pmPSF *psf, bool final);
 bool            psphotBlendFit (pmReadout *readout, psMetadata *config, psArray *sources, pmPSF *psf);
 bool            psphotReplaceUnfit (psArray *sources);
@@ -56,4 +56,5 @@
 pmModel        *pmSourceMagnitudes (pmSource *source, pmPSF *psf, float apRadius);
 float           pmSourceCrossProduct (pmSource *Mi, pmSource *Mj);
+float           pmSourceCrossWeight (pmSource *Mi, pmSource *Mj);
 psArray        *pmSourceContour_EAM (psImage *image, int x, int y, float threshold);
 void            psphotModelGroupInit (void);
@@ -63,8 +64,11 @@
 void            psphotTestArguments (int *argc, char **argv);
 bool            pmCellSetMask (pmCell *cell, psMetadata *recipe);
+bool            psphotBackgroundNames (psMetadata *arguments);
+bool            psphotSkyReplace (pmReadout *readout, psImage *background);
 
 // functions to set the correct source pixels
 bool 		psphotInitRadiusPSF (psMetadata *config, pmModelType type);
 bool 		psphotCheckRadiusPSF (pmReadout *readout, pmSource *source, pmModel *model);
+bool 		psphotCheckRadiusPSFBlend (pmReadout *readout, pmSource *source, pmModel *model, float dR);
 bool 		psphotInitRadiusEXT (psMetadata *config, pmModelType type);
 bool 		psphotCheckRadiusEXT (pmReadout *readout, pmSource *source, pmModel *model);
Index: /trunk/psphot/src/psphotArguments.c
===================================================================
--- /trunk/psphot/src/psphotArguments.c	(revision 6480)
+++ /trunk/psphot/src/psphotArguments.c	(revision 6481)
@@ -42,4 +42,18 @@
 	psArgumentRemove (N, argc, argv);
 	psMetadataAddStr (config->arguments, PS_LIST_TAIL, "RESID_IMAGE", PS_META_REPLACE, "", argv[N]);
+	psArgumentRemove (N, argc, argv);
+    }
+
+    // background image (output) is used by psphotOutput
+    if ((N = psArgumentGet (*argc, argv, "-background"))) {
+	psArgumentRemove (N, argc, argv);
+	psMetadataAddStr (config->arguments, PS_LIST_TAIL, "BACKGROUND_IMAGE", PS_META_REPLACE, "", argv[N]);
+	psArgumentRemove (N, argc, argv);
+    }
+
+    // background model (output) is used by psphotOutput
+    if ((N = psArgumentGet (*argc, argv, "-backmodel"))) {
+	psArgumentRemove (N, argc, argv);
+	psMetadataAddStr (config->arguments, PS_LIST_TAIL, "BACKGROUND_MODEL", PS_META_REPLACE, "", argv[N]);
 	psArgumentRemove (N, argc, argv);
     }
@@ -89,4 +103,11 @@
 	psArgumentRemove (N, argc, argv);
 	psMetadataAddStr (config->options, PS_LIST_TAIL, "PSF_INPUT_FILE", PS_META_REPLACE, "", argv[N]);
+	psArgumentRemove (N, argc, argv);
+    }
+
+    // PSF : optional psf file to be loaded
+    if ((N = psArgumentGet (*argc, argv, "-psfout"))) {
+	psArgumentRemove (N, argc, argv);
+	psMetadataAddStr (config->options, PS_LIST_TAIL, "PSF_OUTPUT_FILE", PS_META_REPLACE, "", argv[N]);
 	psArgumentRemove (N, argc, argv);
     }
Index: /trunk/psphot/src/psphotBasicDeblend.c
===================================================================
--- /trunk/psphot/src/psphotBasicDeblend.c	(revision 6480)
+++ /trunk/psphot/src/psphotBasicDeblend.c	(revision 6481)
@@ -86,5 +86,5 @@
 	threshold = FRACTION * source->moments->Peak;
 	// threshold is no less than NSIGMA above the local median sigma?
-	threshold = PS_MAX (threshold, NSIGMA*source->moments->dSky);
+	threshold = PS_MAX (threshold, NSIGMA*sqrt(source->moments->dSky));
 
 	// generate a basic contour (set of x,y coordinates at-or-below flux level)
Index: /trunk/psphot/src/psphotBlendFit.c
===================================================================
--- /trunk/psphot/src/psphotBlendFit.c	(revision 6480)
+++ /trunk/psphot/src/psphotBlendFit.c	(revision 6481)
@@ -3,4 +3,6 @@
 // XXX I don't like this name
 bool psphotBlendFit (pmReadout *readout, psMetadata *config, psArray *sources, pmPSF *psf) { 
+
+    bool status;
 
     psTimerStart ("psphot");
@@ -9,10 +11,17 @@
     sources = psArraySort (sources, psphotSortBySN);
     
+    // S/N limit to perform full non-linear fits
+    float FIT_SN_LIM = psMetadataLookupF32 (&status, config, "FULL_FIT_SN_LIM");
+
     psphotInitLimitsPSF (config);
     psphotInitLimitsEXT (config);
     psphotInitRadiusPSF (config, psf->type);
 
+    // option to limit analysis to a specific region
+    char *region = psMetadataLookupStr (&status, config, "ANALYSIS_REGION");
+    psRegion AnalysisRegion = psRegionForImage (readout->image, psRegionFromString (region));
+    if (psRegionIsNaN (AnalysisRegion)) psAbort ("psphot", "analysis region mis-defined");
+
     for (int i = 0; i < sources->n; i++) {
-
 	pmSource *source = sources->data[i];
 
@@ -22,6 +31,15 @@
 	if (source->type == PM_SOURCE_SATURATED) continue;
 
+	// limit selection to some SN limit
+	if (source->moments == NULL) continue;
+	if (source->moments->SN < FIT_SN_LIM) continue;
+
 	// if model is NULL, we don't have a starting guess
 	if (source->modelPSF == NULL) continue;
+
+	if (source->moments->x < AnalysisRegion.x0) continue;
+	if (source->moments->y < AnalysisRegion.y0) continue;
+	if (source->moments->x > AnalysisRegion.x1) continue;
+	if (source->moments->y > AnalysisRegion.y1) continue;
 
 	// replace object in image
@@ -29,5 +47,5 @@
 
 	psTrace ("psphot.blend", 5, "trying source at %f, %f\n", source->moments->x, source->moments->y);
-
+	
 	// try fitting PSFs, then try extended sources
 	if (psphotFitBlend (readout, source)) continue;
Index: /trunk/psphot/src/psphotChoosePSF.c
===================================================================
--- /trunk/psphot/src/psphotChoosePSF.c	(revision 6480)
+++ /trunk/psphot/src/psphotChoosePSF.c	(revision 6481)
@@ -11,4 +11,13 @@
 
     psTimerStart ("psphot");
+
+    // check if a PSF model is supplied by the user
+    char *psfFile = psMetadataLookupStr (&status, config, "PSF_INPUT_FILE");
+    if (status) {
+	int Nfail;
+	psMetadata *psfData = psMetadataConfigParse (NULL, &Nfail, psfFile, FALSE);
+	pmPSF *psf = pmPSFfromMD (psfData);
+	return psf;
+    }
 
     // examine PSF sources in S/N order (brightest first)
Index: /trunk/psphot/src/psphotEnsemblePSF.c
===================================================================
--- /trunk/psphot/src/psphotEnsemblePSF.c	(revision 6480)
+++ /trunk/psphot/src/psphotEnsemblePSF.c	(revision 6481)
@@ -1,7 +1,8 @@
 # include "psphot.h"
+bool pmSourceChisq (pmModel *model, psImage *image, psImage *mask, psImage *weight);
 
 // 2006.02.07 : no leaks!
 // fit all reasonable sources with the linear PSF model
-bool psphotEnsemblePSF (pmReadout *readout, psMetadata *config, psArray *sources, pmPSF *psf) { 
+bool psphotEnsemblePSF (pmReadout *readout, psMetadata *config, psArray *sources, pmPSF *psf, bool final) { 
 
     bool  status;
@@ -9,4 +10,7 @@
     float y;
     float f;
+    float r;
+
+    // psRegion allArray = psRegionSet (0, 0, 0, 0);
 
     psTimerStart ("psphot");
@@ -29,14 +33,9 @@
 
     // option to limit analysis to a specific region
-    bool UseAnalysisRegion = false;
-    psRegion AnalysisRegion = {0, 0, 0, 0};
     char *region = psMetadataLookupStr (&status, config, "ANALYSIS_REGION");
-    if (status) {
-	UseAnalysisRegion = true;
-	AnalysisRegion = psRegionFromString (region);
-	psLogMsg ("psphotEnsemblePSF", 4, "using region %f,%f - %f,%f\n", 
-		  AnalysisRegion.x0, AnalysisRegion.y0, 
-		  AnalysisRegion.x1, AnalysisRegion.y1);
-    }
+    psRegion AnalysisRegion = psRegionFromString (region);
+    AnalysisRegion = psRegionForImage (readout->image, AnalysisRegion);
+    // psRegion AnalysisRegion = psRegionForImage (readout->image, psRegionFromString (region));
+    if (psRegionIsNaN (AnalysisRegion)) psAbort ("psphot", "analysis region mis-defined");
 
     for (int i = 0; i < sources->n; i++) {
@@ -45,14 +44,16 @@
 	// skip non-astronomical objects (very likely defects)
 	// XXX EAM : should we try these anyway?
-	if (inSource->mode &  PM_SOURCE_BLEND) continue;
 	if (inSource->type == PM_SOURCE_DEFECT) continue; 
 	if (inSource->type == PM_SOURCE_SATURATED) continue;
-
-	if (UseAnalysisRegion) {
-	    if (inSource->moments->x < AnalysisRegion.x0) continue;
-	    if (inSource->moments->y < AnalysisRegion.y0) continue;
-	    if (inSource->moments->x > AnalysisRegion.x1) continue;
-	    if (inSource->moments->y > AnalysisRegion.y1) continue;
-	}
+	if (final) {
+	    if (inSource->mode &  PM_SOURCE_SUBTRACTED) continue;
+	} else {
+	    if (inSource->mode &  PM_SOURCE_BLEND) continue;
+	}
+
+	if (inSource->moments->x < AnalysisRegion.x0) continue;
+	if (inSource->moments->y < AnalysisRegion.y0) continue;
+	if (inSource->moments->x > AnalysisRegion.x1) continue;
+	if (inSource->moments->y > AnalysisRegion.y1) continue;
 
 	pmSource *otSource = pmSourceAlloc ();
@@ -73,6 +74,7 @@
 	    // peak-up on peak (for non-sat objects)
 
-	    int ix = inSource->peak->x;
-	    int iy = inSource->peak->y;
+	    // ix,iy must land on inSource->pixels 
+	    int ix = PS_MAX (inSource->pixels->col0 + 1, PS_MIN (inSource->pixels->col0 + inSource->pixels->numCols - 2, inSource->peak->x));
+	    int iy = PS_MAX (inSource->pixels->row0 + 1, PS_MIN (inSource->pixels->row0 + inSource->pixels->numRows - 2, inSource->peak->y));
 
 	    psPolynomial2D *bicube = psImageBicubeFit (inSource->pixels, ix, iy);
@@ -104,4 +106,5 @@
 	otSource->mask   = psImageCopy (NULL, inSource->mask,   PS_TYPE_U8);
 	otSource->pixels = psImageCopy (NULL, inSource->pixels, PS_TYPE_F32);
+	otSource->weight = psImageCopy (NULL, inSource->weight, PS_TYPE_F32);
 
 	// build the model image 
@@ -117,4 +120,8 @@
 	psImageKeepCircle (mask, x, y, model->radius, "OR", PSPHOT_MASK_MARKED);
 	pmSourceAddModel (flux, mask, model, false, false);
+
+	// calculate nDOF (nPix - 1)
+	// int Nmaskpix = psImageCountPixelMask (mask, allArray, PSPHOT_MASK_SATURATED);
+	// model->nDOF  = mask->numCols*mask->numRows - Nmaskpix - 1;
 
 	// save source in list
@@ -128,4 +135,5 @@
     // fill out the sparse matrix
     psSparse *sparse = psSparseAlloc (models->n, 100);
+    psVector *weight = psVectorAlloc (models->n, PS_TYPE_F32);
     for (int i = 0; i < models->n; i++) {
 	int N = index->data.U32[i];
@@ -133,11 +141,13 @@
 	pmSource *Mi = models->data[i];
 
-	// diagonal element (auto-cross-product)
-	f = pmSourceCrossProduct (Mi, Mi);
-	psSparseMatrixElement (sparse, i, i, f);
+	// scale by diagonal element (auto-cross-product)
+	r = pmSourceCrossProduct (Mi, Mi);
+	weight->data.F32[i] = r;
+
+	psSparseMatrixElement (sparse, i, i, 1.0);
 
 	// find the image x model value
 	f = pmSourceCrossProduct (Fi, Mi);
-	psSparseVectorElement (sparse, i, f);
+	psSparseVectorElement (sparse, i, f / r);
 
 	// loop over all other stars following this one
@@ -153,12 +163,17 @@
 	    // got an overlap; calculate cross-product and add to output array
 	    f = pmSourceCrossProduct (Mi, Mj);
-	    psSparseMatrixElement (sparse, j, i, f); 
+	    psSparseMatrixElement (sparse, j, i, f / r); 
 	}
     }
     psLogMsg ("psphot.emsemble", 4, "built matrix: %f (%d elements)\n", psTimerMark ("psphot"), sparse->Nelem);
+
+    psSparseConstraint constraint;
+    constraint.paramMin   = 0;
+    constraint.paramMax   = 1e8;
+    constraint.paramDelta = 1e8;
 
     // solve for normalization terms (need include local sky?)
     psSparseResort (sparse);
-    psVector *norm = psSparseSolve (NULL, sparse, 3);
+    psVector *norm = psSparseSolve (NULL, constraint, sparse, 3);
 
     // adjust models, set sources and subtract
@@ -168,4 +183,7 @@
 	pmSource *Mi = models->data[i];
 
+	// if we already have a PSF model, free it.
+	psFree (Fi->modelPSF);
+
 	// need to increment counter so we can free models here and sources above
 	Fi->modelPSF = psMemCopy (Mi->modelPSF);
@@ -173,16 +191,52 @@
 	// assign linearly-fitted normalization
 	Fi->modelPSF->params->data.F32[1] = norm->data.F32[i];
+	Fi->modelPSF->dparams->data.F32[1] = sqrt(sqrt(2) * norm->data.F32[i] / (sparse->Bfj->data.F32[i] * weight->data.F32[i]));
+	// XXX EAM : this factor of sqrt(2) makes the errors consistent, but I don't understand it
 
 	// subtract object
 	pmSourceSubModel (Fi->pixels, Fi->mask, Fi->modelPSF, false, false);
 	Fi->mode |= PM_SOURCE_SUBTRACTED;
-	Fi->mode |= PM_SOURCE_TEMPSUB;
-    }
+	if (!final) Fi->mode |= PM_SOURCE_TEMPSUB;
+    }
+
+    // measure chisq for each source
+    for (int i = 0; final && (i < models->n); i++) {
+	int N = index->data.U32[i];
+	pmSource *Fi = sources->data[N];
+	pmModel *model = Fi->modelPSF;
+
+	x = model->params->data.F32[2];
+	y = model->params->data.F32[3];
+
+	psImageKeepCircle (Fi->mask, x, y, model->radius, "OR", PSPHOT_MASK_MARKED);
+	pmSourceChisq (model, Fi->pixels, Fi->mask, Fi->weight);
+	psImageKeepCircle (Fi->mask, x, y, model->radius, "AND", ~PSPHOT_MASK_MARKED);
+    }
+
     psFree (index);
     psFree (sparse);
     psFree (models);
     psFree (norm);
+    psFree (weight);
 
     psLogMsg ("psphot.emsemble", 3, "measure ensemble of PSFs: %f sec\n", psTimerMark ("psphot"));
     return true;
 }
+
+bool pmSourceChisq (pmModel *model, psImage *image, psImage *mask, psImage *weight) {
+
+    double dC = 0.0;
+    int Npix = 0;
+    for (int j = 0; j < image->numRows; j++) {
+	for (int i = 0; i < image->numCols; i++) {
+	    if (mask->data.U8[j][i]) continue;
+	    if (weight->data.F32[j][i] <= 0) continue;
+	    dC += PS_SQR (image->data.F32[j][i]) / weight->data.F32[j][i];
+	    Npix ++;
+	}
+    }
+    model->nDOF = Npix - 1;
+    model->chisq = dC;
+
+    return (true);
+}
Index: /trunk/psphot/src/psphotFindPeaks.c
===================================================================
--- /trunk/psphot/src/psphotFindPeaks.c	(revision 6480)
+++ /trunk/psphot/src/psphotFindPeaks.c	(revision 6481)
@@ -10,7 +10,4 @@
 
     // smooth the image and weight map
-
-    psphotSaveImage (NULL, readout->image, "image.fits");
-
     psTimerStart ("psphot");
 
@@ -22,11 +19,7 @@
     psLogMsg ("psphot", 4, "smooth image: %f sec\n", psTimerMark ("psphot"));
 
-    psphotSaveImage (NULL, smooth_im, "smooth_im.fits");
-
     psImage *smooth_wt = psImageCopy (NULL, readout->weight, PS_TYPE_F32);
     psImageSmooth (smooth_wt, SIGMA, NSIGMA);
     psLogMsg ("psphot", 4, "smooth weight: %f sec\n", psTimerMark ("psphot"));
-
-    psphotSaveImage (NULL, smooth_wt, "smooth_wt.fits");
 
     psTimerStart ("psphot");
Index: /trunk/psphot/src/psphotImageMedian.c
===================================================================
--- /trunk/psphot/src/psphotImageMedian.c	(revision 6480)
+++ /trunk/psphot/src/psphotImageMedian.c	(revision 6481)
@@ -8,4 +8,17 @@
 // use no more than MAX_SAMPLE_PIXELS pixels for each median box
 static int MAX_SAMPLE_PIXELS;
+
+static char *backImage = NULL;
+static char *backModel = NULL;
+
+bool psphotBackgroundNames (psMetadata *arguments) {
+
+    bool status;
+
+    backImage = psMetadataLookupStr (&status, arguments, "BACKGROUND_IMAGE");
+    backModel = psMetadataLookupStr (&status, arguments, "BACKGROUND_MODEL");
+
+    return true;
+}
 
 // generate the median in NxN boxes, clipping heavily
@@ -211,4 +224,9 @@
 	}
     }
+
+    // optionally save background
+    if (backImage != NULL) psphotSaveImage (NULL, background, backImage);
+    if (backModel != NULL) psphotSaveImage (NULL, model, backModel);
+
     psLogMsg ("psphot", 3, "subtracted background model: %f sec\n", psTimerMark ("psphot"));
     psFree (rnd);
Index: /trunk/psphot/src/psphotModelTest.c
===================================================================
--- /trunk/psphot/src/psphotModelTest.c	(revision 6480)
+++ /trunk/psphot/src/psphotModelTest.c	(revision 6481)
@@ -26,5 +26,5 @@
     // in psfMode, psf sets the model type
     if (psfMode) {
-	char *psfFile = psMetadataLookupStr (&status, arguments, "PSF_INPUT_FILE");
+	char *psfFile = psMetadataLookupStr (&status, recipe, "PSF_INPUT_FILE");
 	if (!status) psAbort ("psphotModelTest", "PSF_INPUT_FILE not supplied");
 	psMetadata *psfData = psMetadataConfigParse (NULL, &Nfail, psfFile, FALSE);
@@ -157,5 +157,5 @@
     if (status) {
 	status = psphotFitSet (source, model, fitset, psfMode);
-	return status;
+	exit (0);
     }
 
Index: /trunk/psphot/src/psphotOutput.c
===================================================================
--- /trunk/psphot/src/psphotOutput.c	(revision 6480)
+++ /trunk/psphot/src/psphotOutput.c	(revision 6481)
@@ -11,4 +11,6 @@
 static char *extNameKey = NULL;
 static char *extRoot    = NULL;
+static char *psfFile    = NULL;
+static char *psfSample  = NULL;
 
 void psphotOutputCleanup () {
@@ -21,4 +23,6 @@
     psFree (extNameKey);
     psFree (extRoot);
+    psFree (psfFile);
+    psFree (psfSample);
     return;
 }
@@ -34,4 +38,7 @@
     outputMode      = psMetadataLookupStr (&status, config->recipe, "OUTPUT_MODE");
     outputFormat    = psMetadataLookupStr (&status, config->recipe, "OUTPUT_FORMAT");
+
+    psfFile         = psMetadataLookupStr (&status, config->recipe, "PSF_OUTPUT_FILE");
+    psfSample       = psMetadataLookupStr (&status, config->recipe, "PSF_SAMPLE_FILE");
 
     // rules to construct output names
@@ -75,4 +82,7 @@
     }
 
+    // register background image-file names
+    psphotBackgroundNames (config->arguments);
+
     // save so freeing config will not drop these references
     psMemCopy (outputRoot);
@@ -83,4 +93,6 @@
     psMemCopy (extNameKey);
     psMemCopy (extRoot);
+    psMemCopy (psfFile);
+    psMemCopy (psfSample);
 
     return;
@@ -165,6 +177,4 @@
     char *outputFile;
 
-    psTimerStart ("psphot");
-
     psMetadata *header = pmReadoutGetHeader (readout);
 
@@ -180,6 +190,4 @@
     } 
 
-    char *psfFile    = psMetadataLookupStr (&status, arguments, "PSF_OUTPUT_FILE");
-    char *psfSample  = psMetadataLookupStr (&status, arguments, "PSF_SAMPLE_FILE");
     char *residImage = psMetadataLookupStr (&status, arguments, "RESID_IMAGE");
 
@@ -236,4 +244,6 @@
     psF32 *PAR, *dPAR;
     float dmag, apResid;
+
+    psTimerStart ("string");
 
     psLine *line = psLineAlloc (104);  // 104 is dophot-defined line length
@@ -274,4 +284,48 @@
     }
     fclose (f);
+    psFree (line);
+    fprintf (stderr, "%f seconds for %d objects with psLine\n", psTimerMark ("string"), (int)sources->n);
+
+    psTimerStart ("string");
+
+    f = fopen ("test.obj", "w");
+    if (f == NULL) {
+	psLogMsg ("WriteSourceOBJ", 3, "can't open output file for output %s\n", "test.obj");
+	return false;
+    }
+
+    char *string;
+    // write sources with models 
+    for (int i = 0; i < sources->n; i++) {
+	pmSource *source = (pmSource *) sources->data[i];
+	pmModel *model = pmModelSelect (source);
+	if (model == NULL) continue;
+
+	PAR = model->params->data.F32;
+	dPAR = model->dparams->data.F32;
+
+	dmag = dPAR[1] / PAR[1];
+	type = pmSourceDophotType (source);
+	apResid = source->apMag - source->fitMag;
+
+	string = NULL;
+	psStringAppend (&string, "%3d",   type);
+	psStringAppend (&string, "%8.2f", PAR[2]);
+	psStringAppend (&string, "%8.2f", PAR[3]);
+	psStringAppend (&string, "%8.3f", source->fitMag);
+	psStringAppend (&string, "%6.3f", dmag);
+	psStringAppend (&string, "%9.2f", PAR[0]);
+	psStringAppend (&string, "%9.3f", PAR[4]);
+	psStringAppend (&string, "%9.3f", PAR[5]);
+	psStringAppend (&string, "%7.2f", PAR[6]);
+	psStringAppend (&string, "%8.3f", 99.999);
+	psStringAppend (&string, "%8.3f", source->apMag);
+	psStringAppend (&string, "%8.2f\n", apResid);
+	fwrite (string, 1, strlen(string), f);
+	psFree (string);
+    }
+    fclose (f);
+    fprintf (stderr, "%f seconds for %d objects with psString\n", psTimerMark ("string"), (int)sources->n);
+
     return true;
 }
@@ -501,5 +555,5 @@
     }
 
-   // write sources with models first 
+    // write sources with models first 
     for (i = 0; i < sources->n; i++) {
 	pmSource *source = (pmSource *) sources->data[i];
@@ -525,5 +579,5 @@
 	    fprintf (f, "%9.6f ", dPAR[j]);
 	}
-	fprintf (f, ": %7.4f %2d %#5x %7.3f %7.1f %7.2f %4d %2d\n", 
+	fprintf (f, ": %8.4f %2d %#5x %7.3f %7.1f %7.2f %4d %2d\n", 
 		 source[0].apMag, source[0].type, source[0].mode, 
 		 log10(model[0].chisq/model[0].nDOF), 
@@ -617,7 +671,7 @@
 	    
 	if (source->moments == NULL) {
-	  moment = empty;
+	    moment = empty;
 	} else {
-	  moment = source->moments;
+	    moment = source->moments;
 	}
 
Index: /trunk/psphot/src/psphotRadiusChecks.c
===================================================================
--- /trunk/psphot/src/psphotRadiusChecks.c	(revision 6480)
+++ /trunk/psphot/src/psphotRadiusChecks.c	(revision 6481)
@@ -24,4 +24,21 @@
     // set the fit radius based on the object flux limit and the model
     model->radius = modelRadiusPSF (model->params, PSF_FIT_NSIGMA*moments->dSky) + PSF_FIT_PADDING;
+    if (isnan(model->radius)) psAbort ("apply_psf_model", "error in radius");
+	
+    if (source->mode &  PM_SOURCE_SATSTAR) {
+	model->radius *= 2;
+    }
+
+    bool status = psphotRedefinePixels (source, readout, model->params->data.F32[2], model->params->data.F32[3], model->radius);
+    return status;
+}
+
+bool psphotCheckRadiusPSFBlend (pmReadout *readout, pmSource *source, pmModel *model, float dR) {
+
+    pmMoments *moments = source->moments;
+    if (moments == NULL) return false;
+
+    // set the fit radius based on the object flux limit and the model
+    model->radius = modelRadiusPSF (model->params, PSF_FIT_NSIGMA*moments->dSky) + dR + PSF_FIT_PADDING;
     if (isnan(model->radius)) psAbort ("apply_psf_model", "error in radius");
 	
Index: /trunk/psphot/src/psphotReadout.c
===================================================================
--- /trunk/psphot/src/psphotReadout.c	(revision 6480)
+++ /trunk/psphot/src/psphotReadout.c	(revision 6481)
@@ -17,7 +17,9 @@
 
     // construct sources and measure basic stats
+    // limit moments analysis by S/N?
     sources = psphotSourceStats (readout, config, peaks);
 
     // classify sources based on moments, brightness
+    // faint sources not classified?
     psphotRoughClass (sources, config);
 
@@ -28,23 +30,32 @@
     psf = psphotChoosePSF (config, sources);
 
-    // select analysis method
-    char *FITMODE = psMetadataLookupStr (&status, config, "FITMODE");
-    if (!strcasecmp(FITMODE, "ENSEMBLE")) {
-	psphotEnsemblePSF (readout, config, sources, psf);
-    }
+    // linear PSF fit to peaks
+    psphotEnsemblePSF (readout, config, sources, psf, FALSE);
 
-    if (!strcasecmp(FITMODE, "BLEND")) {
-	psphotEnsemblePSF (readout, config, sources, psf);
-	psphotBlendFit (readout, config, sources, psf);
-	psphotReplaceUnfit (sources);
-	psphotApResid (readout, sources, config, psf);
-    }
+    // non-linear PSF and EXT fit to brighter sources
+    psphotBlendFit (readout, config, sources, psf);
 
+    // replace fitted sources
+    psphotReplaceUnfit (sources);
+
+    // find remaining peaks after first source subtraction pass
+    // psphotFindPeaks ();
+
+    // linear PSF fit to remaining peaks
+    psphotEnsemblePSF (readout, config, sources, psf, TRUE);
+
+    // measure aperture photometry corrections
+    psphotApResid (readout, sources, config, psf);
+
+    // calculate source magnitudes
     psphotMagnitudes (config, sources, psf);
 
+    // update the header with stats results
     psMetadata *header = pmReadoutGetHeader (readout);
     psphotUpdateHeader (header, config);
 
-    // psphotSkyReplace (readout, skymodel);
+    // XXX make this an option?
+    // replace background in residual image
+    psphotSkyReplace (readout, skymodel);
 
     // need to do something with the sources, psf, and sky
@@ -66,4 +77,10 @@
 // XXX Deprecate or allow these versions?
 # if (0)
+    // select analysis method
+    char *FITMODE = psMetadataLookupStr (&status, config, "FITMODE");
+    if (!strcasecmp(FITMODE, "ENSEMBLE")) {
+	psphotEnsemblePSF (readout, config, sources, psf);
+    }
+
     if (!strcasecmp(FITMODE, "FULL")) {
 	psphotEnsemblePSF (readout, config, sources, psf);
Index: /trunk/psphot/src/psphotSkyReplace.c
===================================================================
--- /trunk/psphot/src/psphotSkyReplace.c	(revision 6481)
+++ /trunk/psphot/src/psphotSkyReplace.c	(revision 6481)
@@ -0,0 +1,21 @@
+# include "psphot.h"
+
+// generate the median in NxN boxes, clipping heavily
+// linear interpolation to generate full-scale model
+bool psphotSkyReplace (pmReadout *readout, psImage *background) { 
+
+    psImage *image = readout->image;
+    psImage *mask  = readout->mask;
+
+    // replace the background model
+    for (int j = 0; j < image->numRows; j++) {
+	for (int i = 0; i < image->numCols; i++) {
+	    if (!mask->data.U8[j][i]) {
+		image->data.F32[j][i] += background->data.F32[j][i];
+	    }
+	}
+    }
+    
+    return true;
+}
+
Index: /trunk/psphot/src/psphotSourceFits.c
===================================================================
--- /trunk/psphot/src/psphotSourceFits.c	(revision 6480)
+++ /trunk/psphot/src/psphotSourceFits.c	(revision 6481)
@@ -126,8 +126,24 @@
     psTrace ("psphot.blend", 5, "blob as DBL: %f %f\n", ONE->params->data.F32[2], ONE->params->data.F32[3]);
 
-    source->modelPSF = (pmModel *) DBL->data[0];
-    source->mode |=  PM_SOURCE_SUBTRACTED;
-    source->mode &= ~PM_SOURCE_TEMPSUB;
-
+    // drop old model, save new second model...
+    psFree (source->modelPSF);
+    source->modelPSF = psMemCopy (DBL->data[0]);
+    source->mode    |= PM_SOURCE_SUBTRACTED;
+    source->mode    |= PM_SOURCE_PAIR;
+    source->mode    &= ~PM_SOURCE_TEMPSUB;
+
+    // copy most data from the primary source (modelEXT, blends stay NULL)
+    pmSource *newSrc = pmSourceAlloc ();
+    newSrc->peak     = psMemCopy (source->peak);
+    newSrc->pixels   = psMemCopy (source->pixels);
+    newSrc->weight   = psMemCopy (source->weight);
+    newSrc->mask     = psMemCopy (source->mask);
+    newSrc->moments  = psMemCopy (source->moments);
+    newSrc->modelPSF = psMemCopy (DBL->data[1]);
+    newSrc->type     = source->type;
+    newSrc->mode     = source->mode;
+    psArrayAdd (sources, 100, newSrc);
+    psFree (newSrc);
+    psFree (DBL);
     return true;
 }
@@ -154,5 +170,5 @@
     // save the PSF model from the Ensemble fit
     PSF = source->modelPSF;
-    psphotCheckRadiusPSF (readout, source, PSF);
+    psphotCheckRadiusPSFBlend (readout, source, PSF, 8.0);
 
     modelSet = psArrayAlloc (2);
@@ -160,16 +176,16 @@
     DBL = pmModelCopy (PSF);
     DBL->params->data.F32[1] *= 0.5;
-    DBL->params->data.F32[2] += dx;
-    DBL->params->data.F32[3] += dy;
+    DBL->params->data.F32[2] = source->moments->x + dx;
+    DBL->params->data.F32[3] = source->moments->y + dy;
     modelSet->data[0] = DBL;
 
     DBL = pmModelCopy (PSF);
     DBL->params->data.F32[1] *= 0.5;
-    DBL->params->data.F32[2] -= dx;
-    DBL->params->data.F32[3] -= dy;
+    DBL->params->data.F32[2] = source->moments->x - dx;
+    DBL->params->data.F32[3] = source->moments->y - dy;
     modelSet->data[1] = DBL;
 
-    x = PSF->params->data.F32[2];
-    y = PSF->params->data.F32[3];
+    x = source->moments->x;
+    y = source->moments->y;
 
     // fit EXT (not PSF) model (set/unset the pixel mask)
@@ -214,7 +230,4 @@
     // save the PSF model from the Ensemble fit
     pmModel *PSF = pmModelCopy (source->modelPSF);
-
-    // extend source radius as needed
-    psphotCheckRadiusPSF (readout, source, PSF);
 
     x = PSF->params->data.F32[2];
@@ -229,15 +242,11 @@
     psArrayAdd (sourceSet, 16, source);
 
-    // counter to track the blend elements used in the fit
-    psVector *index  = psVectorAlloc (source->blends->n + 1, PS_TYPE_S32);
-    index->data.S32[0] = -1; // first element corresponds to the primary source
-
+    // we need to include all blends in the fit (unless primary is saturated?)
+    dR = 0;
     for (int i = 0; i < source->blends->n; i++) {
 	pmSource *blend = source->blends->data[i];
 
-	// is this object close enough to include in the fit
-	// XXX this test is bogus: need to think about this
-	dR = hypot (blend->peak->x - x, blend->peak->y - y);
-	if (dR > PSF->radius - 3) continue;
+	// find the blend which is furthest from source
+	dR = PS_MAX (dR, hypot (blend->peak->x - x, blend->peak->y - y));
 
 	// create the model and guess parameters for this blend
@@ -256,8 +265,10 @@
 
 	// add this blend to the list
-	index->data.S32[modelSet->n] = i;
 	psArrayAdd (modelSet, 16, model);
 	psArrayAdd (sourceSet, 16, blend);
     }
+
+    // extend source radius as needed
+    psphotCheckRadiusPSFBlend (readout, source, PSF, dR);
 
     // fit PSF model (set/unset the pixel mask)
@@ -273,6 +284,5 @@
 	// if we skip this one, free the corresponding blend entry model
 	if (!psphotEvalPSF (src, model)) {
-	    int n = index->data.S32[i];
-	    pmSource *blend = source->blends->data[n];
+	    pmSource *blend = source->blends->data[i - 1];
 	    psFree (blend->modelPSF);
 	    blend->modelPSF = NULL;
@@ -285,5 +295,4 @@
 	src->mode &= ~PM_SOURCE_TEMPSUB;
     }
-    psFree (index);
     psFree (modelSet);
     psFree (sourceSet);