Index: trunk/Ohana/src/relphot/src/ImageOps.c
===================================================================
--- trunk/Ohana/src/relphot/src/ImageOps.c	(revision 27479)
+++ trunk/Ohana/src/relphot/src/ImageOps.c	(revision 27480)
@@ -1,15 +1,35 @@
 # include "relphot.h"
 
+# define USE_IMAGE_ID 1
+
+static off_t       **bin;     // link from catalog,measure to image	     
+static off_t       **clist;   // catalog which supplied measurement on image 
+static off_t       **mlist;   // measure reference for measurement on image  
+static off_t        *Nlist;   // number of measurements on image	     
+static off_t        *NLIST;   // allocated number of measurements on image   
+
+static Image        *image;   // list of available images  
+static off_t        Nimage;   // number of available images
+
+// if we search by image ID, we sort (imageIDs, imageIdx) by imageIDs to get a sorted
+// index
+
+# if USE_IMAGE_ID
+static off_t        *imageIDs; // list of all image IDs
+static off_t        *imageIdx; // list of index for image IDs 
+# else
 static unsigned int *start;
 static unsigned int *stop;
-static off_t       **bin;
-
-static off_t       **clist;
-static off_t       **mlist;
-static off_t        *Nlist;
-static off_t        *NLIST;
-
-static Image        *image;
-static off_t        Nimage;
+# endif
+
+Image *getimages (off_t *N) {
+
+  *N = Nimage;
+  return (image);
+}
+
+Image *getimage (off_t N) {
+  return (&image[N]);
+}
 
 void initImages (Image *input, off_t N) {
@@ -20,4 +40,14 @@
   Nimage = N;
 
+# if USE_IMAGE_ID
+  ALLOCATE (imageIDs, off_t, Nimage);
+  ALLOCATE (imageIdx, off_t, Nimage);
+
+  for (i = 0; i < Nimage; i++) {
+    imageIdx[i] = i;
+    imageIDs[i] = image[i].imageID;
+  }
+  llsortpair (imageIDs, imageIdx, Nimage);
+# else
   ALLOCATE (start,   unsigned, Nimage);
   ALLOCATE (stop,    unsigned, Nimage);
@@ -27,4 +57,32 @@
     stop[i]  = image[i].tzero + MAX(1.01*image[i].trate*image[i].NY, 1);
   }
+# endif
+}
+
+off_t getImageByID (off_t ID) {
+
+  // we have a pair of vectors (imageIDs, imageIdx) sorted by imageIDs
+  // use bisection to find the specified image ID
+
+# if USE_IMAGE_ID
+  off_t Nlo, Nhi, N;
+
+  Nlo = 0; Nhi = Nimage;
+  while (Nhi - Nlo > 10) {
+    N = 0.5*(Nlo + Nhi);
+    if (imageIDs[N] < ID) {
+      Nlo = MAX(N, 0);
+    } else {
+      Nhi = MIN(N + 1, Nimage);
+    }
+  }
+
+  for (N = Nlo; N < Nhi; N++) {
+    if (imageIDs[N] == ID)
+      return (imageIdx[N]);
+  }
+# endif
+
+  return (-1);
 }
 
@@ -67,84 +125,4 @@
   free (mlist);
 }
-
-/* **** use this to mimic the Mosaic version ***/
-# if (0)
-void initImageGrid (Image *image, off_t Nimage) {
-
-  /* find max dR, dD range for all mosaics */
-  /* define mosaic.coords to cover dR, dD */
-  /* send results to initGridBins */
-
-  off_t i, j, m, NX, NY;
-  int dXmax, dYmax;
-  double dS, dX, dY;
-  double R, D, Rmin, Rmax, Dmin, Dmax;
-  double Mcal, dMcal, Xm;
-
-  dXmax = dYmax = 0.0;
-  for (i = 0; i < Nmosaic; i++) {
-    Dmin = Rmin =  360.0;
-    Dmax = Rmax = -360.0;
-    dS = 0.0;
-    Mcal = dMcal = Xm = 0;
-    for (j = 0; j < Nimlist[i]; j++) {
-      m = imlist[i][j];
-      NX = image[m].NX;
-      NY = image[m].NY;
-      dS += hypot (image[m].coords.cdelt1*image[m].coords.pc1_1, image[m].coords.cdelt1*image[m].coords.pc2_1);
-      XY_to_RD (&R, &D, 0.0, 0.0, &image[m].coords);
-      Rmin = MIN (Rmin, R);
-      Rmax = MAX (Rmax, R);
-      Dmin = MIN (Dmin, D);
-      Dmax = MAX (Dmax, D);
-      XY_to_RD (&R, &D, (double) NX, 0.0, &image[m].coords);
-      Rmin = MIN (Rmin, R);
-      Rmax = MAX (Rmax, R);
-      Dmin = MIN (Dmin, D);
-      Dmax = MAX (Dmax, D);
-      XY_to_RD (&R, &D, (double) NX, (double) NY, &image[m].coords);
-      Rmin = MIN (Rmin, R);
-      Rmax = MAX (Rmax, R);
-      Dmin = MIN (Dmin, D);
-      Dmax = MAX (Dmax, D);
-      XY_to_RD (&R, &D, 0.0, (double) NY, &image[m].coords);
-      Rmin = MIN (Rmin, R);
-      Rmax = MAX (Rmax, R);
-      Dmin = MIN (Dmin, D);
-      Dmax = MAX (Dmax, D);
-      Mcal  += image[m].Mcal;
-      dMcal += image[m].dMcal;
-      Xm    += image[m].Xm;
-      /* we are using mosaic.Mcal, not image.Mcal. reset image.Mcal */
-      image[m].Mcal  = 0.0;
-      image[m].dMcal = NAN;
-      image[m].Xm    = NAN_S_SHORT;
-    }
-    dS /= Nimlist[i];
-    strcpy (mosaic[i].coords.ctype, "RA---TAN");
-    mosaic[i].coords.crval1 = Rmin;
-    mosaic[i].coords.crval2 = Dmin;
-    mosaic[i].coords.crpix1 = 0.0;
-    mosaic[i].coords.crpix2 = 0.0;
-    mosaic[i].coords.cdelt1 = dS;
-    mosaic[i].coords.cdelt2 = dS;
-    mosaic[i].coords.pc1_1  = 1.0;
-    mosaic[i].coords.pc2_2  = 1.0;
-    mosaic[i].coords.pc1_2  = 0.0;
-    mosaic[i].coords.pc2_1  = 0.0;
-    RD_to_XY (&dX, &dY, Rmax, Dmax, &mosaic[i].coords);
-
-    mosaic[i].Mcal  = Mcal / Nimlist[i];
-    mosaic[i].dMcal = dMcal / Nimlist[i];
-    mosaic[i].Xm    = Xm / Nimlist[i];
-  }
-  if (!USE_GRID) return;
-
-  dXmax = MAX (dXmax, dX);
-  dYmax = MAX (dYmax, dY);
-  initGrid (dXmax, dYmax);
-  return;
-}
-# endif
 
 /* select all image equivalent to the current photcode */
@@ -163,12 +141,93 @@
 }
 
-/* XXX convert this to use the values of measure->Xccd,Yccd */
-/* XXX need to apply the conversion for ELIXIR and LONEOS formats */
-void matchImage (Catalog *catalog, off_t meas, int cat) {
-
-  off_t i;
+int findCCD (off_t idx, off_t meas, int cat, Measure *measure) {
+
   int ccdnum;
   double X, Y;
   char *pname, *filter, *p, base[256];
+
+  /* identify the ccd on the basis of the photcode name */
+  pname = GetPhotcodeNamebyCode (image[idx].photcode);
+  filter = photcode[0].name;
+  sprintf (base, "%s.%s.", MOSAICNAME, filter);
+  if (strncmp (pname, base, strlen (base))) return (FALSE);
+  p = pname + strlen(base);
+
+  /* p is pointing at the DETECTOR ID part of the photcode: CAMERA.FILTER.DETECTOR 
+     for now, we will let this be of the form SSSSnn where SSS is an arbitrary string
+     and nn is an integer -> ccdnum */
+      
+  /// XXX we have an inconsistency here wrt dvo.layout and dvo.photcode
+
+  while (!isdigit(*p) && *p) p++;
+  if (*p == 0) return (FALSE);
+  ccdnum = atoi (p);
+
+  /* ccdnum is an integer, but not necessarily a sequence number.  when the camera layout is
+     constructed, there will be null values for undefined ccdnums */
+
+  // old code to add this measurement to the grid cell for this chip 
+  // ave = measure[0].averef;
+  // ra  = catalog[cat].average[ave].R - measure[0].dR / 3600.0;
+  // dec = catalog[cat].average[ave].D - measure[0].dD / 3600.0;
+  // RD_to_XY (&X, &Y, ra, dec, &image[i].coords);
+
+  // XXX we can now use these values (but need to be careful about old formats)
+  X = measure[0].Xccd;
+  Y = measure[0].Yccd;
+  setGridMeasure (meas, cat, X, Y, ccdnum);
+
+  return (TRUE);
+}
+
+# if USE_IMAGE_ID
+void matchImage (Catalog *catalog, off_t meas, int cat) {
+
+  off_t idx, ID;
+  int status;
+  Measure *measure;
+  
+  measure = &catalog[cat].measure[meas];
+
+  ID = measure[0].imageID;
+  idx = getImageByID (ID);
+  if (idx == -1) {
+    fprintf (stderr, "can't match detection to image?\n");
+    abort();
+  }
+
+  if (USE_GRID) {
+    status = findCCD (idx, meas, cat, measure);
+    if (!status) {
+      fprintf (stderr, "failed to determine CCD for %lld, %d\n", (long long) meas, cat);
+      abort();
+    }
+  }
+
+  // index for (catalog, measure) -> image
+  bin[cat][meas] = idx;
+
+  // index for image, Nentry -> catalog
+  clist[idx][Nlist[idx]] = cat;
+
+  // index for image, Nentry -> measure
+  mlist[idx][Nlist[idx]] = meas;
+  Nlist[idx] ++;
+
+  if (Nlist[idx] == NLIST[idx]) {
+    NLIST[idx] += 100;
+    REALLOCATE (clist[idx], off_t, NLIST[idx]);
+    REALLOCATE (mlist[idx], off_t, NLIST[idx]);
+  }	
+
+  return;
+}
+
+# else
+// this is the time-based match
+void matchImage (Catalog *catalog, off_t meas, int cat) {
+
+  off_t i;
+  int status;
   Measure *measure;
   
@@ -181,40 +240,15 @@
     
     if (USE_GRID) {
-
-      /* identify the ccd on the basis of the photcode name */
-      pname = GetPhotcodeNamebyCode (image[i].photcode);
-      filter = photcode[0].name;
-      sprintf (base, "%s.%s.", MOSAICNAME, filter);
-      if (strncmp (pname, base, strlen (base))) continue;
-      p = pname + strlen(base);
-
-      /* p is pointing at the DETECTOR ID part of the photcode: CAMERA.FILTER.DETECTOR 
-	 for now, we will let this be of the form SSSSnn where SSS is an arbitrary string
-	 and nn is an integer -> ccdnum */
-      
-      /// XXX we have an inconsistency here wrt dvo.layout and dvo.photcode
-
-      while (!isdigit(*p) && *p) p++;
-      if (*p == 0) continue;
-      ccdnum = atoi (p);
-
-      /* ccdnum is an integer, but not necessarily a sequence number.  when the camera layout is
-	 constructed, there will be null values for undefined ccdnums */
-
-      // old code to add this measurement to the grid cell for this chip 
-      // ave = measure[0].averef;
-      // ra  = catalog[cat].average[ave].R - measure[0].dR / 3600.0;
-      // dec = catalog[cat].average[ave].D - measure[0].dD / 3600.0;
-      // RD_to_XY (&X, &Y, ra, dec, &image[i].coords);
-
-      // XXX we can now use these values (but need to be careful about old formats)
-      X = measure[0].Xccd;
-      Y = measure[0].Yccd;
-      setGridMeasure (meas, cat, X, Y, ccdnum);
-    }
-
+      status = findCCD(i, meas, cat);
+      if (!status) continue;
+    }
+
+    // index for (catalog, measure) -> image
     bin[cat][meas] = i;
 
+    // index for image, Nentry -> catalog
     clist[i][Nlist[i]] = cat;
+
+    // index for image, Nentry -> measure
     mlist[i][Nlist[i]] = meas;
     Nlist[i] ++;
@@ -229,4 +263,5 @@
   /*  fprintf (stderr, "can't find source image for this measurement: %d (%d)\n", measure[0].t, measure[0].photcode); */
 }
+# endif
 
 off_t getImageEntry (off_t meas, int cat) {
@@ -585,12 +620,2 @@
 }
 
-Image *getimages (off_t *N) {
-
-  *N = Nimage;
-  return (image);
-}
-
-Image *getimage (off_t N) {
-  return (&image[N]);
-}
-
Index: trunk/Ohana/src/relphot/src/MosaicOps.c
===================================================================
--- trunk/Ohana/src/relphot/src/MosaicOps.c	(revision 27479)
+++ trunk/Ohana/src/relphot/src/MosaicOps.c	(revision 27480)
@@ -1,11 +1,14 @@
 # include "relphot.h"
 
+// array of mosaic definition structures
 static off_t  Nmosaic;
 static Mosaic *mosaic;
 
-static off_t  *Nimlist;
+// list of all images associated with a mosaic
+static off_t   *Nimlist;
 static off_t   **imlist; /* mosaic -> image[] */
 static off_t   **bin;    /* catalog, measure -> mosaic */
 
+// list of mosaic associated with each image  
 static int     **clist;  /* mosaic -> catalog[] */
 static off_t   **mlist;  /* mosiac -> measure[] */
@@ -29,4 +32,9 @@
   ALLOCATE (Nimlist, off_t,   NMOSAIC);
   ALLOCATE (NIMLIST, off_t,   NMOSAIC);
+
+  /* a 'mosaic' in relphot is (unlike relastro) a virtual concept: there is no 
+   * entry in the image table that represents this mosaic.  Instead, it is an
+   * internal construct that defines a group of related images 
+   */
 
   /* generate list of unique mosaics */
Index: trunk/Ohana/src/relphot/src/relphot.c
===================================================================
--- trunk/Ohana/src/relphot/src/relphot.c	(revision 27479)
+++ trunk/Ohana/src/relphot/src/relphot.c	(revision 27480)
@@ -1,4 +1,10 @@
 # include "relphot.h"
 # define USE_DIRECT 0
+
+# define MARKTIME(MSG,...) { \
+  float dtime; \
+  gettimeofday (&stop, (void *) NULL); \
+  dtime = DTIME (stop, start); \
+  fprintf (stderr, MSG, __VA_ARGS__); }
 
 int main (int argc, char **argv) {
@@ -7,6 +13,8 @@
   Catalog *catalog;
   FITS_DB db;
+  struct timeval start, stop;
+  SkyList *skylist = NULL;
 
-  SkyList *skylist = NULL;
+  gettimeofday (&start, (void *) NULL);
 
   /* get configuration info, args */
@@ -39,4 +47,5 @@
     if (!dvo_image_load (&db, VERBOSE, FALSE)) Shutdown ("can't read image catalog %s", db.filename);
   }
+  MARKTIME("load image data: %f sec\n", dtime);
 
   /* load regions and images based on specified sky patch */
@@ -45,4 +54,5 @@
   // XXX this is fairly lame: argv[1] is photcode if UserPatchSelect is true
   skylist = load_images (&db, argv[1], &UserPatch, UserPatchSelect);
+  MARKTIME("load images: %f sec\n", dtime);
 
   /* unlock, if we can (else, unlocked below) */
@@ -51,4 +61,5 @@
   /* load catalog data from region files */
   catalog = load_catalogs (skylist, &Ncatalog);
+  MARKTIME("load catalog data: %f sec\n", dtime);
   
   /* add in a loop over the catalogs calling dvo_catalog_chipcoords */
@@ -56,4 +67,6 @@
   /* match measurements with images, mosaics */
   initImageBins  (catalog, Ncatalog);
+  MARKTIME("make image bins: %f sec\n", dtime);
+
   initMosaicBins (catalog, Ncatalog);
   initGridBins   (catalog, Ncatalog);
@@ -61,5 +74,9 @@
 
   findImages (catalog, Ncatalog);
+  MARKTIME("set up image indexes: %f sec\n", dtime);
+
   findMosaics (catalog, Ncatalog);  /* also sets Grid values */
+  MARKTIME("set up mosaic indexes: %f sec\n", dtime);
+
   SAVEPLOT = FALSE;
 
Index: trunk/Ohana/src/relphot/src/select_images.c
===================================================================
--- trunk/Ohana/src/relphot/src/select_images.c	(revision 27479)
+++ trunk/Ohana/src/relphot/src/select_images.c	(revision 27480)
@@ -13,9 +13,12 @@
 } SkyRegionCoords;
 
+void dsortindex (double *X, off_t *Y, int N);
+off_t getRegionStartByRA (double R, double *Rref, off_t Nregions);
+
 Image *select_images (SkyList *skylist, Image *timage, off_t Ntimage, off_t **LineNumber, off_t *Nimage) {
   
   Image *image;
-  off_t i, j, k, m;
-  off_t *line_number, nimage, NIMAGE;
+  off_t i, j, k, m, nStart, iSky, nimage, NIMAGE;
+  off_t *line_number;
   int InRange, ecode, found;
   double Ri[5], Di[5], Xi[5], Yi[5], dx, dy;
@@ -23,4 +26,7 @@
   SkyRegionCoords *skycoords;
   
+  double *RmaxSky;
+  off_t *index;
+
   if (skylist[0].Nregions < 1) {
     *Nimage = 0;
@@ -38,4 +44,7 @@
 
   ALLOCATE (skycoords, SkyRegionCoords, skylist[0].Nregions);
+
+  ALLOCATE (RmaxSky, double, skylist[0].Nregions);
+  ALLOCATE (index, off_t, skylist[0].Nregions);
 
   /* compare with each region file */
@@ -56,4 +65,7 @@
     skycoords[i].Yc[4] = skycoords[i].Yc[0];    
 
+    RmaxSky[i] = skylist[0].regions[i][0].Rmax;
+    index[i] = i;
+
     dx = 0.02*(skycoords[i].Xc[2] - skycoords[i].Xc[0]);
     dy = 0.02*(skycoords[i].Yc[2] - skycoords[i].Yc[0]);
@@ -65,4 +77,6 @@
   }
 
+  dsortindex (RmaxSky, index, skylist[0].Nregions);
+
   if (VERBOSE) fprintf (stderr, "finding images\n");
   BuildChipMatch (timage, Ntimage);
@@ -86,5 +100,8 @@
     }
     
-    if (!FindMosaicForImage (timage, Ntimage, i)) continue;
+    if (!FindMosaicForImage (timage, Ntimage, i)) {
+      fprintf (stderr, "cannot find mosaic for %lld\n", (long long) i);
+      continue;
+    }
 
     /* define image corners */
@@ -96,11 +113,18 @@
     found = FALSE;
 
-    /* transform to ra,dec */
+    /* transform corners to ra,dec */
+    double RminImage = 360.0;
     for (j = 0; j < 5; j++) {
       XY_to_RD (&Ri[j], &Di[j], Xi[j], Yi[j], &timage[i].coords);
-    }
+      RminImage = MIN(RminImage, Ri[j]);
+    }
+
+    // RA(nStart) is guaranteed to be < RminImage:
+    nStart = getRegionStartByRA (RminImage, RmaxSky, skylist[0].Nregions);
 
     /* compare with each region file */
-    for (m = 0; (m < skylist[0].Nregions) && !found; m++) { 
+    for (iSky = 0; (iSky < skylist[0].Nregions) && !found; iSky++) { 
+
+      m = index[iSky];
 
       /* we make positional comparisons in the projection of catalog */
@@ -238,2 +262,35 @@
 }
 
+void dsortindex (double *X, off_t *Y, int N) {
+
+# define SWAPFUNC(A,B){ double tmpf; off_t tmpi; \
+  tmpf = X[A]; X[A] = X[B]; X[B] = tmpf; \
+  tmpi = Y[A]; Y[A] = Y[B]; Y[B] = tmpi; \
+}
+# define COMPARE(A,B)(X[A] < X[B])
+
+  OHANA_SORT (N, COMPARE, SWAPFUNC);
+
+# undef SWAPFUNC
+# undef COMPARE
+
+}
+
+off_t getRegionStartByRA (double R, double *Rref, off_t Nregions) {
+
+  // use bisection to find the overlapping mosaic
+
+  off_t Nlo, Nhi, N;
+
+  // find the last mosaic before start
+  Nlo = 0; Nhi = Nregions;
+  while (Nhi - Nlo > 10) {
+    N = 0.5*(Nlo + Nhi);
+    if (Rref[N] < R) {
+      Nlo = MAX(N, 0);
+    } else {
+      Nhi = MIN(N, Nregions);
+    }
+  }
+  return (Nlo);
+}