Index: /trunk/Ohana/src/relphot/src/RepairWarpMeasures.c
===================================================================
--- /trunk/Ohana/src/relphot/src/RepairWarpMeasures.c	(revision 39482)
+++ /trunk/Ohana/src/relphot/src/RepairWarpMeasures.c	(revision 39482)
@@ -0,0 +1,132 @@
+# include "relphot.h"
+
+// XXX need to load_images first and generate some lookup tables
+int RepairWarpMeasures (Catalog *catalog) {
+
+  off_t Nimage;
+  Image *image = getimages (&Nimage, NULL);
+
+  myAssert (!catalog->Nmeasure || catalog->measure, "programming error");
+
+  int NfixChipID = 0, NfixStackID = 0, NfixWarpID = 0, NfixWarpImageID = 0, NmissWarp = 0, NmissStack = 0, NbadWarp = 0;
+
+  int onePercent = catalog->Nmeasure / 100;
+
+  for (off_t j = 0; j < catalog->Nmeasure; j++) {
+    if (j % onePercent == 0) fprintf (stderr, ".");
+
+    Measure *measure = &catalog->measure[j];
+
+    // repair extID for non-warps:
+    if (isGPC1chip(measure->photcode))  {
+      double mjd = ohana_sec_to_mjd (measure->t);
+      int ccdnum = measure->photcode % 100;
+      uint64_t extID = CreatePSPSDetectionID (mjd, ccdnum, measure->detID);
+      if (extID != measure->extID) {
+	measure->extID = extID;
+	NfixChipID ++;
+      }
+      continue;
+    }
+
+    // repair extID for non-warps:
+    if (isGPC1stack(measure->photcode))  {
+      int im = getImageByID (measure->imageID);
+      if (im < 0) {
+	if (NmissStack < 10) fprintf (stderr, "missing stack exposure: %f %f : %d %d\n", measure->R, measure->D, measure->imageID, measure->photcode);
+	NmissStack ++;
+	continue;
+      }
+      uint64_t extID = CreatePSPSStackDetectionID (35, image[im].externID, measure->detID);
+      if (extID != measure->extID) {
+	measure->extID = extID;
+	NfixStackID ++;
+      }
+      continue;
+    }
+
+    // we are only going to repair warp detections
+    if (!isGPC1warp(measure->photcode)) continue;
+
+    // get the associated average value:
+    int averef = measure->averef;
+    Average *average = &catalog->average[averef];
+
+    // double-check for consistency
+    myAssert (average->objID == measure->objID, "objID mismatch: %f %f : %d %d\n", average->R, average->D, average->objID, measure->objID);
+
+    // warp coordinates to confirm warp
+    double X = measure->Xccd;
+    double Y = measure->Yccd;
+
+    // check if this detection is far from its correct location
+    double Rwrp = measure->R;
+    double Dwrp = measure->D;
+
+    double Rave = average->R;
+    double Dave = average->D;
+      
+    float csdec = cos(Dave * RAD_DEG);
+
+    // find the ra,dec displacement in arcsec:
+    double dR = 3600.0*(Rwrp - Rave)*csdec;
+    double dD = 3600.0*(Dwrp - Dave);
+
+    // skip detections which are within a small distance of the expected location
+    // NOTE: this actually works surprisingly well near the pole
+    double dPos = hypot(dR,dD);
+    if (dPos < 2.0) { 
+      int im = getImageByID (measure->imageID);
+      uint64_t extID = CreatePSPSStackDetectionID (34, image[im].externID, measure->detID);
+      if (extID != measure->extID) {
+	measure->extID = extID;
+	NfixWarpID ++;
+      }
+      continue;
+    }
+
+    // find the corrected warp ID:
+    int warpSeq = GetWarpSeq (image, measure->t, measure->photcode, Rave, Dave, X, Y);
+    if (warpSeq < 0) {
+      if (NmissWarp < 10) fprintf (stderr, "missing warp exposure: %f %f : %d %d\n", Rave, Dave, measure->t, measure->photcode);
+      NmissWarp ++;
+      continue;
+    }
+
+    // assert on coords being TAN?
+    Coords *imcoords = &image[warpSeq].coords;
+
+    double Rnew, Dnew;
+    XY_to_RD (&Rnew, &Dnew, X, Y, imcoords);
+
+    Rnew = ohana_normalize_angle_to_midpoint(Rnew, 180.0);
+    Dnew = ohana_normalize_angle_to_midpoint(Dnew,   0.0);
+
+    // find the ra,dec displacement in arcsec:
+    double dRnew = 3600.0*(Rnew - Rave)*csdec;
+    double dDnew = 3600.0*(Dnew - Dave);
+
+    // detections should now be repaired; fail if not
+    dPos = hypot(dRnew,dDnew);
+    if (dPos > 5.0) {
+      fprintf (stderr, "measurement still far from average location: %f %f vs %f %f : %d %d\n", Rave, Dave, Rnew, Dnew, measure->t, measure->photcode);
+      NbadWarp ++;
+    }
+
+    measure->imageID = image[warpSeq].imageID;
+    measure->R = Rnew;
+    measure->D = Dnew;
+
+    uint64_t extID = CreatePSPSStackDetectionID (34, image[warpSeq].externID, measure->detID);
+    if (extID != measure->extID) {
+      measure->extID = extID;
+      NfixWarpID ++;
+    }
+    NfixWarpImageID ++;
+  }
+
+  fprintf (stderr, "\n");
+  fprintf (stderr, "repaired %s : warp image ID %d : ext ID chip %d stack %d warp %d : missed warp %d stack %d : bad warp: %d\n", catalog->filename, NfixWarpImageID, NfixChipID, NfixStackID, NfixWarpID, NmissWarp, NmissStack, NbadWarp);
+
+  return (TRUE);
+}
Index: /trunk/Ohana/src/relphot/src/WarpImageMaps.c
===================================================================
--- /trunk/Ohana/src/relphot/src/WarpImageMaps.c	(revision 39482)
+++ /trunk/Ohana/src/relphot/src/WarpImageMaps.c	(revision 39482)
@@ -0,0 +1,279 @@
+# include "relphot.h"
+
+// we have a list of warp images, each with obstime & photcode.  I need to make a list of the obstime groups
+
+# define D_NOBSTIMES 10000
+
+typedef struct {
+  int obstime;
+  unsigned short photcode;
+  int Nwarps;
+  int NWARPS;
+  off_t *warps;
+  double *Rmin;
+  double *Rmax;
+  double *Dmin;
+  double *Dmax;
+  int    *onBoundary;
+} WarpGroup;
+
+WarpGroup *warpgroup = NULL;
+int Nwarpgroup = 0;
+int NWARPGROUP = 0;
+
+myIndexType *warpObstimeIndex = NULL;
+
+int FindWarpGroups (void) {
+
+  off_t Nimage;
+  Image *image = getimages (&Nimage, NULL);
+
+  int *obstimes = NULL;
+  int Nobstimes = 0;
+  int NOBSTIMES = D_NOBSTIMES;
+  ALLOCATE (obstimes, int, NOBSTIMES);
+
+  // first generate a list of obstime values
+  for (off_t i = 0; i < Nimage; i++) {
+    if (!isGPC1warp(image[i].photcode)) continue;
+
+    obstimes[Nobstimes] = image[i].tzero;
+    Nobstimes ++;
+
+    CHECK_REALLOCATE (obstimes, int, NOBSTIMES, Nobstimes, D_NOBSTIMES);
+  }
+  if (!Nobstimes) { 
+    free (obstimes);
+    return FALSE;
+  }
+
+  // sort the list:
+  isort (obstimes, Nobstimes);
+
+  // find the unique obstimes
+  int *uniqtimes = NULL;
+  int *uniqcount = NULL;
+  int Nuniqtimes = 0;
+  ALLOCATE (uniqtimes, int, Nobstimes);
+  ALLOCATE (uniqcount, int, Nobstimes);
+
+  // generate a uniq set of obstimes
+  for (int i = 0; i < Nobstimes; Nuniqtimes ++) {
+    uniqtimes[Nuniqtimes] = obstimes[i];
+    int Ndup = 0;
+    int lastValue = uniqtimes[Nuniqtimes];
+    while ((i < Nobstimes) && (obstimes[i] == lastValue)) {
+      i++;
+      Ndup ++;
+      uniqcount[Nuniqtimes] = Ndup;
+    }
+  }
+  REALLOCATE (uniqtimes, int, Nuniqtimes);
+  REALLOCATE (uniqcount, int, Nuniqtimes);
+  
+  // now I need to assign all warp images to a warpgroup
+  // I have 2 options to go from warp obstime to warpgroup:
+  // a) use a bisection lookup [may be slow]
+  // b) use an index on obstime [require ~5*3.14*1e7*4 bytes ~ 600MB for index]
+
+  // generate an index on obstime
+  warpObstimeIndex = myIndexInit ();
+  warpObstimeIndex->minID = obstimes[0];
+  warpObstimeIndex->maxID = obstimes[Nobstimes-1];
+
+  myIndexSetRange (warpObstimeIndex);
+
+  // generate the warp groups (create the index as we go)
+  Nwarpgroup = Nuniqtimes;
+  ALLOCATE (warpgroup, WarpGroup, Nwarpgroup);
+  
+  for (int i = 0; i < Nuniqtimes; i++) {
+    warpgroup[i].obstime = uniqtimes[i];
+    warpgroup[i].photcode = 0; // not yet set
+    warpgroup[i].Nwarps = 0;
+    warpgroup[i].NWARPS = 100;
+    ALLOCATE (warpgroup[i].warps, off_t, warpgroup[i].NWARPS);
+    myIndexSetEntry (warpObstimeIndex, uniqtimes[i], i);
+  }
+
+  // assign all warps to one of the warpgroups
+  for (off_t i = 0; i < Nimage; i++) {
+    if (!isGPC1warp(image[i].photcode)) continue;
+
+    int seq = myIndexGetEntry (warpObstimeIndex, image[i].tzero);
+    myAssert (warpgroup[seq].obstime == image[i].tzero, "oops");
+
+    if (!warpgroup[seq].photcode) {
+      warpgroup[seq].photcode = image[i].photcode;
+    } else {
+      myAssert (warpgroup[seq].photcode == image[i].photcode, "oops");
+    }
+
+    int N = warpgroup[seq].Nwarps;
+    warpgroup[seq].warps[N] = i;
+    warpgroup[seq].Nwarps ++;
+    
+    CHECK_REALLOCATE (warpgroup[seq].warps, off_t, warpgroup[seq].NWARPS, warpgroup[seq].Nwarps, 100);
+  }
+
+  // we now have the warps assigned to groups.  now we need to generate the grid data to find the warp assignments
+  for (int i = 0; i < Nwarpgroup; i++) {
+
+    myAssert (warpgroup[i].Nwarps == uniqcount[i], "failure");
+    ALLOCATE (warpgroup[i].Rmin, double, warpgroup[i].Nwarps);
+    ALLOCATE (warpgroup[i].Rmax, double, warpgroup[i].Nwarps);
+    ALLOCATE (warpgroup[i].Dmin, double, warpgroup[i].Nwarps);
+    ALLOCATE (warpgroup[i].Dmax, double, warpgroup[i].Nwarps);
+    ALLOCATE (warpgroup[i].onBoundary, int, warpgroup[i].Nwarps);
+
+    for (int j = 0; j < warpgroup[i].Nwarps; j++) {
+
+      off_t N = warpgroup[i].warps[j];
+      
+      int Nx = image[N].NX;
+      int Ny = image[N].NY;
+
+      warpgroup[i].onBoundary[j] = FALSE;
+
+      double R, D;
+      XY_to_RD (&R, &D, 0.0, 0.0, &image[N].coords);
+      R = ohana_normalize_angle_to_midpoint(R, 180.0);
+      warpgroup[i].Rmin[j] = R;
+      warpgroup[i].Rmax[j] = R;
+      warpgroup[i].Dmin[j] = D;
+      warpgroup[i].Dmax[j] = D;
+
+      // XXX need to worry about 0,360 boundary
+      XY_to_RD (&R, &D, Nx, 0.0, &image[N].coords);
+      R = ohana_normalize_angle_to_midpoint(R, 180.0);
+      warpgroup[i].Rmin[j] = MIN (warpgroup[i].Rmin[j], R);
+      warpgroup[i].Rmax[j] = MAX (warpgroup[i].Rmax[j], R);
+      warpgroup[i].Dmin[j] = MIN (warpgroup[i].Dmin[j], D);
+      warpgroup[i].Dmax[j] = MAX (warpgroup[i].Dmax[j], D);
+
+      // XXX need to worry about 0,360 boundary
+      XY_to_RD (&R, &D, 0.0, Ny, &image[N].coords);
+      R = ohana_normalize_angle_to_midpoint(R, 180.0);
+      warpgroup[i].Rmin[j] = MIN (warpgroup[i].Rmin[j], R);
+      warpgroup[i].Rmax[j] = MAX (warpgroup[i].Rmax[j], R);
+      warpgroup[i].Dmin[j] = MIN (warpgroup[i].Dmin[j], D);
+      warpgroup[i].Dmax[j] = MAX (warpgroup[i].Dmax[j], D);
+
+      // XXX need to worry about 0,360 boundary
+      XY_to_RD (&R, &D, Nx, Ny, &image[N].coords);
+      R = ohana_normalize_angle_to_midpoint(R, 180.0);
+      warpgroup[i].Rmin[j] = MIN (warpgroup[i].Rmin[j], R);
+      warpgroup[i].Rmax[j] = MAX (warpgroup[i].Rmax[j], R);
+      warpgroup[i].Dmin[j] = MIN (warpgroup[i].Dmin[j], D);
+      warpgroup[i].Dmax[j] = MAX (warpgroup[i].Dmax[j], D);
+
+      // bump Rmin,Rmax,Dmin,Dmax 10 arcsec worth of padding
+      double dR = 10.0/3600.0 / cos (RAD_DEG*warpgroup[i].Dmin[j]);
+      double dD = 10.0/3600.0;
+
+      // at north pole, force test of all nearby skycells
+      if (warpgroup[i].Dmax[j] > 89.8) {
+	warpgroup[i].Rmin[j] =   0.0;
+	warpgroup[i].Rmax[j] = 360.0;
+	warpgroup[i].Dmax[j] =  90.0;
+	warpgroup[i].Dmin[j] -=  dD;
+	continue;
+      } 
+
+      if (warpgroup[i].Rmax[j] - warpgroup[i].Rmin[j] > 270.0) {
+	// Rmin and Rmax are in the range 0 - 360.  For images at the 0,360 boundary,
+	// "Rmax" is the lower edge, and "Rmin" is the upper edge.  we need to flip them 
+	// and then break the 0-360 range:
+	double tmp = warpgroup[i].Rmin[j];
+	warpgroup[i].Rmin[j] = warpgroup[i].Rmax[j] - 360.0;
+	warpgroup[i].Rmax[j] = tmp;
+	warpgroup[i].onBoundary[j] = TRUE;
+      }
+
+      warpgroup[i].Rmin[j] -=  dR;
+      warpgroup[i].Rmax[j] +=  dR;
+      warpgroup[i].Dmin[j] -=  dD;
+      warpgroup[i].Dmax[j] +=  dD;
+    }
+  }
+
+  free (uniqtimes); 
+  free (uniqcount); 
+  free (obstimes);
+  return TRUE;
+}
+
+int GetWarpSeq (Image *image, int obstime, unsigned short photcode, double Rave, double Dave, float X, float Y) {
+
+  if (!warpObstimeIndex) return -1;
+
+  int seq = myIndexGetEntry (warpObstimeIndex, obstime);
+  myAssert (seq > -1, "ooops");
+  myAssert (warpgroup[seq].photcode == photcode, "oops");
+  
+  double dPosMin = NAN;
+  int nPosMin = -1;
+
+  // we now have the warp group, but which is the correct warp?
+  for (int i = 0; i < warpgroup[seq].Nwarps; i++) {
+    if (warpgroup[seq].onBoundary[i]) {
+      int inRange1 = (Rave >= warpgroup[seq].Rmin[i]) && (Rave <= warpgroup[seq].Rmax[i]);
+      int inRange2 = (Rave >= warpgroup[seq].Rmin[i] + 360.0) && (Rave <= warpgroup[seq].Rmax[i] + 360.0);
+      if (!inRange1 && !inRange2) continue;
+    } else {
+      if (Rave < warpgroup[seq].Rmin[i]) continue;
+      if (Rave > warpgroup[seq].Rmax[i]) continue;
+    }
+    if (Dave < warpgroup[seq].Dmin[i]) continue;
+    if (Dave > warpgroup[seq].Dmax[i]) continue;
+
+    // this is a possible image: check the coordinates:
+    int N = warpgroup[seq].warps[i];
+
+    // project Rave,Dave to image pixels
+    double Xtst, Ytst;
+    RD_to_XY (&Xtst, &Ytst, Rave, Dave, &image[N].coords);
+
+    // find the pixel offset
+    double dX = (Xtst - X);
+    double dY = (Ytst - Y);
+    
+    // skip detections which are within a small distance from the expected location
+    double dPos = hypot(dX,dY);
+    if (dPos < 20.0) return N; // 20.0 pixels = 5.0 arcsec
+
+    // check for multiple possible matches??
+
+    // record the min pos so we can see how things fail
+    if (isnan(dPosMin)) {
+      dPosMin = dPos;
+      nPosMin = i;
+    } else {
+      if (dPos < dPosMin) {
+	dPosMin = dPos;
+	nPosMin = i;
+      }
+    }
+  }
+  
+  // if we did not find any matched, give up
+  if (nPosMin < 0) return -1;
+
+  // return closest image
+  return warpgroup[seq].warps[nPosMin];
+}
+
+void FreeWarpGroups (void) {
+
+  for (int i = 0; i < Nwarpgroup; i++) {
+    free (warpgroup[i].warps);
+    free (warpgroup[i].Rmin);
+    free (warpgroup[i].Rmax);
+    free (warpgroup[i].Dmin);
+    free (warpgroup[i].Dmax);
+    free (warpgroup[i].onBoundary);
+  }    
+
+  free (warpgroup);
+  myIndexFree (warpObstimeIndex);
+}
Index: /trunk/Ohana/src/relphot/src/myIndex.c
===================================================================
--- /trunk/Ohana/src/relphot/src/myIndex.c	(revision 39482)
+++ /trunk/Ohana/src/relphot/src/myIndex.c	(revision 39482)
@@ -0,0 +1,77 @@
+# include "relphot.h"
+
+myIndexType *myIndexInit () {
+
+  myIndexType *myIndex;
+  ALLOCATE (myIndex, myIndexType, 1);
+
+  myIndex->minID = MAX_INT;
+  myIndex->maxID = 0;
+
+  myIndex->NINDEX = 1000;
+  myIndex->Nindex = 0;
+
+  ALLOCATE (myIndex->index, int, myIndex->NINDEX);
+
+  return myIndex;
+}
+
+int myIndexFree (myIndexType *myIndex) {
+
+  if (!myIndex) return TRUE;
+  free (myIndex->index);
+  free (myIndex);
+  return TRUE;
+}
+
+// set minID and maxID externally
+int myIndexUpdateLimits (myIndexType *myIndex, int value) {
+
+  myIndex->maxID = MAX(value, myIndex->maxID);
+  myIndex->minID = MIN(value, myIndex->minID);
+
+  return TRUE;
+}
+
+// once minID and maxID are set, this function defines the range and inits
+int myIndexSetRange (myIndexType *myIndex) {
+
+  myIndex->Nindex = myIndex->maxID - myIndex->minID + 1;
+  myIndex->NINDEX = myIndex->Nindex;
+  
+  REALLOCATE (myIndex->index, int, myIndex->NINDEX);
+
+  int i;
+  for (i = 0; i < myIndex->Nindex; i++) {
+    myIndex->index[i] = -1;
+  }
+
+  return TRUE;
+}
+
+// once minID and maxID are set, this function defines the range and inits
+int myIndexSetEntry (myIndexType *myIndex, int value, int entry) {
+
+  int n = value - myIndex->minID;
+
+  myAssert (n >= 0, "impossible!");
+
+  if (myIndex->index[n] != -1) {
+    fprintf (stderr, "skipping duplicate ID value: %d\n", value);
+    return -1;
+  }
+
+  myIndex->index[n] = entry;
+  return n;
+}
+
+// once minID and maxID are set, this function defines the range and inits
+int myIndexGetEntry (myIndexType *myIndex, int value) {
+
+  if (value < myIndex->minID) return -1; // not in this index
+  if (value > myIndex->maxID) return -1; // not in this index
+
+  int n = value - myIndex->minID;
+
+  return myIndex->index[n];
+}
Index: /trunk/Ohana/src/relphot/src/psps_ids.c
===================================================================
--- /trunk/Ohana/src/relphot/src/psps_ids.c	(revision 39482)
+++ /trunk/Ohana/src/relphot/src/psps_ids.c	(revision 39482)
@@ -0,0 +1,43 @@
+# include "relphot.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 CreatePSPSStackDetectionID(int sourceID, int imageID, int detID) {
+  // sourceID : ID of database + table that tracked the image (< 0x7f = 127)
+  // imageID : external ID of the image which provided the detections (< 0x1000.0000 ~ 2.7e8)
+  // detID : detection sequence in image (< 0x1000.0000 ~ 2.7e8)
+
+  assert (detID    < 0x10000000);
+  assert (imageID  < 0x10000000);
+  assert (sourceID < 0x7f);
+  
+  uint64_t detectid = ((uint64_t)sourceID << 56) + ((uint64_t)imageID << 28) + (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;
+}