Index: /trunk/ppStack/src/ppStackSources.c
===================================================================
--- /trunk/ppStack/src/ppStackSources.c	(revision 18352)
+++ /trunk/ppStack/src/ppStackSources.c	(revision 18353)
@@ -15,4 +15,5 @@
     psRegion *bounds;                   // Bounding box for sources
     psArray *sources;                   // Sources within region
+    psVector *indices;                  // Indices from tree to sources
     psTree *tree;                       // kd tree with sources
 } ppStackSourceList;
@@ -42,4 +43,5 @@
     psFree(list->bounds);
     psFree(list->sources);
+    psFree(list->indices);
     psFree(list->tree);
 }
@@ -48,5 +50,6 @@
 // Returns number of valid sources
 static long stackSourcesParse(psRegion **bounds, // Region to update with bounding box
-                              psVector **x, psVector **y, // Vectors to return
+                              psVector **x, psVector **y, // Coordinate vectors to return
+                              psVector **indices, // Source indices to return
                               const psArray *sources // Input sources
     )
@@ -54,18 +57,19 @@
     psAssert(bounds, "Must be given a region for bounding box");
     psAssert(x && y, "Must be given vectors");
+    psAssert(indices, "Must be given indices");
     psAssert(sources, "Must be given sources");
 
-    long num = sources->n;              // Number of sources
-    *x = psVectorRecycle(*x, num, PS_TYPE_F32);
-    *y = psVectorRecycle(*y, num, PS_TYPE_F32);
+    long numSources = sources->n;              // Number of sources
+    *x = psVectorRecycle(*x, numSources, PS_TYPE_F32);
+    *y = psVectorRecycle(*y, numSources, PS_TYPE_F32);
+    *indices = psVectorRecycle(*indices, numSources, PS_TYPE_U32);
     float xMin = INFINITY, xMax = -INFINITY, yMin = INFINITY, yMax = -INFINITY; // Bounds of sources
-    long numValid = 0;                  // Number of valid sources
-    for (long i = 0; i < num; i++) {
+    long num = 0;                       // Number of valid sources
+    for (long i = 0; i < numSources; i++) {
         pmSource *source = sources->data[i]; // Source of interest
-        if (!source) {
+        if (!source || source->mode & SOURCE_MASK) {
             continue;
         }
 
-        numValid++;
         float xSrc, ySrc;               // Coordinates of source
         coordsFromSource(&xSrc, &ySrc, source);
@@ -74,5 +78,13 @@
         if (ySrc < yMin) yMin = ySrc;
         if (ySrc > yMax) yMax = ySrc;
-    }
+
+        (*x)->data.F32[num] = xSrc;
+        (*y)->data.F32[num] = ySrc;
+        (*indices)->data.U32[num] = i;
+        num++;
+    }
+    (*x)->n = num;
+    (*y)->n = num;
+    (*indices)->n = num;
 
     if (*bounds) {
@@ -85,5 +97,8 @@
     }
 
-    return numValid;
+    psTrace("ppStack.sources", 8, "%ld sources: bounds are [%.2f:%.2f,%.2f:%.2f]\n",
+            num, xMin, xMax, yMin, yMax);
+
+    return num;
 }
 
@@ -100,7 +115,8 @@
     list->sources = psMemIncrRefCounter(sources);
     list->bounds = NULL;
+    list->indices = NULL;
 
     psVector *x = NULL, *y = NULL;      // Source coordinates
-    stackSourcesParse(&list->bounds, &x, &y, sources);
+    stackSourcesParse(&list->bounds, &x, &y, &list->indices, sources);
     list->tree = psTreePlant(2, SOURCES_MAX_LEAF, x, y); // kd Tree with sources
     psFree(x);
@@ -123,8 +139,11 @@
 
     list->sources = oldSources = psArrayRealloc(oldSources, numOld + numNew);
+    for (long i = 0; i < numNew; i++) {
+        psArrayAdd(oldSources, 1, newSources->data[i]);
+    }
     psFree(list->tree);
 
     psVector *x = NULL, *y = NULL;      // Source coordinates
-    stackSourcesParse(&list->bounds, &x, &y, list->sources);
+    stackSourcesParse(&list->bounds, &x, &y, &list->indices, list->sources);
     list->tree = psTreePlant(2, SOURCES_MAX_LEAF, x, y); // kd Tree with sources
     psFree(x);
@@ -158,13 +177,15 @@
     psArray *sources = *sourcesPtr;     // Sources to compare with list
     psVector *x = NULL, *y = NULL;      // Coordinates of sources
+    psVector *indices = NULL;           // Indices for sources
     psRegion *bounds = NULL;            // Bounding box for sources
-    long numSources = sources->n;       // Number of sources
-    stackSourcesParse(&bounds, &x, &y, sources);
+    long numSources = stackSourcesParse(&bounds, &x, &y, &indices, sources); // Number of (good) sources
 
     if (bounds->x0 > list->bounds->x1 || bounds->x1 < list->bounds->x0 ||
         bounds->x0 > list->bounds->x1 || bounds->x1 < list->bounds->x0) {
         // Bounds don't overlap, so the sources don't either
+        psTrace("ppStack.sources", 7, "Bounds don't overlap\n");
         psFree(x);
         psFree(y);
+        psFree(indices);
         psFree(bounds);
         return 0;
@@ -179,8 +200,4 @@
     psArray *outside = psArrayAllocEmpty(numSources); // Sources that don't correlate
     for (long i = 0; i < numSources; i++) {
-        pmSource *source = sources->data[i]; // New source
-        if (!source) {
-            continue;
-        }
         coords->data.F32[0] = x->data.F32[i];
         coords->data.F32[1] = y->data.F32[i];
@@ -189,10 +206,11 @@
             numOut++;
             if (outside) {
-                psArrayAdd(outside, 1, source);
+                psArrayAdd(outside, 1, sources->data[indices->data.U32[i]]);
             }
             continue;
         }
 
-        pmSource *listSource = list->sources->data[match]; // Source in list
+        pmSource *listSource = list->sources->data[list->indices->data.U32[match]]; // Source in list
+        pmSource *source = sources->data[indices->data.U32[i]]; // Source of interest
         float listMag = listSource->psfMag; // Magnitude of source in list
         float sourceMag = source->psfMag; // Magnitude of source
@@ -207,4 +225,6 @@
     psFree(x);
     psFree(y);
+
+    psTrace("ppStack.sources", 7, "%ld sources (vs %d) overlap list %d\n", numIn, minOverlap, listIndex);
 
     if (numIn > minOverlap) {
@@ -220,4 +240,5 @@
         }
         float meanDiff = stats->clippedMean; // Mean magnitude difference
+        psTrace("ppStack.sources", 7, "Magnitude difference is %f\n", meanDiff);
         psFree(stats);
 
@@ -235,4 +256,6 @@
 
             // Suck sources from this list into the first one we found.
+            psTrace("ppStack.sources", 4, "Bridge found: Merging lists.");
+
             ppStackSourceListExtend(*merge, list->sources);
             psFree(list);
@@ -245,4 +268,5 @@
             }
             // Put the sources that didn't correlate into the list
+            psTrace("ppStack.sources", 4, "Extending list");
             ppStackSourceListExtend(list, outside);
             *merge = list;
@@ -306,4 +330,5 @@
     if (!merge) {
         // The source list is disjoint, so throw them into a new list
+        psTrace("ppStack.sources", 4, "Disjoint source: Creating new list");
         ppStackSourceList *list = ppStackSourceListAlloc(sources);
         psArrayAdd(lists, 1, list);
@@ -371,4 +396,5 @@
             if (numIn == 0) {
                 // It doesn't match anything at all.  Merge it in to the first list.
+                psTrace("ppStack.sources", 4, "Disjoint list: Merging into first");
                 ppStackSourceListExtend(firstList, list->sources);
             }
@@ -380,6 +406,10 @@
         FILE *srcFile = fopen("sources.txt", "w");
         for (long i = 0; i < firstList->sources->n; i++) {
+            pmSource *source = firstList->sources->data[i];
+            if (source->mode & SOURCE_MASK) {
+                continue;
+            }
             float x, y;
-            coordsFromSource(&x, &y, firstList->sources->data[i]);
+            coordsFromSource(&x, &y, source);
             fprintf(srcFile, "%f\t%f\n", x, y);
         }