Index: /trunk/psphot/src/psphotBasicDeblend.c
===================================================================
--- /trunk/psphot/src/psphotBasicDeblend.c	(revision 5703)
+++ /trunk/psphot/src/psphotBasicDeblend.c	(revision 5704)
@@ -6,70 +6,67 @@
 { 
 
-    // sources must be sorted by Y
-    // index must be sorted by SN
+    // source analysis is done in S/N order (brightest first)
+    sources = psArraySort (sources, psphotSortByY);
 
-    // search for sources which fall within current box
+    psVector *SN = psVectorAlloc (sources->n, PS_DATA_F32);
+    for (int i = 0; i < SN->n; i++) {
+      source = sources->data[i];
+      SN->data.F32[i] = source->moments->SN;
+    }
+    psVector *index = psVectorSortIndex (NULL, SN);
+    
+    psArray *overlap = psArrayAlloc (100);
 
-    // generate source contour (1/4 peak counts)
-    
-    // find sources from overlap which fall in contour
-    // mark source
+    for (int i = 0; i < sources->n; i++) {
+      N = index->data.U32[i];
+      source = sources->data[N];
+      overlap->n = 0;
 
-    // select the candidate PSF stars (pointers to original sources)
-    for (int i = 0; (i < sources->n) && (stars->n < NSTARS); i++) {
-	pmSource *source = sources->data[i];
-	if (source->type != PM_SOURCE_PSFSTAR) continue;
-	psArrayAdd (stars, 200, source);
-    }
-    psTrace (".psphot.pspsf", 3, "selected candidate %d PSF objects\n", stars->n);
+      // search backwards for overlapping sources
+      for (int j = N - 1; j >= 0; j--) {
+	testSource = sources->data[j];
+	if (testSource->moments->x < source->pixels->col0) continue;
+	if (testSource->moments->x >= source->pixels->col0 + source->pixels->numCols) continue;
+	if (testSource->moments->y < source->pixels->row0) break;
+	if (testSource->moments->y >= source->pixels->row0 + source->pixels->numRows) {
+	  fprintf (stderr, "warning: invalid condition\n");
+	  continue;
+	}
+	psArrayAdd (overlap, 100, testSource);
+      }
 
-    // get the fixed PSF fit radius
-    // XXX EAM : check that PSF_FIT_RADIUS < SKY_OUTER_RADIUS
-    float RADIUS = psMetadataLookupF32 (&status, config, "PSF_FIT_RADIUS");
+      // search forwards for overlapping sources
+      for (int j = N + 1; j < sources->n; j++) {
+	testSource = sources->data[j];
+	if (testSource->moments->x < source->pixels->col0) continue;
+	if (testSource->moments->x >= source->pixels->col0 + source->pixels->numCols) continue;
+	if (testSource->moments->y < source->pixels->row0) {
+	  fprintf (stderr, "warning: invalid condition\n");
+	  continue;
+	}
+	if (testSource->moments->y >= source->pixels->row0 + source->pixels->numRows) break;
+	psArrayAdd (overlap, 100, testSource);
+      }
 
-    // get the list pointers for the PSF_MODEL entries
-    psMetadataItem *mdi = psMetadataLookup (config, "PSF_MODEL");
-    if (mdi == NULL) psAbort ("psphotChoosePSF", "missing PSF_MODEL selection");
+      // generate source contour (1/4 peak counts)
+      if (overlap->n > 0) {
+	// XXX EAM : make the 1/4 user-defined.
+	threshold = 0.25 * (source->moments->Peak - source->moments->Sky) + source->moments->Sky); 
+	psArray *contour = pmSourceContour_EAM (source->pixels, source->moments->x, source->moments->x, threshold);
 
-    psList *list = (psList *) mdi->data.list;
-    psListIterator *iter = psListIteratorAlloc (list, PS_LIST_HEAD, FALSE); 
-
-    // set up an array to store the results
-    psArray *models = psArrayAlloc (list->n);
-
-    // try each model option listed in config
-    for (int i = 0; i < models->n; i++) { 
-
-	item = psListGetAndIncrement (iter);
-	modelName = item->data.V;
-
-	models->data[i] = pmPSFtryModel (stars, modelName, RADIUS);
-	psFree (modelName);
-	psFree (item);
-    }
-    psFree (iter);
-    // psFree (list); XXX EAM - is list freed with iter?
-    psFree (stars);
-
-    // select the best of the models
-    // here we are using the clippedStdev on the metric as the indicator
-    try = models->data[0];
-    int bestN = 0;
-    float bestM = try->psf->dApResid;
-    for (int i = 1; i < models->n; i++) {
-	try = models->data[i];
-	float M = try->psf->dApResid;
-	if (M < bestM) {
-	    bestM = M;
-	    bestN = i;
-	}
+	// XXX EAM : the order is wrong: need to check j, j+1 contour points
+	psVector *xv = contour->data[0];
+	psVector *yv = contour->data[1];
+	for (int j = 0; j < xv->n; j++) {
+	  for (int k = 0; k < overlap->n; k++) {
+	    testSource = overlap->data[k];
+	    if (fabs(yv->data.F32[j] - testSource->moments->y) > 0.5) continue;
+	    if (xv->data.F32[j] > testSource->moments->x) continue;
+	    if (xv->data.F32[j] < testSource->moments->x) continue;
+	    // XXX EAM : mark source with flag (blended)
+	  }
+	}  
     }
 
-    // keep only the selected model:
-    try = models->data[bestN];
-    pmPSF *psf = psMemCopy(try->psf);
-    psFree (models);				 // keep only the pmPSF resulting from this analysis
-
-    modelName = pmModelGetType (psf->type);
     psLogMsg ("psphot.pspsf", 3, "selected psf model %s, ApResid: %f +/- %f\n", modelName, psf->ApResid, psf->dApResid);
 
Index: /trunk/psphot/src/psphotEnsemblePSF.c
===================================================================
--- /trunk/psphot/src/psphotEnsemblePSF.c	(revision 5703)
+++ /trunk/psphot/src/psphotEnsemblePSF.c	(revision 5704)
@@ -96,4 +96,6 @@
 	    pmSource *Mj = models->data[j];
 
+	    // XXX double check that this is working!  should not break here
+	    check me here;
 	    if (Mi->pixels->col0 + Mi->pixels->numCols < Mj->pixels->col0) break;
 	    if (Mj->pixels->col0 + Mj->pixels->numCols < Mi->pixels->col0) continue;