Index: trunk/psphot/src/psphotSourceFits.c
===================================================================
--- trunk/psphot/src/psphotSourceFits.c	(revision 6427)
+++ 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);