Index: trunk/ppTranslate/src/ppMops.c
===================================================================
--- trunk/ppTranslate/src/ppMops.c	(revision 28211)
+++ trunk/ppTranslate/src/ppMops.c	(revision 28212)
@@ -20,17 +20,15 @@
     }
 
-    ppMopsDetections *merged = ppMopsMerge(detections); // Merged detections
-    psFree(detections);
-    if (!merged) {
+    if (!ppMopsPurgeDuplicates(detections)) {
         psErrorStackPrint(stderr, "Unable to merge detections");
         exit(PS_EXIT_SYS_ERROR);
     }
 
-    if (!ppMopsWrite(merged, args)) {
+    if (!ppMopsWrite(detections, args)) {
         psErrorStackPrint(stderr, "Unable to write detections");
         exit(PS_EXIT_SYS_ERROR);
     }
 
-    psFree(merged);
+    psFree(detections);
     psFree(args);
 
Index: trunk/ppTranslate/src/ppMops.h
===================================================================
--- trunk/ppTranslate/src/ppMops.h	(revision 28211)
+++ trunk/ppTranslate/src/ppMops.h	(revision 28212)
@@ -51,5 +51,5 @@
 
 /// Purge duplicate detections
-ppMopsDetections *ppMopsPurgeDuplicates(const psArray *detections);
+bool ppMopsPurgeDuplicates(const psArray *detections);
 
 /// Write detections
Index: trunk/ppTranslate/src/ppMopsDetections.c
===================================================================
--- trunk/ppTranslate/src/ppMopsDetections.c	(revision 28211)
+++ trunk/ppTranslate/src/ppMopsDetections.c	(revision 28212)
@@ -38,4 +38,5 @@
     det->mjd = NAN;
     det->seeing = NAN;
+    det->num = 0;
     det->header = NULL;
     det->table = NULL;
Index: trunk/ppTranslate/src/ppMopsMerge.c
===================================================================
--- trunk/ppTranslate/src/ppMopsMerge.c	(revision 28211)
+++ trunk/ppTranslate/src/ppMopsMerge.c	(revision 28212)
@@ -22,29 +22,44 @@
     )
 {
-    float dx = detections->x->data.F32[index] - detections->naxis1->data.S32[index] / 2.0;
-    float dy = detections->y->data.F32[index] - detections->naxis2->data.S32[index] / 2.0;
+    float dx = detections->x->data.F32[index] - detections->naxis1 / 2.0;
+    float dy = detections->y->data.F32[index] - detections->naxis2 / 2.0;
     return PS_SQR(dx) + PS_SQR(dy);
 }
 
 
-ppMopsDetections *ppMopsPurgeDuplicates(const psArray *detections)
+bool ppMopsPurgeDuplicates(const psArray *detections)
 {
     PS_ASSERT_ARRAY_NON_NULL(detections, NULL);
 
     int numInputs = detections->n;                // Number of inputs
-    psTrace("ppMops.merge", 1, "Checking detections from %ld inputs\n", numInputs);
-
-    psArray *dupes = psArrayAlloc(numInputs); // Vector of duplicate bits for each input
+    psTrace("ppMops.merge", 1, "Checking detections from %d inputs\n", numInputs);
+
+    long total = 0;                                // Total number of sources
+    psArray *duplicates = psArrayAlloc(numInputs); // Vector of duplicate bits for each input
+    psVector *dupNum = psVectorAlloc(numInputs, PS_TYPE_U32); // Number of duplicates for each input
+    psVectorInit(dupNum, 0);
     for (int i = 0; i < numInputs; i++) {
         ppMopsDetections *det = detections->data[i]; // Detections from
-
-
-
-
-
-        dupes->data[i] = psVector
-
-    ppMopsDetections *merged = NULL;    // Merged list
-    int num = 1;                                                         // Number of merged files
+        psVector *dupes = duplicates->data[i] = psVectorAlloc(det->num, PS_TYPE_U8);
+        psVectorInit(dupes, 0);
+        total += det->num;
+    }
+
+    psVector *raMerged = psVectorAllocEmpty(total, PS_TYPE_F64); // Merged RAs
+    psVector *decMerged = psVectorAllocEmpty(total, PS_TYPE_F64); // Merged Decs
+#if 0
+    psVector *xMerged = psVectorAllocEmpty(total, PS_TYPE_F32);   // Merged x coords
+    psVector *yMerged = psVectorAllocEmpty(total, PS_TYPE_F32);   // Merged y coords
+#endif
+    psVector *sourceMerged = psVectorAllocEmpty(total, PS_TYPE_U16); // Source image of merged sources
+    psVector *indexMerged = psVectorAllocEmpty(total, PS_TYPE_U32); // Index of merged sources
+    long num = 0;                                                   // Number of merged sources
+
+    const char *raBoresight = NULL, *decBoresight = NULL; // Boresight coordinates
+    const char *filter = NULL;                    // Filter name
+    float airmass = NAN, exptime = NAN;           // Exposure details
+    double posangle = NAN, alt = NAN, az = NAN; // Telescope pointing
+    double mjd = NAN;                           // Time of exposure
+
     for (int i = 0; i < detections->n; i++) {
         ppMopsDetections *det = detections->data[i]; // Detections of interest
@@ -56,69 +71,91 @@
             continue;
         }
-        num++;
-        if (!merged) {
-            psTrace("ppMops.merge", 3, "Accepting %ld detections from input %d\n", det->num, i);
-            merged = psMemIncrRefCounter(det);
-            continue;
-        }
-        psTrace("ppMops.merge", 3, "Merging %ld detections from input %d\n", det->num, i);
-
-        // XXX compare exposure properties
-        if (strcmp(merged->raBoresight, det->raBoresight) != 0) {
-            psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Exposure RA values differ: %s vs %s",
-                    merged->raBoresight, det->raBoresight);
-            return NULL;
-        }
-        if (strcmp(merged->decBoresight, det->decBoresight) != 0) {
-            psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Exposure Dec values differ: %s vs %s",
-                    merged->decBoresight, det->decBoresight);
-            return NULL;
-        }
-        if (strcmp(merged->filter, det->filter) != 0) {
-            psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Exposure filter values differ: %s vs %s",
-                    merged->filter, det->filter);
-            return NULL;
-        }
-
-        if (fabsf(merged->airmass - det->airmass) > AIRMASS_TOL) {
-            psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Exposure airmass values differ: %f vs %f",
-                    merged->airmass, det->airmass);
-            return NULL;
-        }
-        if (fabsf(merged->exptime - det->exptime) > EXPTIME_TOL) {
-            psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Exposure exposure time values differ: %f vs %f",
-                    merged->exptime, det->exptime);
-            return NULL;
-        }
-        if (fabs(merged->posangle - det->posangle) > POSANGLE_TOL) {
-            psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Exposure position angle values differ: %f vs %f",
-                    merged->posangle, det->posangle);
-            return NULL;
-        }
-        if (fabs(merged->alt - det->alt) > BORESIGHT_TOL) {
-            psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Exposure altitude values differ: %lf vs %lf",
-                    merged->alt, det->alt);
-            return NULL;
-        }
-        if (fabs(merged->az - det->az) > BORESIGHT_TOL) {
-            psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Exposure azimuth values differ: %lf vs %lf",
-                    merged->az, det->az);
-            return NULL;
-        }
-        if (fabs(merged->mjd - det->mjd) > MJD_TOL) {
-            psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Exposure MJD values differ: %lf vs %lf",
-                    merged->mjd, det->mjd);
-            return NULL;
-        }
-
-        merged->seeing += det->seeing;  // Taking average
-
-        psTree *tree = psTreePlant(2, LEAF_SIZE, PS_TREE_SPHERICAL, merged->ra, merged->dec); // kd tree
-        if (!tree) {
-            psError(PS_ERR_UNEXPECTED_NULL, false, "Unable to generate kd tree");
-            psFree(merged);
-            return NULL;
-        }
-
+
+        // Check exposure characteristics
+        if (num == 0) {
+            raBoresight = det->raBoresight;
+            decBoresight = det->decBoresight;
+            filter = det->filter;
+            airmass = det->airmass;
+            exptime = det->exptime;
+            posangle = det->posangle;
+            alt = det->alt;
+            az = det->az;
+            continue;
+        } else {
+            if (strcmp(raBoresight, det->raBoresight) != 0) {
+                psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Exposure RA values differ: %s vs %s",
+                        raBoresight, det->raBoresight);
+                return false;
+            }
+            if (strcmp(decBoresight, det->decBoresight) != 0) {
+                psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Exposure Dec values differ: %s vs %s",
+                        decBoresight, det->decBoresight);
+                return false;
+            }
+            if (strcmp(filter, det->filter) != 0) {
+                psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Exposure filter values differ: %s vs %s",
+                        filter, det->filter);
+                return false;
+            }
+            if (fabsf(airmass - det->airmass) > AIRMASS_TOL) {
+                psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Exposure airmass values differ: %f vs %f",
+                        airmass, det->airmass);
+                return false;
+            }
+            if (fabsf(exptime - det->exptime) > EXPTIME_TOL) {
+                psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Exposure exposure time values differ: %f vs %f",
+                        exptime, det->exptime);
+                return false;
+            }
+            if (fabs(posangle - det->posangle) > POSANGLE_TOL) {
+                psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Exposure position angle values differ: %f vs %f",
+                        posangle, det->posangle);
+                return false;
+            }
+            if (fabs(alt - det->alt) > BORESIGHT_TOL) {
+                psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Exposure altitude values differ: %lf vs %lf",
+                        alt, det->alt);
+                return false;
+            }
+            if (fabs(az - det->az) > BORESIGHT_TOL) {
+                psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Exposure azimuth values differ: %lf vs %lf",
+                        az, det->az);
+                return false;
+            }
+            if (fabs(mjd - det->mjd) > MJD_TOL) {
+                psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Exposure MJD values differ: %lf vs %lf",
+                        mjd, det->mjd);
+                return false;
+            }
+        }
+
+        psTrace("ppMops.merge", 3, "Accepting %ld detections from input %d\n", det->num, i);
+        memcpy(&raMerged->data.F64[num], det->ra, det->num * PSELEMTYPE_SIZEOF(PS_TYPE_F64));
+        memcpy(&decMerged->data.F64[num], det->dec, det->num * PSELEMTYPE_SIZEOF(PS_TYPE_F64));
+        for (long j = 0, k = num; j < det->num; j++, k++) {
+            sourceMerged->data.U16[k] = i;
+            indexMerged->data.U32[k] = j;
+        }
+        num += det->num;
+    }
+
+    psTrace("ppMops.merge", 3, "Generating kd-tree from %ld sources\n", num);
+    psTree *tree = psTreePlant(2, LEAF_SIZE, PS_TREE_SPHERICAL, raMerged, decMerged); // kd tree
+    if (!tree) {
+        psError(PS_ERR_UNEXPECTED_NULL, false, "Unable to generate kd tree");
+        return false;
+    }
+
+    for (int i = 0; i < detections->n; i++) {
+        ppMopsDetections *det = detections->data[i]; // Detections of interest
+        if (!det) {
+            continue;
+        } else if (det->num == 0) {
+            continue;
+        }
+        psTrace("ppMops.merge", 3, "Checking %ld detections from input %d\n", det->num, i);
+
+        psVector *dupes = duplicates->data[i];            // Duplicates list
         psVector *coords = psVectorAlloc(2, PS_TYPE_F64); // Coordinates of interest
         for (int j = 0; j < det->num; j++) {
@@ -130,14 +167,8 @@
                 psFree(coords);
                 psFree(tree);
-                psFree(merged);
-                return NULL;
-            }
-            if (indices->n == 0) {
-                psTrace("ppMops.merge", 9, "No matches for source %d in input %d\n", j, i);
-                psFree(indices);
-                ppMopsDetectionsCopySingle(merged, det, j);
-                continue;
+                return false;
             }
             psTrace("ppMops.merge", 5, "%ld matches for source %d from input %d\n", indices->n, j, i);
+            psAssert(indices->n > 0, "Expect at least one match for source %d in input %d", j, i);
 
             // Which one do we keep?
@@ -145,6 +176,16 @@
             long bestIndex = -1;           // Index with best distance
             for (int k = 0; k < indices->n; k++) {
-                long index = indices->data.S64[k]; // Index of point
-                float distance = mergeDistance(merged, index); // Distance to centre of image
+                long mergeIndex = indices->data.S64[k]; // Index of point in merged list
+                int source = sourceMerged->data.U16[mergeIndex]; // Source image
+                if (source == i) {
+                    continue;
+                }
+                long index = indexMerged->data.U32[mergeIndex];  // Index in source
+                psVector *dupes = duplicates->data[source];     // Duplicates list
+                if (dupes->data.U8[index]) {
+                    continue;
+                }
+
+                float distance = mergeDistance(detections->data[source], index); // Distance to centre of image
                 if (distance < bestDistance) {
                     bestDistance = distance;
@@ -154,29 +195,78 @@
 
             float distance = mergeDistance(det, j); // Distance to centre of image
-            if (distance < bestDistance) {
+            if (!isfinite(bestDistance) || distance < bestDistance) {
                 psTrace("ppMops.merge", 6, "New source clobbers old sources\n");
-                // Blow away existing sources
                 for (int k = 0; k < indices->n; k++) {
-                    long index = indices->data.S64[k]; // Index of point
-                    merged->mask->data.U8[index] = 0xFF;
-                }
-                ppMopsDetectionsCopySingle(merged, det, j);
+                    long mergeIndex = indices->data.S64[k]; // Index of point
+                    int source = sourceMerged->data.U16[mergeIndex]; // Source image
+                    if (source == i) {
+                        continue;
+                    }
+                    long index = indexMerged->data.U32[mergeIndex];  // Index in source
+                    psVector *dupes = duplicates->data[source];      // Duplicates list
+                    if (!dupes->data.U8[index]) {
+                        dupes->data.U8[index] = 0xFF;
+                        dupNum->data.U32[source]++;
+                    }
+                }
             } else {
                 psTrace("ppMops.merge", 6, "Old sources clobber new source\n");
-            }
-            psFree(indices);
-        }
-
-        psTrace("ppMops.merge", 3, "Done merging input %d, %ld merged sources\n", i, merged->num);
-
-        psFree(tree);
-        ppMopsDetectionsPurge(merged);
-    }
-
-    psTrace("ppMops.merge", 2, "%ld sources in merged detections list\n", merged->num);
-
-    merged->seeing /= num;
-
-    return merged;
+                dupes->data.U8[j] = 0xFF;
+                dupNum->data.U32[i]++;
+            }
+        }
+        psFree(coords);
+    }
+    psFree(tree);
+    psFree(raMerged);
+    psFree(decMerged);
+    psFree(sourceMerged);
+    psFree(indexMerged);
+
+    // Remove duplicates
+    for (int i = 0; i < detections->n; i++) {
+        ppMopsDetections *det = detections->data[i]; // Detections of interest
+        if (!det) {
+            continue;
+        } else if (det->num == 0) {
+            continue;
+        }
+        psTrace("ppMops.merge", 3, "Purging %d duplicates from input %d\n", dupNum->data.U32[i], i);
+
+#define VECTOR_PURGE_CASE(TYPE)                                      \
+        case PS_TYPE_##TYPE:                                         \
+          for (int i = 0, j = 0; i < vector->n; i++) {               \
+              if (!dupes->data.U8[i]) {                              \
+                  if (i == j) {                                      \
+                      j++;                                           \
+                      continue;                                      \
+                  }                                                  \
+                  vector->data.TYPE[j] = vector->data.TYPE[i];       \
+              }                                                      \
+          }                                                          \
+          break;
+
+        psVector *dupes = duplicates->data[i]; // Duplicates
+        psArray *table = psListToArray(det->table->list); // Table of data
+        for (int j = 0; j < table->n; j++) {
+            psVector *vector = table->data[j]; // Vector to purge
+            switch (vector->type.type) {
+                VECTOR_PURGE_CASE(U8);
+                VECTOR_PURGE_CASE(U16);
+                VECTOR_PURGE_CASE(U32);
+                VECTOR_PURGE_CASE(U64);
+                VECTOR_PURGE_CASE(S8);
+                VECTOR_PURGE_CASE(S16);
+                VECTOR_PURGE_CASE(S32);
+                VECTOR_PURGE_CASE(S64);
+                VECTOR_PURGE_CASE(F32);
+                VECTOR_PURGE_CASE(F64);
+              default:
+                psAbort("Unrecognised vector type: %x", vector->type.type);
+            }
+        }
+    }
+
+    return true;
 }
 
Index: trunk/ppTranslate/src/ppMopsRead.c
===================================================================
--- trunk/ppTranslate/src/ppMopsRead.c	(revision 28211)
+++ trunk/ppTranslate/src/ppMopsRead.c	(revision 28212)
@@ -41,6 +41,7 @@
             continue;
         }
-        ppMopsDetections *det = ppMopsDetectionsAlloc(size);
+        ppMopsDetections *det = ppMopsDetectionsAlloc();
         det->header = header;
+        det->num = size;
 
         det->raBoresight = psMemIncrRefCounter(psMetadataLookupStr(NULL, header, "FPA.RA"));
@@ -60,5 +61,5 @@
         det->naxis2 = psMetadataLookupS32(NULL, header, "IMNAXIS2");
 
-        det->table = psFitsReadTableAllColumns(fits); // Table of interest
+        psMetadata *table = det->table = psFitsReadTableAllColumns(fits); // Table of interest
         if (!table) {
             psError(psErrorCodeLast(), false, "Unable to read table %d", i);
@@ -66,6 +67,4 @@
         }
         psFitsClose(fits);
-
-        psTrace("ppMops.read", 2, "Read %ld rows from %s\n", numGood, (const char*)inNames->data[i]);
 
         det->ra = psMemIncrRefCounter(psMetadataLookupVector(NULL, table, "RA_PSF"));
@@ -78,4 +77,6 @@
         }
 
+        psTrace("ppMops.read", 2, "Read %ld rows from %s\n", det->ra->n, (const char*)inNames->data[i]);
+
         detections->data[i] = det;
     }
Index: trunk/ppTranslate/src/ppMopsWrite.c
===================================================================
--- trunk/ppTranslate/src/ppMopsWrite.c	(revision 28211)
+++ trunk/ppTranslate/src/ppMopsWrite.c	(revision 28212)
@@ -36,5 +36,5 @@
         psMetadataAddStr(header, PS_LIST_TAIL, "OBSCODE", 0, "IAU Observatory code", OBSERVATORY_CODE);
 
-        if (!psFitsWriteBlank(fits, header)) {
+        if (!psFitsWriteBlank(fits, header, NULL)) {
             psError(psErrorCodeLast(), false, "Unable to write primary header");
             return false;