Index: trunk/ppTranslate/src/ppMopsRead.c
===================================================================
--- trunk/ppTranslate/src/ppMopsRead.c	(revision 32406)
+++ trunk/ppTranslate/src/ppMopsRead.c	(revision 34279)
@@ -8,201 +8,346 @@
 #include "ppMops.h"
 
+static psMetadata* sort_by_increasing_idet(psMetadata *fittedTrail, long length, long fillLength);
+
 /*
   ppMopsRead possibly modifies the args->version if the user did not
   set it explicitely.
- */
+*/
 psArray *ppMopsRead(ppMopsArguments *args)
 {
-    psTrace("ppMops.read", 1, "Reading input detections\n");
-
-    psArray *inNames = args->input;          // Input names
-    long num = inNames->n;                   // Number of inputs
-    psArray *detections = psArrayAlloc(num); // Array of detections, to return
-    for (int i = 0; i < num; i++) {
-        const char *name = inNames->data[i];
-
-        psFits *fits = psFitsOpen(name,  "r"); // FITS file
-        if (!fits) {
-            psError(PS_ERR_IO, false, "Unable to open input %d", i);
-            return false;
-        }
-
-        psMetadata *header = psFitsReadHeader(NULL, fits); // Primary header
-        if (!header) {
-            psError(PS_ERR_IO, false, "Unable to read header %d", i);
-            return false;
-        }
-
-        psS64 diffSkyfileId = psMetadataLookupS64(NULL, header, "IMAGEID"); // Identifier for image
-        if (diffSkyfileId == 0) {
-            psError(PS_ERR_UNEXPECTED_NULL, false, "Unable to find identifier for image %d", i);
-            return false;
-        }
-
-        if (!psFitsMoveExtName(fits, "SkyChip.psf")) {
-            psError(PS_ERR_IO, false, "Unable to move to HDU with detections");
-            return false;
-        }
-	int skyChipPsfVersion = ppMopsGetSkyChipPsfVersion(fits);
-	if (args->version == 0) {
-	  psTrace("ppMops.read", 1, "Changing args->version to %d\n", skyChipPsfVersion);
-	  args->version = (unsigned short) skyChipPsfVersion;
+  psTrace("ppMops.read", 1, "Reading input detections\n");
+
+  psArray *inNames = args->input;          // Input names
+  long num = inNames->n;                   // Number of inputs
+  psArray *detections = psArrayAlloc(num); // Array of detections, to return
+  for (int i = 0; i < num; i++) {
+    const char *name = inNames->data[i];
+
+    psFits *fits = psFitsOpen(name,  "r"); // FITS file
+    if (!fits) {
+      psError(PS_ERR_IO, false, "Unable to open input %d", i);
+      return false;
+    }
+
+    psMetadata *header = psFitsReadHeader(NULL, fits); // Primary header
+    if (!header) {
+      psError(PS_ERR_IO, false, "Unable to read header %d", i);
+      return false;
+    }
+
+    psS64 diffSkyfileId = psMetadataLookupS64(NULL, header, "IMAGEID"); // Identifier for image
+    if (diffSkyfileId == 0) {
+      psError(PS_ERR_UNEXPECTED_NULL, false, "Unable to find identifier for image %d", i);
+      return false;
+    }
+
+    if (!psFitsMoveExtName(fits, "SkyChip.psf")) {
+      psError(PS_ERR_IO, false, "Unable to move to HDU with detections");
+      return false;
+    }
+    int skyChipPsfVersion = ppMopsGetSkyChipPsfVersion(fits);
+    if (args->version == 0) {
+      psTrace("ppMops.read", 1, "Changing args->version to %d\n", skyChipPsfVersion);
+      args->version = (unsigned short) skyChipPsfVersion;
+    }
+    if (skyChipPsfVersion == 0) {
+      // Try to read with the user specified version?
+      skyChipPsfVersion = args->version;
+    }
+    /* Display a warning message if there are version
+       inconsistencies between the file and the flag (note that
+       those inconsistencies might be wanted) */
+    if (skyChipPsfVersion != args->version) {
+      if (skyChipPsfVersion > args->version) {
+	psWarning("The FITS data will be downgraded from PS1_DV%d to PS1_DV%d\n",
+		  skyChipPsfVersion, args->version);
+      } else { // Necessarily: skyChipPsfVersion > args->version
+	psWarning("The FITS data will be upgraded from PS1_DV%d to PS1_DV%d (new values set to default 0, NaN...)\n",
+		  skyChipPsfVersion, args->version);
+      }
+    }
+
+    long size = psFitsTableSize(fits); // Size of table
+    if (size <= 0) {
+      psErrorStackPrint(stderr, "Unable to determine size of table %d", i);
+      psErrorClear();
+      psWarning("Ignoring input %d", i);
+      psFree(header);
+      psFitsClose(fits);
+      continue;
+    }
+    ppMopsDetections *det = ppMopsDetectionsAlloc(size);
+    detections->data[i] = det;
+    det->component = psStringNCopy(name, strrchr(name, '.') - name); // Strip off extension
+    det->num = size;
+    det->diffSkyfileId = diffSkyfileId;
+
+    psTrace("ppMops.read", 3, "Reading %ld rows from %s\n", size, (const char*)inNames->data[i]);
+
+    det->raBoresight = psMemIncrRefCounter(psMetadataLookupStr(NULL, header, "FPA.RA"));
+    det->decBoresight = psMemIncrRefCounter(psMetadataLookupStr(NULL, header, "FPA.DEC"));
+    det->filter = psMemIncrRefCounter(psMetadataLookupStr(NULL, header, "FPA.FILTER"));
+    det->airmass = psMetadataLookupF32(NULL, header, "AIRMASS");
+    det->exptime = psMetadataLookupF32(NULL, header, "EXPTIME");
+    det->posangle = psMetadataLookupF64(NULL, header, "FPA.POSANGLE");
+    det->alt = psMetadataLookupF64(NULL, header, "FPA.ALT");
+    det->az = psMetadataLookupF64(NULL, header, "FPA.AZ");
+    det->mjd = psMetadataLookupF64(NULL, header, "MJD-OBS") + det->exptime / 2.0 / 3600 / 24;
+
+    det->seeing = (float) 0.5 * (psMetadataLookupF32(NULL, header, "FWHM_MAJ") +
+				 psMetadataLookupF32(NULL, header, "FWHM_MIN"));
+
+    det->naxis1 = psMetadataLookupS32(NULL, header, "IMNAXIS1"); // Number of columns
+    det->naxis2 = psMetadataLookupS32(NULL, header, "IMNAXIS2"); // Number of rows
+
+    psFree(header);
+
+    psMetadata *table = psFitsReadTableAllColumns(fits); // Table of interest
+    if (!table) {
+      psError(PS_ERR_IO, false, "Unable to read table %d", i);
+      return NULL;
+    }
+    //Sort table by idet if necessary
+    psMetadata* sortedTable = sort_by_increasing_idet(table, det->num, det->num);
+    if (sortedTable != table) {
+      psFree(table);
+      table = sortedTable;
+    } // else det->table is sorted
+    det->table = table;
+    psTrace("ppMops.read", 9, "PSF table:\n%s\n", psMetadataConfigFormat(det->table));
+    psTrace("ppMops.read", 9, "End of PSF table\n");
+
+    //fittedTrail extension (Supported for releases 2 and above)
+    if (skyChipPsfVersion >= 2) {
+      if (!psFitsMoveExtName(fits, "SkyChip.xfit")) {
+	psTrace("ppMops.read", 3, "No fitted trails extension");
+      } else {
+	psTrace("ppMops.read", 3, "Fitted trails extension found\n");
+	det->fittedTrailsSize = psFitsTableSize(fits);
+	if (det->fittedTrailsSize <= 0) {
+	  psErrorStackPrint(stderr, "Unable to determine size of fitted trails extension table %d", i);
+	  psTrace("ppMops.read", 3, "No entry in fitted trails extension!!!!\n");
+	  psErrorClear();
+	  det->fittedTrails = NULL;
+	} else {
+	  det->fittedTrails = psFitsReadTableAllColumns(fits); // Table of interest
+	  if (!det->fittedTrails) {
+	    psError(PS_ERR_IO, false, "Unable to read fittedTrails table %d", i);
+	    return NULL;
+	  }
+	  psMetadata* sortedFittedTrails = sort_by_increasing_idet(det->fittedTrails, det->fittedTrailsSize, det->num);
+	  if (sortedFittedTrails != det->fittedTrails) {
+	    psFree(det->fittedTrails);
+	    det->fittedTrails = sortedFittedTrails;
+	  } // else det->fittedTrails is sorted
+	  psTrace("ppMops.read", 9, "Fitted trail:\n%s\n", psMetadataConfigFormat(det->fittedTrails));
+	  psTrace("ppMops.read", 9, "End of Fitted trail\n");
 	}
-	if (skyChipPsfVersion == 0) {
-	  // Try to read with the user specified version?
-	  skyChipPsfVersion = args->version;
-	}
-	/* Display a warning message if there are version
-	   inconsistencies between the file and the flag (note that
-	   those inconsistencies might be wanted) */
-	if (skyChipPsfVersion != args->version) {
-	  if (skyChipPsfVersion > args->version) {
-	    psWarning("The FITS data will be downgraded from PS1_DV%d to PS1_DV%d\n",
-		      skyChipPsfVersion, args->version);
-	  } else { // Necessarily: skyChipPsfVersion > args->version
-	    psWarning("The FITS data will be upgraded from PS1_DV%d to PS1_DV%d (new values set to default 0, NaN...)\n",
-		      skyChipPsfVersion, args->version);	    
-	  }
-	}
-
-        long size = psFitsTableSize(fits); // Size of table
-        if (size <= 0) {
-            psErrorStackPrint(stderr, "Unable to determine size of table %d", i);
-            psErrorClear();
-            psWarning("Ignoring input %d", i);
-            psFree(header);
-            psFitsClose(fits);
-            continue;
-        }
-        ppMopsDetections *det = ppMopsDetectionsAlloc(size);
-        detections->data[i] = det;
-        det->component = psStringNCopy(name, strrchr(name, '.') - name); // Strip off extension
-        det->num = size;
-        det->diffSkyfileId = diffSkyfileId;
-
-        psTrace("ppMops.read", 3, "Reading %ld rows from %s\n", size, (const char*)inNames->data[i]);
-
-        det->raBoresight = psMemIncrRefCounter(psMetadataLookupStr(NULL, header, "FPA.RA"));
-        det->decBoresight = psMemIncrRefCounter(psMetadataLookupStr(NULL, header, "FPA.DEC"));
-        det->filter = psMemIncrRefCounter(psMetadataLookupStr(NULL, header, "FPA.FILTER"));
-        det->airmass = psMetadataLookupF32(NULL, header, "AIRMASS");
-        det->exptime = psMetadataLookupF32(NULL, header, "EXPTIME");
-        det->posangle = psMetadataLookupF64(NULL, header, "FPA.POSANGLE");
-        det->alt = psMetadataLookupF64(NULL, header, "FPA.ALT");
-        det->az = psMetadataLookupF64(NULL, header, "FPA.AZ");
-        det->mjd = psMetadataLookupF64(NULL, header, "MJD-OBS") + det->exptime / 2.0 / 3600 / 24;
-
-        det->seeing = (float) 0.5 * (psMetadataLookupF32(NULL, header, "FWHM_MAJ") +
-				     psMetadataLookupF32(NULL, header, "FWHM_MIN"));
-
-        det->naxis1 = psMetadataLookupS32(NULL, header, "IMNAXIS1"); // Number of columns
-        det->naxis2 = psMetadataLookupS32(NULL, header, "IMNAXIS2"); // Number of rows
-
-        psFree(header);
-
-        psMetadata *table = psFitsReadTableAllColumns(fits); // Table of interest
-        if (!table) {
-            psError(PS_ERR_IO, false, "Unable to read table %d", i);
-            return NULL;
-        }
-        det->table = table;
-        psFitsClose(fits);
-        if (args->version == 0) {
-          if (skyChipPsfVersion < 2) {
-            // XXX: TODO: Do we need to add dummy vectors for the missing columns?
-           }
-        }
-
-        psVector *ra = psMetadataLookupVector(NULL, table, "RA_PSF");
-        psVector *dec = psMetadataLookupVector(NULL, table, "DEC_PSF");
-
-        det->raErr = psVectorAlloc(size, PS_TYPE_F64);
-        det->decErr = psVectorAlloc(size, PS_TYPE_F64);
-        det->mask = psVectorAlloc(size, PS_TYPE_U8);
-
-        // convert ra and dec to radians for use in the purge duplicates function
-        det->ra = (psVector*)psBinaryOp(NULL, ra, "*", psScalarAlloc(DEG_TO_RAD(1.0), PS_TYPE_F64));
-        det->dec = (psVector*)psBinaryOp(NULL, dec, "*", psScalarAlloc(DEG_TO_RAD(1.0), PS_TYPE_F64));
-        det->x = psMemIncrRefCounter(psMetadataLookupVector(NULL, table, "X_PSF"));
-        det->y = psMemIncrRefCounter(psMetadataLookupVector(NULL, table, "Y_PSF"));
-        if (!det->ra || !det->dec || !det->x || !det->y) {
-            psError(PS_ERR_UNEXPECTED_NULL, true, "Unable to find all of RA, Dec, X and Y columns");
-            return NULL;
-        }
-
-        // Add our new vectors to the table so that duplicates and masked items may be purged
-        psMetadataAddVector(table, PS_LIST_HEAD, "DEC_ERR", 0, NULL, det->decErr);
-        psMetadataAddVector(table, PS_LIST_HEAD, "RA_ERR", 0, NULL, det->raErr);
-
-        psTrace("ppMops.read", 2, "Read %ld rows from %s\n", det->num, det->component);
-
-        psVector *mag    = psMetadataLookupVector(NULL, table, "PSF_INST_MAG");
-        psVector *magErr = psMetadataLookupVector(NULL, table, "PSF_INST_MAG_SIG");
-        psVector *xErrV  = psMetadataLookupVector(NULL, table, "X_PSF_SIG");
-        psVector *yErrV  = psMetadataLookupVector(NULL, table, "Y_PSF_SIG");
-        psVector *scaleV = psMetadataLookupVector(NULL, table, "PLTSCALE");
-        psVector *angleV = psMetadataLookupVector(NULL, table, "POSANGLE");
-        psVector *flagsV = psMetadataLookupVector(NULL, table, "FLAGS");
-
-        double plateScale = 0.0;        // Plate scale
-        long numGood = 0;               // Number of good rows
-        for (long row = 0; row < size; row++) {
-
-            psU32 flags = flagsV->data.U32[row]; // psFitsTableGetU32(NULL, table, row, "FLAGS");
-            if (flags & SOURCE_MASK) {
-                psTrace("ppMops.read", 10, "Discarding row %ld from input %d because of flags: %ud", row, i, flags);
-                det->mask->data.U8[row] = 0xFF;
-                continue;
-            }
-
-            // Calculate error in RA, Dec
-
-            double xErr = xErrV->data.F32[row];
-            double yErr = yErrV->data.F32[row];
-            double scale = scaleV->data.F32[row];
-            double angle = angleV->data.F32[row];
-
-            if (!isfinite(det->x->data.F32[row]) || !isfinite(det->y->data.F32[row]) ||
-                !isfinite(det->ra->data.F64[row]) || !isfinite(det->dec->data.F64[row]) ||
-                !isfinite(mag->data.F32[row]) || !isfinite(magErr->data.F32[row]) ||
-                !isfinite(xErr) || !isfinite(yErr) || !isfinite(scale) || !isfinite(angle)) {
-                psTrace("ppMops.read", 10,
-                        "Discarding row %ld from input %d because of non-finite values: "
-                        "%f %f %lf %lf %f %f %f %f %f %f",
-                        row, i,
-                        det->x->data.F32[row], det->y->data.F32[row],
-                        det->ra->data.F64[row], det->dec->data.F64[row],
-                        mag->data.F32[row], magErr->data.F32[row],
-                        xErr, yErr, scale, angle);
-                det->mask->data.U8[row] = 0xFF;
-                continue;
-            }
-
-            // XXX Not at all sure I've got the angles around the right way here...
-            double cosAngle = cos(angle), sinAngle = sin(angle);
-            double cosAngle2 = PS_SQR(cosAngle), sinAngle2 = PS_SQR(sinAngle);
-            double xErr2 = PS_SQR(xErr), yErr2 = PS_SQR(yErr);
-            double errScale = scale / 3600.0;
-            det->raErr->data.F64[row] = errScale * sqrt(cosAngle2 * xErr2 + sinAngle2 * yErr2);
-            det->decErr->data.F64[row] = errScale * sqrt(sinAngle2 * xErr2 + cosAngle2 * yErr2);
-
-            det->mask->data.U8[row] = 0;
-            plateScale += scale;
-            numGood++;
-        }
-        det->seeing *= ((float) plateScale) / ((float) numGood);
-
-        // Are we using numGood for anything outside of this function?
-        det->numGood = numGood;
-
-        if (isfinite(args->zp) && numGood > 0) {
-            psBinaryOp(mag, mag, "+", psScalarAlloc(args->zp, PS_TYPE_F32));
-        }
-
-        psTrace("ppMops.read", 2, "Read %ld good rows from %s\n", numGood, (const char*)name);
-    }
-
-    psTrace("ppMops.read", 1, "Done reading input detections\n");
-
-    return detections;
-}
+      }
+    }
+
+    psFitsClose(fits);
+
+    if (args->version == 0) {
+      if (skyChipPsfVersion < 2) {
+	// XXX: TODO: Do we need to add dummy vectors for the missing columns?
+      }
+    }
+
+    psVector *ra = psMetadataLookupVector(NULL, table, "RA_PSF");
+    psVector *dec = psMetadataLookupVector(NULL, table, "DEC_PSF");
+
+    det->raErr = psVectorAlloc(size, PS_TYPE_F64);
+    det->decErr = psVectorAlloc(size, PS_TYPE_F64);
+    det->mask = psVectorAlloc(size, PS_TYPE_U8);
+
+    // convert ra and dec to radians for use in the purge duplicates function
+    det->ra = (psVector*)psBinaryOp(NULL, ra, "*", psScalarAlloc(DEG_TO_RAD(1.0), PS_TYPE_F64));
+    det->dec = (psVector*)psBinaryOp(NULL, dec, "*", psScalarAlloc(DEG_TO_RAD(1.0), PS_TYPE_F64));
+    det->x = psMemIncrRefCounter(psMetadataLookupVector(NULL, table, "X_PSF"));
+    det->y = psMemIncrRefCounter(psMetadataLookupVector(NULL, table, "Y_PSF"));
+    if (!det->ra || !det->dec || !det->x || !det->y) {
+      psError(PS_ERR_UNEXPECTED_NULL, true, "Unable to find all of RA, Dec, X and Y columns");
+      return NULL;
+    }
+
+    // Add our new vectors to the table so that duplicates and masked items may be purged
+    psMetadataAddVector(table, PS_LIST_HEAD, "DEC_ERR", 0, NULL, det->decErr);
+    psMetadataAddVector(table, PS_LIST_HEAD, "RA_ERR", 0, NULL, det->raErr);
+
+    psTrace("ppMops.read", 2, "Read %ld rows from %s\n", det->num, det->component);
+
+    psVector *mag    = psMetadataLookupVector(NULL, table, "PSF_INST_MAG");
+    psVector *magErr = psMetadataLookupVector(NULL, table, "PSF_INST_MAG_SIG");
+    psVector *xErrV  = psMetadataLookupVector(NULL, table, "X_PSF_SIG");
+    psVector *yErrV  = psMetadataLookupVector(NULL, table, "Y_PSF_SIG");
+    psVector *scaleV = psMetadataLookupVector(NULL, table, "PLTSCALE");
+    psVector *angleV = psMetadataLookupVector(NULL, table, "POSANGLE");
+    psVector *flagsV = psMetadataLookupVector(NULL, table, "FLAGS");
+
+    double plateScale = 0.0;        // Plate scale
+    long numGood = 0;               // Number of good rows
+    for (long row = 0; row < size; row++) {
+
+      psU32 flags = flagsV->data.U32[row]; // psFitsTableGetU32(NULL, table, row, "FLAGS");
+      if (flags & SOURCE_MASK) {
+	psTrace("ppMops.read", 10, "Discarding row %ld from input %d because of flags: %ud", row, i, flags);
+	det->mask->data.U8[row] = 0xFF;
+	continue;
+      }
+
+      // Calculate error in RA, Dec
+
+      double xErr = xErrV->data.F32[row];
+      double yErr = yErrV->data.F32[row];
+      double scale = scaleV->data.F32[row];
+      double angle = angleV->data.F32[row];
+
+      if (!isfinite(det->x->data.F32[row]) || !isfinite(det->y->data.F32[row]) ||
+	  !isfinite(det->ra->data.F64[row]) || !isfinite(det->dec->data.F64[row]) ||
+	  !isfinite(mag->data.F32[row]) || !isfinite(magErr->data.F32[row]) ||
+	  !isfinite(xErr) || !isfinite(yErr) || !isfinite(scale) || !isfinite(angle)) {
+	psTrace("ppMops.read", 10,
+		"Discarding row %ld from input %d because of non-finite values: "
+		"%f %f %lf %lf %f %f %f %f %f %f",
+		row, i,
+		det->x->data.F32[row], det->y->data.F32[row],
+		det->ra->data.F64[row], det->dec->data.F64[row],
+		mag->data.F32[row], magErr->data.F32[row],
+		xErr, yErr, scale, angle);
+	det->mask->data.U8[row] = 0xFF;
+	continue;
+      }
+
+      // XXX Not at all sure I've got the angles around the right way here...
+      double cosAngle = cos(angle), sinAngle = sin(angle);
+      double cosAngle2 = PS_SQR(cosAngle), sinAngle2 = PS_SQR(sinAngle);
+      double xErr2 = PS_SQR(xErr), yErr2 = PS_SQR(yErr);
+      double errScale = scale / 3600.0;
+      det->raErr->data.F64[row] = errScale * sqrt(cosAngle2 * xErr2 + sinAngle2 * yErr2);
+      det->decErr->data.F64[row] = errScale * sqrt(sinAngle2 * xErr2 + cosAngle2 * yErr2);
+      det->mask->data.U8[row] = 0;
+      plateScale += scale;
+      numGood++;
+    }
+    det->seeing *= ((float) plateScale) / ((float) numGood);
+
+    // Are we using numGood for anything outside of this function?
+    det->numGood = numGood;
+
+    if (isfinite(args->zp) && numGood > 0) {
+      psBinaryOp(mag, mag, "+", psScalarAlloc(args->zp, PS_TYPE_F32));
+    }
+
+    psTrace("ppMops.read", 2, "Read %ld good rows from %s\n", numGood, (const char*)name);
+  }
+
+  psTrace("ppMops.read", 1, "Done reading input detections\n");
+
+  return detections;
+}
+
+static long* sorted_idets_indexes(psVector *idets);
+static int compare (const void *a, const void *b);
+static psVector* fill_F32_vector(psVector* in, int size, long* indexes, int indexes_length);
+
+// We want the entries in table/fittedTrails to be sorted by increasing idet
+static psMetadata* sort_by_increasing_idet(psMetadata *aTable, long length, long fillLength) {
+  if (aTable == NULL) {
+    return aTable;
+  }
+  psVector *idets = psMetadataLookupVector(NULL, aTable, "IPP_IDET");
+  if (!idets) {
+    psError(PS_ERR_PROGRAMMING, true, "Failed to find IPP_IDET column");
+    return NULL;
+  }
+  bool is_sorted = true;
+  int count = 0;
+  long previous_value = -1;
+  while (is_sorted && (count<length) ){
+    is_sorted = (previous_value < idets->data.S64[count]);
+    previous_value = idets->data.S64[count];
+    count += 1;
+  }
+  if (is_sorted && (length==fillLength)) {
+    psTrace("ppMops.read", 5, "Values are already sorted by idet and the table is complete\n");
+    return aTable;
+  }
+  //otherwise sort (even if it should not happen) and fill
+  long* sortedIndexVector = sorted_idets_indexes(idets);
+  psMetadata* sortedTable = psMetadataAlloc();
+  char* column_name = "X_EXT";
+  psMetadataAddVector( sortedTable, PS_LIST_TAIL, column_name, 0, NULL, 
+		       fill_F32_vector(psMetadataLookupVector(NULL, aTable, column_name), fillLength, sortedIndexVector, idets->n) );
+  column_name = "Y_EXT";
+  psMetadataAddVector( sortedTable, PS_LIST_TAIL, column_name, 0, NULL, 
+		       fill_F32_vector(psMetadataLookupVector(NULL, aTable, column_name), fillLength, sortedIndexVector, idets->n) );
+  column_name = "X_EXT_SIG";
+  psMetadataAddVector( sortedTable, PS_LIST_TAIL, column_name, 0, NULL, 
+		       fill_F32_vector(psMetadataLookupVector(NULL, aTable, column_name), fillLength, sortedIndexVector, idets->n) );
+  column_name = "Y_EXT_SIG";
+  psMetadataAddVector( sortedTable, PS_LIST_TAIL, column_name, 0, NULL, 
+		       fill_F32_vector(psMetadataLookupVector(NULL, aTable, column_name), fillLength, sortedIndexVector, idets->n) );
+  column_name = "EXT_INST_MAG";
+  psMetadataAddVector( sortedTable, PS_LIST_TAIL, column_name, 0, NULL, 
+		       fill_F32_vector(psMetadataLookupVector(NULL, aTable, column_name), fillLength, sortedIndexVector, idets->n) );
+  column_name = "EXT_INST_MAG_SIG";
+  psMetadataAddVector( sortedTable, PS_LIST_TAIL, column_name, 0, NULL, 
+		       fill_F32_vector(psMetadataLookupVector(NULL, aTable, column_name), fillLength, sortedIndexVector, idets->n) );
+  column_name = "NPARAMS";
+  psMetadataAddVector( sortedTable, PS_LIST_TAIL, column_name, 0, NULL, 
+		       fill_F32_vector(psMetadataLookupVector(NULL, aTable, column_name), fillLength, sortedIndexVector, idets->n) );
+  column_name = "EXT_WIDTH_MAJ";
+  psMetadataAddVector( sortedTable, PS_LIST_TAIL, column_name, 0, NULL, 
+		       fill_F32_vector(psMetadataLookupVector(NULL, aTable, column_name), fillLength, sortedIndexVector, idets->n) );
+  column_name = "EXT_WIDTH_MIN";
+  psMetadataAddVector( sortedTable, PS_LIST_TAIL, column_name, 0, NULL, 
+		       fill_F32_vector(psMetadataLookupVector(NULL, aTable, column_name), fillLength, sortedIndexVector, idets->n) );
+  column_name = "EXT_THETA";
+  psMetadataAddVector( sortedTable, PS_LIST_TAIL, column_name, 0, NULL, 
+		       fill_F32_vector(psMetadataLookupVector(NULL, aTable, column_name), fillLength, sortedIndexVector, idets->n) );
+  column_name = "EXT_WIDTH_MAJ_ERR";
+  psMetadataAddVector( sortedTable, PS_LIST_TAIL, column_name, 0, NULL, 
+		       fill_F32_vector(psMetadataLookupVector(NULL, aTable, column_name), fillLength, sortedIndexVector, idets->n) );
+  column_name = "EXT_WIDTH_MIN_ERR";
+  psMetadataAddVector( sortedTable, PS_LIST_TAIL, column_name, 0, NULL, 
+		       fill_F32_vector(psMetadataLookupVector(NULL, aTable, column_name), fillLength, sortedIndexVector, idets->n) );
+  column_name = "EXT_THETA_ERR";
+  psMetadataAddVector( sortedTable, PS_LIST_TAIL, column_name, 0, NULL, 
+		       fill_F32_vector(psMetadataLookupVector(NULL, aTable, column_name), fillLength, sortedIndexVector, idets->n) );
+  return sortedTable;
+}
+
+static psVector* fill_F32_vector(psVector* in, int length, long* indexes, int indexes_length) {
+  if (in == NULL) {
+    return NULL;
+  }
+  psVector* out = psVectorAlloc(length, in->type.type);
+  int current_index = 0;
+  for (int i = 0; i < indexes_length; i++) {
+    while ( (current_index < indexes[i]) && (current_index < length) ) {
+      out->data.F32[current_index] = NAN;
+      current_index += 1;
+    }
+    out->data.F32[current_index] = in->data.F32[i];
+    current_index += 1;
+  }
+  while ( current_index < length ) {
+    out->data.F32[current_index] = NAN;
+    current_index += 1;
+  }
+  return out;
+}
+
+static int compare (const void *a, const void *b) {
+  return ( *((long*) a)-*((long*) b));
+}
+
+static long* sorted_idets_indexes(psVector *idets) {
+  long* sortedIndexes = (long*) malloc(idets->n * sizeof(long));
+  for (int i=0;i<idets->n;i++) {
+    sortedIndexes[i] = idets->data.S64[i];
+  }
+  qsort(sortedIndexes, idets->n, sizeof(long), (void *) compare);
+  return sortedIndexes;
+}