Index: trunk/Ohana/src/relphot/Makefile
===================================================================
--- trunk/Ohana/src/relphot/Makefile	(revision 29938)
+++ trunk/Ohana/src/relphot/Makefile	(revision 30616)
@@ -29,4 +29,5 @@
 $(SRC)/StarOps.$(ARCH).o	 \
 $(SRC)/args.$(ARCH).o		 \
+$(SRC)/help.$(ARCH).o		 \
 $(SRC)/bcatalog.$(ARCH).o	 \
 $(SRC)/global_stats.$(ARCH).o	 \
Index: trunk/Ohana/src/relphot/include/relphot.h
===================================================================
--- trunk/Ohana/src/relphot/include/relphot.h	(revision 29938)
+++ trunk/Ohana/src/relphot/include/relphot.h	(revision 30616)
@@ -13,4 +13,6 @@
   float Mcal;
   float dMcal;
+  float dMsys;
+  short nFitPhotom;
   short Xm;
   float secz;
@@ -50,4 +52,5 @@
 double STAR_CHISQ;
 double MIN_ERROR;
+double IMFIT_SYS_SIGMA_LIM;
 
 int    VERBOSE;
@@ -103,4 +106,6 @@
 SkyRegion UserPatch;
 int UserPatchSelect;
+
+int USE_BASIC_CHECK;
 
 # ifdef GRID_V1
@@ -215,2 +220,5 @@
 void          write_coords        PROTO((Header *header, Coords *coords));
 int relphot_objects (void);
+
+void relphot_usage (void);
+void relphot_help (int argc, char **argv);
Index: trunk/Ohana/src/relphot/src/ConfigInit.c
===================================================================
--- trunk/Ohana/src/relphot/src/ConfigInit.c	(revision 29938)
+++ trunk/Ohana/src/relphot/src/ConfigInit.c	(revision 30616)
@@ -20,4 +20,7 @@
   GetConfig (config, "MAG_LIM",                "%lf", 0, &MAG_LIM);
   GetConfig (config, "SIGMA_LIM",              "%lf", 0, &SIGMA_LIM);
+
+  if (!ScanConfig (config, "RELPHOT_IMFIT_SYS_SIGMA_LIM", "%lf", 0, &IMFIT_SYS_SIGMA_LIM)) IMFIT_SYS_SIGMA_LIM = 0.01;
+
   GetConfig (config, "STAR_SCATTER",           "%lf", 0, &STAR_SCATTER);
   GetConfig (config, "IMAGE_SCATTER",          "%lf", 0, &IMAGE_SCATTER);
Index: trunk/Ohana/src/relphot/src/ImageOps.c
===================================================================
--- trunk/Ohana/src/relphot/src/ImageOps.c	(revision 29938)
+++ trunk/Ohana/src/relphot/src/ImageOps.c	(revision 30616)
@@ -124,4 +124,6 @@
   free (clist);
   free (mlist);
+  free (Nlist);
+  free (NLIST);
 }
 
@@ -299,6 +301,6 @@
 
   off_t i, j, m, c, n, N, Nmax, mark, bad;
-  int Nfew, Nbad, Nmos, Ngrid, Nrel, Nsys;
-  float Msys, Mrel, Mmos, Mgrid;
+  int Nfew, Nbad, Nmos, Ngrid, Nrel, Nsys, Nbright;
+  float Msys, Mrel, Mmos, Mgrid, McalBright, McalBright2;
   double *list, *dlist;
   StatType stats;
@@ -318,4 +320,8 @@
 
   Nfew = Nbad = Nmos = Ngrid = Nrel = Nsys = 0;
+
+  // counters to measure the bright-end scatter
+  McalBright = McalBright2 = 0.0;
+  Nbright = 0;
 
   for (i = 0; i < Nimage; i++) {
@@ -361,4 +367,9 @@
       list[N] = Msys - Mrel - Mmos - Mgrid;
       dlist[N] = MAX (catalog[c].measure[m].dM, MIN_ERROR);
+      if (catalog[c].measure[m].dM < IMFIT_SYS_SIGMA_LIM) {
+	McalBright += list[N];
+	McalBright2 += SQ(list[N]);
+	Nbright ++;
+      }
       N++;
     }
@@ -378,5 +389,7 @@
     liststats (list, dlist, N, &stats);
     image[i].Mcal  = stats.mean;
-    image[i].dMcal = stats.sigma;
+    image[i].dMcal = stats.error;
+    image[i].dMagSys = stats.sigma;
+    image[i].nFitPhotom = N;
     image[i].Xm    = 100.0*log10(stats.chisq);
   }
Index: trunk/Ohana/src/relphot/src/MosaicOps.c
===================================================================
--- trunk/Ohana/src/relphot/src/MosaicOps.c	(revision 29938)
+++ trunk/Ohana/src/relphot/src/MosaicOps.c	(revision 30616)
@@ -73,4 +73,5 @@
     mosaic[Nmosaic].Mcal  = 0.0;
     mosaic[Nmosaic].dMcal = 0.0;
+    mosaic[Nmosaic].dMsys = 0.0;
     mosaic[Nmosaic].Xm    = 0.0;
     mosaic[Nmosaic].flags  = image[i].flags;
@@ -189,4 +190,6 @@
       image[im].dMcal = mosaic[i].dMcal;
       image[im].Xm = mosaic[i].Xm;
+      image[im].dMagSys = mosaic[i].dMsys;
+      image[im].nFitPhotom = mosaic[i].nFitPhotom;
       image[im].flags |= (mosaic[i].flags & ID_IMAGE_PHOTOM_FEW);
       image[im].flags |= (mosaic[i].flags & ID_IMAGE_PHOTOM_POOR);
@@ -400,5 +403,5 @@
       mark = (N < IMAGE_TOOFEW) || (N < IMAGE_GOOD_FRACTION*Nlist[i]);
       if (mark) {
-	fprintf (stderr, "marked mosaic %s ("OFF_T_FMT"), ("OFF_T_FMT" < %d) || ("OFF_T_FMT" < %f*"OFF_T_FMT")\n", image[imlist[i][0]].name,  i,  N, IMAGE_TOOFEW,  N, IMAGE_GOOD_FRACTION,  Nlist[i]);
+	if (VERBOSE2) { fprintf (stderr, "marked mosaic %s ("OFF_T_FMT"), ("OFF_T_FMT" < %d) || ("OFF_T_FMT" < %f*"OFF_T_FMT")\n", image[imlist[i][0]].name,  i,  N, IMAGE_TOOFEW,  N, IMAGE_GOOD_FRACTION,  Nlist[i]); }
 	mosaic[i].flags |= ID_IMAGE_PHOTOM_FEW;
 	Nfew ++;
@@ -410,5 +413,7 @@
     if (PoorImages) fprintf (stderr, "Mmos: %f %f %d "OFF_T_FMT"\n", stats.mean, stats.sigma, stats.Nmeas,  N);
     mosaic[i].Mcal  = stats.mean;
-    mosaic[i].dMcal = stats.sigma;
+    mosaic[i].dMcal = stats.error;
+    mosaic[i].dMsys = stats.sigma;
+    mosaic[i].nFitPhotom = N;
     mosaic[i].Xm    = 100.0*log10(stats.chisq);
   }
Index: trunk/Ohana/src/relphot/src/args.c
===================================================================
--- trunk/Ohana/src/relphot/src/args.c	(revision 29938)
+++ trunk/Ohana/src/relphot/src/args.c	(revision 30616)
@@ -1,4 +1,3 @@
 # include "relphot.h"
-void usage (void);
 
 int args (int argc, char **argv) {
@@ -53,4 +52,10 @@
   }
 
+  USE_BASIC_CHECK = FALSE;
+  if ((N = get_argument (argc, argv, "-basic-image-search"))) {
+    remove_argument (N, &argc, argv);
+    USE_BASIC_CHECK = TRUE;
+  }
+
   VERBOSE = VERBOSE2 = FALSE;
   if ((N = get_argument (argc, argv, "-v"))) {
@@ -228,32 +233,6 @@
   if (UpdateAverages && (argc == 1)) return TRUE;
   if (UserPatchSelect && (argc == 2)) return TRUE;
-  if (argc != 3) usage ();
+  if (argc != 3) relphot_usage ();
 
   return TRUE;
 }
-
-void usage () {
-  fprintf (stderr, "ERROR: USAGE: relphot (region) (photcode)\n");
-  fprintf (stderr, "       or:    relphot (photcode) -region RA RA DEC DEC\n");
-  fprintf (stderr, "       or:    relphot -averages -region RA RA DEC DEC\n");
-  fprintf (stderr, "  options: \n");
-  fprintf (stderr, "  -time (start) (stop)\n");
-  fprintf (stderr, "  -v : verbose output\n");
-  fprintf (stderr, "  -vv : more verbose output\n");
-  fprintf (stderr, "  -outroot (outroot)\n");
-  fprintf (stderr, "  -plot\n");
-  fprintf (stderr, "  -plotdelay (seconds)\n");
-  fprintf (stderr, "  -statmode (mode)\n");
-  fprintf (stderr, "  -refcode (name) : give extra weight to this photcode\n");
-  fprintf (stderr, "  -n (nloop)\n");
-  fprintf (stderr, "  -reset\n");
-  fprintf (stderr, "  -update\n");
-  fprintf (stderr, "  -params\n");
-  fprintf (stderr, "  -mosaic (mosaic)\n");
-  fprintf (stderr, "  -imfreeze\n");
-  fprintf (stderr, "  -grid\n");
-  fprintf (stderr, "  -area Xmin Xmax Ymin Ymax\n");
-  fprintf (stderr, "  -instmag min max\n");
-  fprintf (stderr, "  \n");
-  exit (2);
-} 
Index: trunk/Ohana/src/relphot/src/bcatalog.c
===================================================================
--- trunk/Ohana/src/relphot/src/bcatalog.c	(revision 29938)
+++ trunk/Ohana/src/relphot/src/bcatalog.c	(revision 30616)
@@ -119,6 +119,6 @@
 
   if (VERBOSE) {
-    fprintf (stderr, "using "OFF_T_FMT" stars ("OFF_T_FMT" measures) of "OFF_T_FMT" for catalog\n", 
-	      subcatalog[0].Naverage,  subcatalog[0].Nmeasure,  i);
+    fprintf (stderr, "using "OFF_T_FMT" stars ("OFF_T_FMT" measures) of "OFF_T_FMT" for catalog %s\n", 
+	     subcatalog[0].Naverage,  subcatalog[0].Nmeasure,  i, catalog[0].filename);
     fprintf (stderr, "rejections: %d code, %d time, %d dophot, %d mag, %d sigma, %d imag, %d few\n", 
 	     Ncode, Ntime, Ndophot, Nmag, Nsigma, Nimag, Nfew);
Index: trunk/Ohana/src/relphot/src/help.c
===================================================================
--- trunk/Ohana/src/relphot/src/help.c	(revision 30616)
+++ trunk/Ohana/src/relphot/src/help.c	(revision 30616)
@@ -0,0 +1,44 @@
+# include "relphot.h"
+
+void relphot_usage (void) {
+  fprintf (stderr, "ERROR: USAGE: relphot (region) (photcode)\n");
+  fprintf (stderr, "       or:    relphot (photcode) -region RA RA DEC DEC\n");
+  fprintf (stderr, "       or:    relphot -averages -region RA RA DEC DEC\n");
+  fprintf (stderr, "  use -h for more usage information\n");
+  exit (2);
+} 
+
+void relphot_help (int argc, char **argv) {
+
+  /* check for help request */
+  if (get_argument (argc, argv, "-help")) goto show_help;
+  if (get_argument (argc, argv, "-h"))    goto show_help;
+  if (argc == 1) relphot_usage();
+  return;
+
+show_help:
+  fprintf (stderr, "ERROR: USAGE: relphot (region) (photcode)\n");
+  fprintf (stderr, "       or:    relphot (photcode) -region RA RA DEC DEC\n");
+  fprintf (stderr, "       or:    relphot -averages -region RA RA DEC DEC\n");
+  fprintf (stderr, "  options: \n");
+  fprintf (stderr, "  -time (start) (stop)\n");
+  fprintf (stderr, "  -v : verbose output\n");
+  fprintf (stderr, "  -vv : more verbose output\n");
+  fprintf (stderr, "  -outroot (outroot)\n");
+  fprintf (stderr, "  -plot\n");
+  fprintf (stderr, "  -plotdelay (seconds)\n");
+  fprintf (stderr, "  -statmode (mode)\n");
+  fprintf (stderr, "  -refcode (name) : give extra weight to this photcode\n");
+  fprintf (stderr, "  -n (nloop)\n");
+  fprintf (stderr, "  -reset\n");
+  fprintf (stderr, "  -update\n");
+  fprintf (stderr, "  -params\n");
+  fprintf (stderr, "  -mosaic (mosaic)\n");
+  fprintf (stderr, "  -imfreeze\n");
+  fprintf (stderr, "  -grid\n");
+  fprintf (stderr, "  -area Xmin Xmax Ymin Ymax\n");
+  fprintf (stderr, "  -instmag min max\n");
+  fprintf (stderr, "  \n");
+  exit (2);
+}
+
Index: trunk/Ohana/src/relphot/src/initialize.c
===================================================================
--- trunk/Ohana/src/relphot/src/initialize.c	(revision 29938)
+++ trunk/Ohana/src/relphot/src/initialize.c	(revision 30616)
@@ -5,4 +5,5 @@
   int N;
 
+  relphot_help (argc, argv);
   ConfigInit (&argc, argv);
   args (argc, argv);
Index: trunk/Ohana/src/relphot/src/load_catalogs.c
===================================================================
--- trunk/Ohana/src/relphot/src/load_catalogs.c	(revision 29938)
+++ trunk/Ohana/src/relphot/src/load_catalogs.c	(revision 30616)
@@ -6,5 +6,5 @@
   Catalog *catalog, tcatalog;
 
-  if (VERBOSE) fprintf (stderr, "loading catalog data\n");
+  if (VERBOSE2) fprintf (stderr, "loading catalog data\n");
 
   ALLOCATE (catalog, Catalog, skylist[0].Nregions);
@@ -22,9 +22,15 @@
     tcatalog.Nsecfilt  = GetPhotcodeNsecfilt ();               // set the desired number in case we need to create the catalog
 
-    if (!dvo_catalog_open (&tcatalog, skylist[0].regions[i], VERBOSE, "r")) {
+    if (!dvo_catalog_open (&tcatalog, skylist[0].regions[i], VERBOSE2, "r")) {
       fprintf (stderr, "ERROR: failure reading catalog %s\n", tcatalog.filename);
       exit (1);
     }
-    if (VERBOSE && !tcatalog.Naves_disk) fprintf (stderr, "no data in %s, skipping\n", tcatalog.filename);
+    if (!tcatalog.Naves_disk) {
+	if (VERBOSE2) { fprintf (stderr, "no data in %s, skipping\n", tcatalog.filename); }
+	dvo_catalog_unlock (&tcatalog);
+	dvo_catalog_free (&tcatalog);
+	continue;
+    }
+
     Nstar_total += tcatalog.Naverage;
     Nmeas_total += tcatalog.Nmeasure;
Index: trunk/Ohana/src/relphot/src/load_images.c
===================================================================
--- trunk/Ohana/src/relphot/src/load_images.c	(revision 29938)
+++ trunk/Ohana/src/relphot/src/load_images.c	(revision 30616)
@@ -1,3 +1,9 @@
 # include "relphot.h"
+
+# define MARKTIME(MSG,...) { \
+  float dtime; \
+  gettimeofday (&stop, (void *) NULL); \
+  dtime = DTIME (stop, start); \
+  fprintf (stderr, MSG, __VA_ARGS__); }
 
 SkyList *load_images (FITS_DB *db, char *regionName, SkyRegion *region, int RegionSelect) {
@@ -6,7 +12,10 @@
   off_t      Nimage, Nsubset, Nchar;
   off_t     *LineNumber;
+  struct timeval start, stop;
 
   SkyTable *sky = NULL;
   SkyList *skylist = NULL;
+
+  gettimeofday (&start, (void *) NULL);
 
   // load the current sky table (layout of all SkyRegions) 
@@ -29,12 +38,18 @@
       exit (2);
   }
+  MARKTIME("read image table: %f sec\n", dtime);
 
   // select the images which overlap the selected sky regions
   subset = select_images (skylist, image, Nimage, &LineNumber, &Nsubset);
+  MARKTIME("selected images: %f sec\n", dtime);
 
   gfits_vtable_from_ftable (&db[0].ftable, &db[0].vtable, LineNumber, Nsubset);
+  MARKTIME("converted ftable to vtable: %f sec\n", dtime);
 
   initImages (subset, Nsubset);
+  MARKTIME("init images: %f sec\n", dtime);
+
   initMosaics (subset, Nsubset);
+  MARKTIME("init mosaics: %f sec\n", dtime);
   
   return (skylist);
Index: trunk/Ohana/src/relphot/src/relphot.c
===================================================================
--- trunk/Ohana/src/relphot/src/relphot.c	(revision 29938)
+++ trunk/Ohana/src/relphot/src/relphot.c	(revision 30616)
@@ -85,5 +85,5 @@
   if (PLOTSTUFF) {
     plot_star_coords (catalog, Ncatalog);
-    plot_mosaic_fields (catalog);
+    // plot_mosaic_fields (catalog);
   }
 
Index: trunk/Ohana/src/relphot/src/select_images.c
===================================================================
--- trunk/Ohana/src/relphot/src/select_images.c	(revision 29938)
+++ trunk/Ohana/src/relphot/src/select_images.c	(revision 30616)
@@ -15,4 +15,10 @@
 void dsortindex (double *X, off_t *Y, int N);
 off_t getRegionStartByRA (double R, double *Rref, off_t Nregions);
+
+# define MARKTIME(MSG,...) { \
+  float dtime; \
+  gettimeofday (&stop, (void *) NULL); \
+  dtime = DTIME (stop, start); \
+  fprintf (stderr, MSG, __VA_ARGS__); }
 
 Image *select_images (SkyList *skylist, Image *timage, off_t Ntimage, off_t **LineNumber, off_t *Nimage) {
@@ -25,5 +31,8 @@
   Coords tcoords;
   SkyRegionCoords *skycoords;
-  
+  struct timeval start, stop;
+  
+  double RmaxSkyRegion, RminSkyRegion, DminSkyRegion, DmaxSkyRegion, RmidSkyRegion;
+
   double *RmaxSky;
   off_t *index;
@@ -35,4 +44,6 @@
     return NULL;
   }
+
+  gettimeofday (&start, (void *) NULL);
 
   // the comparison is made in the catalog local projection. below we set crval1,2
@@ -47,4 +58,9 @@
   ALLOCATE (RmaxSky, double, skylist[0].Nregions);
   ALLOCATE (index, off_t, skylist[0].Nregions);
+
+  RminSkyRegion = +360.0;
+  RmaxSkyRegion = -360.0;
+  DminSkyRegion = +90.0;
+  DmaxSkyRegion = -90.0;
 
   /* compare with each region file */
@@ -75,10 +91,19 @@
     skycoords[i].Xc[3] -= dx; skycoords[i].Yc[3] += dy;
     skycoords[i].Xc[4] -= dx; skycoords[i].Yc[4] -= dy;
-  }
+
+    RminSkyRegion = MIN(RminSkyRegion, skylist[0].regions[i][0].Rmin);
+    RmaxSkyRegion = MAX(RmaxSkyRegion, skylist[0].regions[i][0].Rmax);
+    DminSkyRegion = MIN(DminSkyRegion, skylist[0].regions[i][0].Dmin);
+    DmaxSkyRegion = MAX(DmaxSkyRegion, skylist[0].regions[i][0].Dmax);
+  }
+  RmidSkyRegion = 0.5*(RminSkyRegion + RmaxSkyRegion);
+  MARKTIME("create sky region coords: %f sec\n", dtime);
 
   dsortindex (RmaxSky, index, skylist[0].Nregions);
+  MARKTIME("sort sky coords: %f sec\n", dtime);
 
   if (VERBOSE) fprintf (stderr, "finding images\n");
   BuildChipMatch (timage, Ntimage);
+  MARKTIME("build chip match: %f sec\n", dtime);
 
   nimage = 0;
@@ -100,4 +125,5 @@
     }
     
+    // this adds 1.3 sec for 3M images
     if (!FindMosaicForImage (timage, Ntimage, i)) {
       fprintf (stderr, "cannot find mosaic for "OFF_T_FMT"\n", i);
@@ -113,16 +139,40 @@
     found = FALSE;
 
-    /* transform corners to ra,dec */
+    /* transform corners to ra,dec -- costs ~3sec for 3M images */
     double RminImage = 360.0;
+    double RmaxImage =   0.0;
+    double DminImage = +90.0;
+    double DmaxImage = -90.0;
     for (j = 0; j < 5; j++) {
       XY_to_RD (&Ri[j], &Di[j], Xi[j], Yi[j], &timage[i].coords);
+      Ri[j] = ohana_normalize_angle_to_midpoint (Ri[j], RmidSkyRegion);
+      
       RminImage = MIN(RminImage, Ri[j]);
-    }
-
-    // RA(nStart) is guaranteed to be < RminImage:
+      RmaxImage = MAX(RmaxImage, Ri[j]);
+      DminImage = MIN(DminImage, Di[j]);
+      DmaxImage = MAX(DmaxImage, Di[j]);
+    }
+    if (RmaxImage - RminImage > 180.0) {
+	double tmp = RminImage;
+	RmaxImage = RminImage;
+	RminImage = tmp - 360.0;
+    }
+    
+    // check that this image is even in range of the searched region
+    if (DminImage > DmaxSkyRegion) continue;
+    if (DmaxImage < DminSkyRegion) continue;
+    
+    // the sky region RA is defined to be 0 - 360.0
+    if (RminImage > RmaxSkyRegion) continue;
+    if (RmaxImage < RminSkyRegion) continue;
+
+    // image overlaps region, keep it
+    if (USE_BASIC_CHECK) goto found_it;
+
+    // RA(nStart) is guaranteed to be < RminImage: -- costs 0.5sec for 3M images
     nStart = getRegionStartByRA (RminImage, RmaxSky, skylist[0].Nregions);
 
     /* compare with each region file */
-    for (iSky = 0; (iSky < skylist[0].Nregions) && !found; iSky++) { 
+    for (iSky = nStart; (iSky < skylist[0].Nregions) && !found; iSky++) { 
 
       m = index[iSky];
@@ -142,8 +192,10 @@
       for (j = 0; (j < 4) && !found; j++) {
 	found = corner_check (&Xi[j], &Yi[j], &skycoords[m].Xc[0], &skycoords[m].Yc[0]);
+	if (found) goto found_it;
       }
       /* check if catalog corner inside image */
       for (j = 0; (j < 4) && !found; j++) {
 	found = corner_check (&skycoords[m].Xc[j], &skycoords[m].Yc[j], &Xi[0], &Yi[0]);
+	if (found) goto found_it;
       }
       /* check if edges cross */
@@ -151,28 +203,33 @@
 	for (k = 0; (k < 4) && !found; k++) {
 	  found = edge_check (&Xi[j], &Yi[j], &skycoords[m].Xc[k], &skycoords[m].Yc[k]);
+	  if (found) goto found_it;
 	}
       }
-      if (!found) continue;
-
-      image[nimage] = timage[i]; 
-      /* always allow 'few' images to succeed, if possible */
-      if (image[nimage].flags & ID_IMAGE_PHOTOM_FEW) { 
-	image[nimage].flags &= ~(ID_IMAGE_PHOTOM_FEW | ID_IMAGE_PHOTOM_POOR);
-      }
-      if (RESET) {
-	assignMcal (&image[nimage], (double *) NULL, -1);
-	image[nimage].dMcal = NAN;
-	image[nimage].flags &= ~ID_IMAGE_PHOTOM_POOR;
-      }
-      line_number[nimage] = i;
-      nimage ++;
-      if (nimage == NIMAGE) {
-	NIMAGE += 100;
-	REALLOCATE (image, Image, NIMAGE);
-	REALLOCATE (line_number, off_t, NIMAGE);
-      }
-    }
-  }
-      
+    }
+    if (!found) continue;
+
+  found_it:
+    image[nimage] = timage[i]; 
+    /* always allow 'few' images to succeed, if possible */
+    if (image[nimage].flags & ID_IMAGE_PHOTOM_FEW) { 
+      image[nimage].flags &= ~(ID_IMAGE_PHOTOM_FEW | ID_IMAGE_PHOTOM_POOR);
+    }
+    if (RESET) {
+      // XXX assignMcal (&image[nimage], (double *) NULL, -1);
+      // XXX this needs to be thought through a bit more
+      image[nimage].Mcal = 0.0;
+      image[nimage].dMcal = NAN;
+      image[nimage].flags &= ~ID_IMAGE_PHOTOM_POOR;
+    }
+    line_number[nimage] = i;
+    nimage ++;
+    if (nimage == NIMAGE) {
+      NIMAGE += 100;
+      REALLOCATE (image, Image, NIMAGE);
+      REALLOCATE (line_number, off_t, NIMAGE);
+    }
+  }
+  MARKTIME("finish image selection: %f sec\n", dtime);
+
   if (VERBOSE) fprintf (stderr, "found "OFF_T_FMT" images\n", nimage);
 
@@ -295,2 +352,21 @@
   return (Nlo);
 }
+
+off_t getRegionStopByRA (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);
+}