Index: /trunk/Ohana/src/relphot/Makefile
===================================================================
--- /trunk/Ohana/src/relphot/Makefile	(revision 39512)
+++ /trunk/Ohana/src/relphot/Makefile	(revision 39513)
@@ -33,4 +33,5 @@
 $(SRC)/StarOps.$(ARCH).o	 \
 $(SRC)/WarpImageMaps.$(ARCH).o	 \
+$(SRC)/StackImageMaps.$(ARCH).o	 \
 $(SRC)/RepairWarpMeasures.$(ARCH).o	 \
 $(SRC)/myIndex.$(ARCH).o	 \
@@ -88,4 +89,5 @@
 $(SRC)/StarOps.$(ARCH).o	 \
 $(SRC)/WarpImageMaps.$(ARCH).o	 \
+$(SRC)/StackImageMaps.$(ARCH).o	 \
 $(SRC)/RepairWarpMeasures.$(ARCH).o	 \
 $(SRC)/myIndex.$(ARCH).o	 \
Index: /trunk/Ohana/src/relphot/include/relphot.h
===================================================================
--- /trunk/Ohana/src/relphot/include/relphot.h	(revision 39512)
+++ /trunk/Ohana/src/relphot/include/relphot.h	(revision 39513)
@@ -634,2 +634,7 @@
 uint64_t CreatePSPSStackDetectionID(int sourceID, int imageID, int detID);
 uint64_t CreatePSPSDetectionID(double tobs, int ccdid, int detID);
+
+int MakeStackIndex (void);
+void FreeStackGroups (void);
+int GetStackSeq (Image *image, double Rstk, double Dstk, unsigned short photcode, float X, float Y);
+void sort_by_ra (double *R, double *D, int *I, int *S, int N);
Index: /trunk/Ohana/src/relphot/src/RepairWarpMeasures.c
===================================================================
--- /trunk/Ohana/src/relphot/src/RepairWarpMeasures.c	(revision 39512)
+++ /trunk/Ohana/src/relphot/src/RepairWarpMeasures.c	(revision 39513)
@@ -18,4 +18,11 @@
     Measure *measure = &catalog->measure[j];
 
+    // 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);
+
     // repair extID for non-warps:
     if (isGPC1chip(measure->photcode))  {
@@ -34,7 +41,37 @@
       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;
+	// if we don't have a matching image, we are probably assigned to an image on the wrong side of the sky
+	// reconstruct the imageID from the R,D,photcode & X,Y consistency
+	im = GetStackSeq (image, average->R, average->D, measure->photcode, measure->Xccd, measure->Yccd);
+	if (im < 0) {
+	  // still failed, this is bad
+	  if (NmissStack < 10) fprintf (stderr, "missing stack exposure: %f %f : %d %d\n", measure->R, measure->D, measure->imageID, measure->photcode);
+	  NmissStack ++;
+	  continue;
+	}
+	// the old imageID was wrong.  replace with the new one:
+	measure->imageID = image[im].imageID;
+      } else {
+	// we have a matching image, but check that it is correct by re-projecting the pixel position
+	double Xtst, Ytst;
+	RD_to_XY (&Xtst, &Ytst, average->R, average->D, &image[im].coords);
+	
+	// find the pixel offset
+	double dX = (Xtst - measure->Xccd);
+	double dY = (Ytst - measure->Yccd);
+	
+	// if dPos is small, we have the right image
+	double dPos = hypot(dX,dY);
+	if (dPos > 10.0) {
+	  // if dPos is large, we don't the right imageID, reconstruct the imageID from
+	  // the R,D,photcode & X,Y consistency
+	  im = GetStackSeq (image, average->R, average->D, measure->photcode, measure->Xccd, measure->Yccd);
+	  if (im < 0) {
+	    // still failed, this is bad
+	    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);
@@ -48,11 +85,4 @@
     // 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
Index: /trunk/Ohana/src/relphot/src/StackImageMaps.c
===================================================================
--- /trunk/Ohana/src/relphot/src/StackImageMaps.c	(revision 39512)
+++ /trunk/Ohana/src/relphot/src/StackImageMaps.c	(revision 39513)
@@ -4,4 +4,5 @@
 
 # define D_NSTACKS 1000
+# define NGROUPS 5
 
 typedef struct {
@@ -16,7 +17,5 @@
 
 // one group for each photcode
-StackGroupType stackgroup[5];
-
-// myIndexType *warpObstimeIndex = NULL;
+StackGroupType *stackgroup = NULL;
 
 int MakeStackIndex (void) {
@@ -25,5 +24,7 @@
   Image *image = getimages (&Nimage, NULL);
 
-  for (int i = 0; i < 5; i++) {
+  ALLOCATE (stackgroup, StackGroupType, NGROUPS);
+
+  for (int i = 0; i < NGROUPS; i++) {
     stackgroup[i].photcode = 11000 + i*100;
 
@@ -49,7 +50,9 @@
       default: myAbort("impossible photcode");
     }
+    myAssert (Ng < NGROUPS, "oops");
 
     int N = stackgroup[Ng].Nstacks;
-    XY_to_RD (&stackgroup[Ng].Rcenter[N], &stackgroup[Ng].Dcenter[N], 0.5*image[i].NX, 0.5*image[i].NY, image[i].coords);
+    XY_to_RD (&stackgroup[Ng].Rcenter[N], &stackgroup[Ng].Dcenter[N], 0.5*image[i].NX, 0.5*image[i].NY, &image[i].coords);
+    stackgroup[Ng].Rcenter[N] = ohana_normalize_angle (stackgroup[Ng].Rcenter[N]);
 
     stackgroup[Ng].imageSeq[N] = i;
@@ -67,201 +70,50 @@
   }
 
-  
-
-
   // 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);
+  for (int i = 0; i < NGROUPS; i++) { 
+    sort_by_ra (stackgroup[i].Rcenter, stackgroup[i].Dcenter, stackgroup[i].imageID, stackgroup[i].imageSeq, stackgroup[i].Nstacks);
   }
 
-  // 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) {
+// find the stack which yields the given Rstk, Dstk for the given X,Y (and photcode)
+int GetStackSeq (Image *image, double Rstk, double Dstk, unsigned short photcode, float X, float Y) {
 
-  if (!warpObstimeIndex) return -1;
+  if (!stackgroup) return -1;
 
-  int seq = myIndexGetEntry (warpObstimeIndex, obstime);
-  myAssert (seq > -1, "ooops");
-  myAssert (warpgroup[seq].photcode == photcode, "oops");
-  
+  // we have the stack centers; find all stacks within 0.3 degrees of this point
+  Rstk = ohana_normalize_angle (Rstk);
+  double dD = 0.3;
+  double dR = dD / cos(RAD_DEG*Dstk);
+  double Rmin = Rstk - dR;
+  double Rmax = Rstk + dR;
+  double Dmin = Dstk - dD;
+  double Dmax = Dstk + dD;
+
+  int Ng = (photcode / 100) % 10;
+  int N = ohana_bisection_double (stackgroup[Ng].Rcenter, stackgroup[Ng].Nstacks, Rmin);
+
   double dPosMin = NAN;
-  int nPosMin = -1;
+  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;
+  for (; N < stackgroup[Ng].Nstacks; N++) {
+    if (stackgroup[Ng].Dcenter[N] < Dmin) continue;
+    if (stackgroup[Ng].Dcenter[N] > Dmax) continue;
+    if (stackgroup[Ng].Rcenter[N] > Rmax) break;
+    
+    int im = stackgroup[Ng].imageSeq[N];
 
-    // this is a possible image: check the coordinates:
-    int N = warpgroup[seq].warps[i];
-
-    // project Rave,Dave to image pixels
+    // project to image pixels
     double Xtst, Ytst;
-    RD_to_XY (&Xtst, &Ytst, Rave, Dave, &image[N].coords);
-
+    RD_to_XY (&Xtst, &Ytst, Rstk, Dstk, &image[im].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
+    // if dPos is small, we have the right image
     double dPos = hypot(dX,dY);
-    if (dPos < 20.0) return N; // 20.0 pixels = 5.0 arcsec
+    if (dPos < 10.0) return im; // 10 pixels == 2.0 arcsec
 
     // check for multiple possible matches??
@@ -270,32 +122,44 @@
     if (isnan(dPosMin)) {
       dPosMin = dPos;
-      nPosMin = i;
+      nPosMin = im;
     } else {
       if (dPos < dPosMin) {
 	dPosMin = dPos;
-	nPosMin = i;
+	nPosMin = im;
       }
     }
   }
   
-  // if we did not find any matched, give up
-  if (nPosMin < 0) return -1;
-
-  // return closest image
-  return warpgroup[seq].warps[nPosMin];
+  // we did not find the right image, return what we found (-1 or closest match)
+  return nPosMin;
 }
 
-void FreeWarpGroups (void) {
+void FreeStackGroups (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);
+  if (!stackgroup) return;
+
+  for (int i = 0; i < NGROUPS; i++) {
+    free (stackgroup[i].Rcenter);
+    free (stackgroup[i].Dcenter);
+    free (stackgroup[i].imageID);
+    free (stackgroup[i].imageSeq);
   }    
+  free (stackgroup);
+}
 
-  free (warpgroup);
-  myIndexFree (warpObstimeIndex);
+void sort_by_ra (double *R, double *D, int *I, int *S, int N) {
+
+# define SWAPFUNC(A,B){ double dtmp; int itmp; 	\
+    dtmp = R[A]; R[A] = R[B]; R[B] = dtmp;	\
+    dtmp = D[A]; D[A] = D[B]; D[B] = dtmp;	\
+    itmp = I[A]; I[A] = I[B]; I[B] = itmp;	\
+    itmp = S[A]; S[A] = S[B]; S[B] = itmp;	\
+  }
+# define COMPARE(A,B)(R[A] < R[B])
+  
+  OHANA_SORT (N, COMPARE, SWAPFUNC);
+  
+# undef SWAPFUNC
+# undef COMPARE
+  
 }
Index: /trunk/Ohana/src/relphot/src/relphot_objects.c
===================================================================
--- /trunk/Ohana/src/relphot/src/relphot_objects.c	(revision 39512)
+++ /trunk/Ohana/src/relphot/src/relphot_objects.c	(revision 39513)
@@ -40,5 +40,8 @@
   if (SYNTH_ZERO_POINTS) SynthZeroPointsLoad (SYNTH_ZERO_POINTS);
 
-  if (REPAIR_WARPS) FindWarpGroups ();
+  if (REPAIR_WARPS) { 
+    FindWarpGroups ();
+    MakeStackIndex ();
+  }
 
   // load data from each region file, only use bright stars