Index: /branches/eam_branches/20091201/Ohana/src/addstar/Makefile
===================================================================
--- /branches/eam_branches/20091201/Ohana/src/addstar/Makefile	(revision 26710)
+++ /branches/eam_branches/20091201/Ohana/src/addstar/Makefile	(revision 26711)
@@ -81,5 +81,6 @@
 $(SRC)/ReadImageHeader.$(ARCH).o \
 $(SRC)/UpdateImageIDs.$(ARCH).o \
-$(SRC)/update_coords.$(ARCH).o
+$(SRC)/update_coords.$(ARCH).o \
+$(SRC)/psps_ids.$(ARCH).o
 
 ADDSTARD = \
@@ -118,5 +119,6 @@
 $(SRC)/opening_angle.$(ARCH).o \
 $(SRC)/replace_match.$(ARCH).o \
-$(SRC)/update_coords.$(ARCH).o
+$(SRC)/update_coords.$(ARCH).o \
+$(SRC)/psps_ids.$(ARCH).o
 
 ADDSTART = \
@@ -157,5 +159,6 @@
 $(SRC)/NewImage_Thread.$(ARCH).o \
 $(SRC)/NewReflist_Thread.$(ARCH).o \
-$(SRC)/NewRefcat_Thread.$(ARCH).o 
+$(SRC)/NewRefcat_Thread.$(ARCH).o \
+$(SRC)/psps_ids.$(ARCH).o
 
 ADDSTARC = \
@@ -206,5 +209,6 @@
 $(SRC)/ConfigInit.$(ARCH).o \
 $(SRC)/Shutdown.$(ARCH).o \
-$(SRC)/SetSignals.$(ARCH).o
+$(SRC)/SetSignals.$(ARCH).o \
+$(SRC)/psps_ids.$(ARCH).o
 
 SEDSTAR = \
Index: /branches/eam_branches/20091201/Ohana/src/addstar/include/addstar.h
===================================================================
--- /branches/eam_branches/20091201/Ohana/src/addstar/include/addstar.h	(revision 26710)
+++ /branches/eam_branches/20091201/Ohana/src/addstar/include/addstar.h	(revision 26711)
@@ -289,4 +289,8 @@
 PhotCode *pmm_get_photcode (char *emulsion, char *filter);
 
+#define PSPS_ID TRUE
+uint64_t CreatePSPSDetectionID(double tobs, int ccdid, int detID);
+uint64_t CreatePSPSObjectID(double ra, double dec);
+
 // this is a gnu extension?? caution!
 void *memrchr(const void *s, int c, size_t n);
Index: /branches/eam_branches/20091201/Ohana/src/addstar/src/FilterStars.c
===================================================================
--- /branches/eam_branches/20091201/Ohana/src/addstar/src/FilterStars.c	(revision 26710)
+++ /branches/eam_branches/20091201/Ohana/src/addstar/src/FilterStars.c	(revision 26711)
@@ -1,4 +1,4 @@
 # include "addstar.h"
-# define EXTERNAL_ID FALSE
+# define EXTERNAL_ID TRUE
 
 /* exclude some detections based on various filter options:
@@ -88,14 +88,13 @@
     }
 
-# if (0) 
     if (PSPS_ID) {
       double mjd;
       mjd = ohana_sec_to_mjd (image[0].tzero);
-      // XXX this is clearly wrong : what does PSPS want?
-      stars[N].measure.extID = PSPS_create_detectid (mjd, imageID);
+      stars[N].measure.extID = CreatePSPSDetectionID(mjd, image[0].ccdnum, stars[N].measure.detID);
+    } else {
+      stars[N].measure.extID = 0;
     }
-# endif
 
-    stars[N].measure.imageID = imageID; // this value is update in UpdateImageIDs
+    stars[N].measure.imageID = imageID; // this value is updated in UpdateImageIDs
 
     N ++;
Index: /branches/eam_branches/20091201/Ohana/src/addstar/src/LoadData.c
===================================================================
--- /branches/eam_branches/20091201/Ohana/src/addstar/src/LoadData.c	(revision 26710)
+++ /branches/eam_branches/20091201/Ohana/src/addstar/src/LoadData.c	(revision 26711)
@@ -10,4 +10,5 @@
   char *name;
   int i, j, Nvalid, Nhead, Ndata, Nskip;
+  uint32_t parentID = UINT32_MAX;
   Stars *inStars;
 
@@ -58,8 +59,14 @@
     // skip the table if there is no data segment (eg, mosaic WRP image)
     if (!strcmp(headerSets[i].extdata, "NONE")) {
+      if (!strcmp(headerSets[i].exthead, "PHU") && (Nimages > 1)) {
+        // This image is the parent of subsequent images
+        parentID = Nvalid;
+        images[0][Nvalid].parentID = UINT32_MAX;
+      }
       Nvalid++;
       CHECK_REALLOCATE (images[0], Image, NVALID, Nvalid, 10);
       continue;
     }
+    images[0][Nvalid].parentID = parentID;
 
     // advance the pointer to the start of the corresponding table block
Index: anches/eam_branches/20091201/Ohana/src/addstar/src/PSPS_IDs.c
===================================================================
--- /branches/eam_branches/20091201/Ohana/src/addstar/src/PSPS_IDs.c	(revision 26710)
+++ 	(revision )
@@ -1,71 +1,0 @@
-# include "addstar.h"
-
-uint64_7 PSPS_create_detectID (double MJDobs, int imageID) {
-
-  // MJD has a PS lifetime range of : 3288 (2009/01/01) to 9132 (2025/01/01) (< 2^14)
-
-  static double t0 = 51544.5; // MJD @ 2000/01/01,12:00:00
-
-  uint64_t iTime;
-  uint64_t detectID;
-
-  iTime = floor(tobs - t0);
-
-  // iTime takes the upper 16 bits, imageID can take up to lower 48 bits
-  detectID = (iTime << 48) | imageID;
-
-  return detectID;
-}
-
-uint64_t PSPS_create_objID (double ra, double dec) {
-
-  static double zh = 0.0083333; // zone height = 30 arcsec
-
-  uint64_t objID;
-  uint64_t part1, part2, part3;
-
-  double zid, zresid;
-  int izone;
-
-  zid = (dec + 90.0) / zh; // 0 - 180*60*2 = 21600 (< 15 bits)
-  izone = (int) floor(zid);
-  zresid = zid -  ((float) izone); // 0 - 1.0 
-  part1 = (uint64_t)( izone  * 10000000000000) ;
-  part2 = ((uint64_t)(ra * 1 000 000.)) * 10000 ; // 0 - 360.0*1e6 = 3.6e8 (< 29 bits)
-  part3 = (int) (zresid * 10000.0) ; // 0 - 10000 (1 bit == 30/10000 arcsec = 0.003 arcsec) (< 14 bits)
-
-  objID = part1 + part2 + part3;
-
-  return objID;
-}
-
-uint64_t PSPS_create_objID_eam (double ra, double dec) {
-
-  static double zh = 0.0083333; // zone height = 30 arcsec
-
-  uint64_t objID;
-  uint64_t part1, part2, part3;
-
-  double zid, zresid;
-  int izone, ira;
-
-  zid = (dec + 90.0) / zh; // 0 - 180*60*2 = 21600 (< 15 bits)
-  izone = (int) floor(zid);
-  zresid = zid -  ((float) izone); // 0 - 1.0 
-
-  ira = ra * 1 000 000; // 0 - 360.0*1e6 = 3.6e8 (< 29 bits) 
-
-  part1 = (izone << 48);
-  part2 = (ira << 16);
-  part3 = (int) (zresid * 10000.0) ; // 0 - 10000 (1 bit == 30/10000 arcsec = 0.003 arcsec) (< 14 bits)
-
-  objID = part1 | part2 | part3;
-
-  return objID;
-}
-
-// part1 : 15 bits
-// part2 : 29 bits
-// part3 : 14 bits
-// 15 + 14 + 29 bits = 58 bits...
-
Index: /branches/eam_branches/20091201/Ohana/src/addstar/src/ReadImageHeader.c
===================================================================
--- /branches/eam_branches/20091201/Ohana/src/addstar/src/ReadImageHeader.c	(revision 26710)
+++ /branches/eam_branches/20091201/Ohana/src/addstar/src/ReadImageHeader.c	(revision 26711)
@@ -5,7 +5,7 @@
 int ReadImageHeader (Header *header, Image *image, int photcode) {
 
-  int Nastro, ccdnum, hour, min, Nx, Ny, haveNx, haveNy, sourceID;
+  int Nastro, hour, min, Nx, Ny, haveNx, haveNy, sourceID;
   double tmp, sec, Cerror, ZeroPt, FWHM_X, FWHM_Y;
-  char *c, photname[64], line[80];
+  char *c, photname[64], line[80], ccdnum[64];
   PhotCode *photcodeData = NULL;
 
@@ -174,8 +174,19 @@
   } 
 
-  if (!gfits_scan (header, CCDNumKeyword, "%d", 1, &ccdnum)) {
-    image[0].ccdnum = 0xff;
-  } else {
-    image[0].ccdnum = ccdnum;
+  if (!gfits_scan (header, CCDNumKeyword, "%s", 1, ccdnum)) {
+    fprintf (stderr, "CCDNumKeyword %s not found\n", CCDNumKeyword);
+    return FALSE;
+  } else {
+    // CCDNumKeyword (EXTNAME) is a string yet we need an integer
+    // RULE: ccdnum is the value of the first string of digits in CCDNumKeyword
+    // For example: ccdnum(XY42.hdr) = 42, ccdnum(XY01.hdr) = 1
+    char *p = ccdnum;
+    // get rid of any leading non-digit characters
+    while (*p && !isdigit(*p)) p++;
+    if (*p == 0) {
+      fprintf (stderr, "invalid ccdnum %s\n", ccdnum);
+      return FALSE;
+    }
+    image[0].ccdnum = atoi(p);
   }
 
Index: /branches/eam_branches/20091201/Ohana/src/addstar/src/UpdateImageIDs.c
===================================================================
--- /branches/eam_branches/20091201/Ohana/src/addstar/src/UpdateImageIDs.c	(revision 26710)
+++ /branches/eam_branches/20091201/Ohana/src/addstar/src/UpdateImageIDs.c	(revision 26711)
@@ -40,4 +40,9 @@
   for (i = 0; i < Nimages; i++) {
     images[i].imageID += imageID;
+    if (images[i].parentID == UINT32_MAX) {
+        images[i].parentID = 0;
+    } else {
+        images[i].parentID += imageID;
+    }
   }
 
Index: /branches/eam_branches/20091201/Ohana/src/addstar/src/build_links.c
===================================================================
--- /branches/eam_branches/20091201/Ohana/src/addstar/src/build_links.c	(revision 26710)
+++ /branches/eam_branches/20091201/Ohana/src/addstar/src/build_links.c	(revision 26711)
@@ -1,3 +1,26 @@
 # include "addstar.h"
+
+/* 
+
+There are two modes for the measure table: sorted and unsorted. 
+
+In sorted mode, all measures associated with a given average are in a single block.  The block
+is pointed to by average->measureOffset and the range is average->Nmeasure
+
+In unsorted mode, it is not possible to go directly from average to measure without scanning.
+In this case, it is necessary to use the value measure->averef to find the corresponding
+average entry.  
+
+Note that average->measureOffset and measure->averef are only valid for a given load of the
+data: they refer to the sequence number in the data blocks.
+
+next_meas is a list of the equivalent sequence of the measure block as if it were sorted.
+
+to find the sequence of measurements for a given average:
+n_0 = average->measureOffset
+n_1 = next_meas[n_0]
+n_i = next_meas[n_i-1]
+
+*/
 
 /* build the initial links assuming the table is sorted, 
@@ -31,4 +54,110 @@
 }
 
+/* construct measure links which are valid FOR THIS LOAD
+ * - if we have a full load, we will get links which can
+ *   be used by other programs (eg, relphot, etc)
+ * - if we have a partial load, the links are only valid
+ *   for that partial load
+ */ 
+
+int *build_measure_links (Average *average, int Naverage, Measure *measure, int Nmeasure) {
+
+  int i, m, k, Nm, averef;
+  int *next;
+
+  ALLOCATE (next, int, Nmeasure);
+
+  /* reset the Nm, offset values for average */
+  for (i = 0; i < Naverage; i++) {
+    average[i].measureOffset = -1;
+    average[i].Nmeasure     =  0;
+  }
+
+  for (Nm = 0; Nm < Nmeasure; Nm++) {
+    averef = measure[Nm].averef;
+    m = average[averef].measureOffset;  
+    next[Nm] = -1;
+
+    if (m == -1) { /* no links yet for source */
+      average[averef].measureOffset = Nm;
+      average[averef].Nmeasure     = 1;
+      continue;
+    }
+
+    for (k = 0; next[m] != -1; k++) {
+      m = next[m];
+      if (m >= Nmeasure) {
+	fprintf (stderr, "WARNING: m out of bounds (1)\n");
+      }
+    }
+
+    average[averef].Nmeasure = k + 2;
+    next[m] = Nm;
+    if (m >= Nmeasure) {
+      fprintf (stderr, "WARNING: m out of bounds (2)\n");
+    }
+  }
+  return (next);
+}
+
+/* average[].measureOffset, average[].Nmeasure are valid within an addstar run */
+int add_meas_link (Average *average, int *next, int Nmeasure, int NMEASURE) {
+
+  int k, m;
+
+  /* if we have trouble, check validity of next[m] : m < Nmeasure */
+  m = average[0].measureOffset;  
+
+  for (k = 0; k < average[0].Nmeasure - 1; k++)  {
+    m = next[m];
+    if (m >= NMEASURE) {
+      fprintf (stderr, "WARNING: m out of bounds (3)\n");
+    }
+  }
+
+  /* set up references */
+  next[Nmeasure] = -1;
+  if (Nmeasure >= NMEASURE) {
+    fprintf (stderr, "WARNING: Nmeasure out of bounds (1)\n");
+  }
+
+  if (m == -1) {
+    average[0].measureOffset = Nmeasure;
+  } else {
+    next[m] = Nmeasure;
+    if (m >= NMEASURE) {
+      fprintf (stderr, "WARNING: m out of bounds (4)\n");
+    }
+  }
+
+  return (TRUE);
+}
+
+/* Missing does not carry enough information to reconstruct the links
+   we must always save the missing table, if it exists */
+
+Measure *sort_measure (Average *average, int Naverage, Measure *measure, int Nmeasure, int *next) {
+
+  int i, k, n, N;
+  Measure *tmpmeasure;
+
+  /* fix order of Measure (memory intensive, but fast) */
+  N = 0; 
+  ALLOCATE (tmpmeasure, Measure, Nmeasure);
+  for (i = 0; i < Naverage; i++) {
+    n = average[i].measureOffset;
+    average[i].measureOffset = N;
+    for (k = 0; k < average[i].Nmeasure; k++, N++) {
+      if (n == -1) abort();
+      tmpmeasure[N] = measure[n]; 
+      if (measure[n].averef != i) abort();
+      tmpmeasure[N].averef = i;
+      n = next[n];
+    }
+  }
+  free (measure);
+  return (tmpmeasure);
+}
+
 /* build the initial links assuming the table is sorted */
 int *init_missing_links (Average *average, int Naverage, Missing *missing, int Nmissing) {
@@ -57,37 +186,4 @@
 }
 
-/* average[].measureOffset, average[].Nmeasure are valid within an addstar run */
-int add_meas_link (Average *average, int *next, int Nmeasure, int NMEASURE) {
-
-  int k, m;
-
-  /* if we have trouble, check validity of next[m] : m < Nmeasure */
-  m = average[0].measureOffset;  
-
-  for (k = 0; k < average[0].Nmeasure - 1; k++)  {
-    m = next[m];
-    if (m >= NMEASURE) {
-      fprintf (stderr, "WARNING: m out of bounds (3)\n");
-    }
-  }
-
-  /* set up references */
-  next[Nmeasure] = -1;
-  if (Nmeasure >= NMEASURE) {
-    fprintf (stderr, "WARNING: Nmeasure out of bounds (1)\n");
-  }
-
-  if (m == -1) {
-    average[0].measureOffset = Nmeasure;
-  } else {
-    next[m] = Nmeasure;
-    if (m >= NMEASURE) {
-      fprintf (stderr, "WARNING: m out of bounds (4)\n");
-    }
-  }
-
-  return (TRUE);
-}
-
 int add_miss_link (Average *average, int *next, int Nmissing) {
 
@@ -109,75 +205,4 @@
 }
 
-/* construct measure links which are valid FOR THIS LOAD
- * - if we have a full load, we will get links which can
- *   be used by other programs (eg, relphot, etc)
- * - if we have a partial load, the links are only valid
- *   for that partial load
- */ 
-
-int *build_measure_links (Average *average, int Naverage, Measure *measure, int Nmeasure) {
-
-  int i, m, k, Nm, averef;
-  int *next;
-
-  ALLOCATE (next, int, Nmeasure);
-
-  /* reset the Nm, offset values for average */
-  for (i = 0; i < Naverage; i++) {
-    average[i].measureOffset = -1;
-    average[i].Nmeasure     =  0;
-  }
-
-  for (Nm = 0; Nm < Nmeasure; Nm++) {
-    averef = measure[Nm].averef;
-    m = average[averef].measureOffset;  
-    next[Nm] = -1;
-
-    if (m == -1) { /* no links yet for source */
-      average[averef].measureOffset = Nm;
-      average[averef].Nmeasure     = 1;
-      continue;
-    }
-
-    for (k = 0; next[m] != -1; k++) {
-      m = next[m];
-      if (m >= Nmeasure) {
-	fprintf (stderr, "WARNING: m out of bounds (1)\n");
-      }
-    }
-
-    average[averef].Nmeasure = k + 2;
-    next[m] = Nm;
-    if (m >= Nmeasure) {
-      fprintf (stderr, "WARNING: m out of bounds (2)\n");
-    }
-  }
-  return (next);
-}
-
-/* Missing does not carry enough information to reconstruct the links
-   we must always save the missing table, if it exists */
-
-Measure *sort_measure (Average *average, int Naverage, Measure *measure, int Nmeasure, int *next) {
-
-  int i, k, n, N;
-  Measure *tmpmeasure;
-
-  /* fix order of Measure (memory intensive, but fast) */
-  N = 0; 
-  ALLOCATE (tmpmeasure, Measure, Nmeasure);
-  for (i = 0; i < Naverage; i++) {
-    n = average[i].measureOffset;
-    average[i].measureOffset = N;
-    for (k = 0; k < average[i].Nmeasure; k++, N++) {
-      tmpmeasure[N] = measure[n]; 
-      tmpmeasure[N].averef = i;
-      n = next[n];
-    }
-  }
-  free (measure);
-  return (tmpmeasure);
-}
-
 Missing *sort_missing (Average *average, int Naverage, Missing *missing, int Nmissing, int *next) {
 
@@ -199,2 +224,3 @@
   return (tmpmissing);
 }
+
Index: /branches/eam_branches/20091201/Ohana/src/addstar/src/find_matches.c
===================================================================
--- /branches/eam_branches/20091201/Ohana/src/addstar/src/find_matches.c	(revision 26710)
+++ /branches/eam_branches/20091201/Ohana/src/addstar/src/find_matches.c	(revision 26711)
@@ -271,4 +271,10 @@
     catalog[0].average[Nave].catID     	   = catID;
     catalog[0].average[Nave].flags         = 0;
+    if (PSPS_ID) {
+        catalog[0].average[Nave].extID = CreatePSPSObjectID(catalog[0].average[Nave].R,
+                                                            catalog[0].average[Nave].D);
+    } else {
+        catalog[0].average[Nave].extID         = 0;
+    }
 
     objID ++;
Index: /branches/eam_branches/20091201/Ohana/src/addstar/src/find_matches_closest.c
===================================================================
--- /branches/eam_branches/20091201/Ohana/src/addstar/src/find_matches_closest.c	(revision 26710)
+++ /branches/eam_branches/20091201/Ohana/src/addstar/src/find_matches_closest.c	(revision 26711)
@@ -266,4 +266,11 @@
     catalog[0].average[Nave].catID     	   = catID;
     catalog[0].average[Nave].flags         = 0;
+    if (PSPS_ID) {
+        catalog[0].average[Nave].extID = CreatePSPSObjectID(catalog[0].average[Nave].R,
+                                                            catalog[0].average[Nave].D);
+    } else {
+        catalog[0].average[Nave].extID         = 0;
+    }
+
 
     objID ++;
Index: /branches/eam_branches/20091201/Ohana/src/addstar/src/find_matches_closest_refstars.c
===================================================================
--- /branches/eam_branches/20091201/Ohana/src/addstar/src/find_matches_closest_refstars.c	(revision 26710)
+++ /branches/eam_branches/20091201/Ohana/src/addstar/src/find_matches_closest_refstars.c	(revision 26711)
@@ -291,4 +291,10 @@
     catalog[0].average[Nave].catID     	   = catID;
     catalog[0].average[Nave].flags     	   = 0;
+    if (PSPS_ID) {
+        catalog[0].average[Nave].extID = CreatePSPSObjectID(catalog[0].average[Nave].R,
+                                                            catalog[0].average[Nave].D);
+    } else {
+        catalog[0].average[Nave].extID         = 0;
+    }
 
     objID ++;
Index: /branches/eam_branches/20091201/Ohana/src/addstar/src/find_matches_refstars.c
===================================================================
--- /branches/eam_branches/20091201/Ohana/src/addstar/src/find_matches_refstars.c	(revision 26710)
+++ /branches/eam_branches/20091201/Ohana/src/addstar/src/find_matches_refstars.c	(revision 26711)
@@ -261,4 +261,11 @@
     catalog[0].average[Nave].catID     	   = catID;
     catalog[0].average[Nave].flags     	   = 0;
+    if (PSPS_ID) {
+        catalog[0].average[Nave].extID = CreatePSPSObjectID(catalog[0].average[Nave].R,
+                                                            catalog[0].average[Nave].D);
+    } else {
+        catalog[0].average[Nave].extID         = 0;
+    }
+
 
     objID ++;
@@ -290,6 +297,10 @@
 
       stars[N][0].found = Nmeas;
-      next_meas[Nmeas] = -1;
-      Nmeas ++;
+      next_meas[Nmeas] = -1;  // initial value here, update below
+      Nmeas ++; 
+    }
+
+    for (j = 0; j < NREFSTAR_GROUP - 1; j++) {
+	next_meas[Nmeas - NREFSTAR_GROUP + j] = Nmeas - NREFSTAR_GROUP + j + 1;
     }
     Nave ++;
Index: /branches/eam_branches/20091201/Ohana/src/addstar/src/psps_ids.c
===================================================================
--- /branches/eam_branches/20091201/Ohana/src/addstar/src/psps_ids.c	(revision 26711)
+++ /branches/eam_branches/20091201/Ohana/src/addstar/src/psps_ids.c	(revision 26711)
@@ -0,0 +1,34 @@
+# include "addstar.h"
+
+// Compute IDs for PSPS based on recipe they defined
+
+uint64_t
+CreatePSPSDetectionID(double tobs, int ccdid, int detID)
+{
+    // t0 for detection id is 2007-01-01 00:00:00 utc
+    double  t0 = 54101.0;
+    double diff = floor( 100000. * (tobs - t0) );
+    int itmp = diff;
+
+    // ccdid must be < 100
+    uint64_t detectid = 1000000000*((uint64_t) itmp) + 10000000 * ((uint64_t) ccdid) +
+             ((uint64_t) detID);
+
+    return detectid;
+}
+    
+uint64_t
+CreatePSPSObjectID(double ra, double dec)
+{
+    double zh = 0.0083333;
+    double zid = (dec + 90.) / zh;             // 0 - 180*60*2 = 21600 < 15 bits
+    int izone = (int) floor(zid);
+    double zresid = zid -  ((float) izone);    // 0.0 - 1.0
+
+    uint64_t part1, part2, part3;
+    part1 = (uint64_t)( izone  * 10000000000000LL) ; 
+    part2 = ((uint64_t)(ra * 1000000.)) * 10000 ; // 0 - 360*1e6 = 3.6e8 (< 29 bits)
+    part3 = (int) (zresid * 10000.0) ; // 0 - 10000 (1 bit == 30/10000 arcsec) (< 14 bits)
+
+    return part1 + part2 + part3;
+}
