Index: /branches/eam_branches/ipp-20120627/Ohana/src/addstar/Makefile
===================================================================
--- /branches/eam_branches/ipp-20120627/Ohana/src/addstar/Makefile	(revision 34204)
+++ /branches/eam_branches/ipp-20120627/Ohana/src/addstar/Makefile	(revision 34205)
@@ -17,19 +17,21 @@
 FULL_LDFLAGS  = -lkapa -ldvo -lFITS -lohana $(BASE_LDFLAGS)
 
-addstar     : $(BIN)/addstar.$(ARCH)
-addstard    : $(BIN)/addstard.$(ARCH)
-addstart    : $(BIN)/addstart.$(ARCH)
-addstarc    : $(BIN)/addstarc.$(ARCH)
-mkacc-2mass : $(BIN)/mkacc-2mass.$(ARCH)
-sedstar     : $(BIN)/sedstar.$(ARCH)
-load2mass   : $(BIN)/load2mass.$(ARCH)
-loadwise    : $(BIN)/loadwise.$(ARCH)
-loadsupercos: $(BIN)/loadsupercos.$(ARCH)
-gztest      : $(BIN)/gztest.$(ARCH)
-mkcmf       : $(BIN)/mkcmf.$(ARCH)
-
-all: addstar addstar_client sedstar load2mass skycells mkcmf loadwise loadsupercos dumpskycells
-
-INSTALL = addstar addstar_client sedstar load2mass skycells mkcmf loadwise loadsupercos dumpskycells
+addstar      : $(BIN)/addstar.$(ARCH)
+addstard     : $(BIN)/addstard.$(ARCH)
+addstart     : $(BIN)/addstart.$(ARCH)
+addstarc     : $(BIN)/addstarc.$(ARCH)
+mkacc-2mass  : $(BIN)/mkacc-2mass.$(ARCH)
+sedstar      : $(BIN)/sedstar.$(ARCH)
+load2mass    : $(BIN)/load2mass.$(ARCH)
+loadwise     : $(BIN)/loadwise.$(ARCH)
+dumpskycells : $(BIN)/dumpskycells.$(ARCH)
+findskycell  : $(BIN)/findskycell.$(ARCH)
+loadsupercos : $(BIN)/loadsupercos.$(ARCH)
+gztest       : $(BIN)/gztest.$(ARCH)
+mkcmf        : $(BIN)/mkcmf.$(ARCH)
+
+all: addstar addstar_client sedstar load2mass skycells mkcmf loadwise loadsupercos dumpskycells findskycell
+
+INSTALL = addstar addstar_client sedstar load2mass skycells mkcmf loadwise loadsupercos dumpskycells findskycell
 
 # I need to fix the client/server version of addstar now that I have dropped Stars
@@ -299,4 +301,8 @@
 $(SRC)/SetSignals.$(ARCH).o \
 
+FINDSKYCELL = \
+$(SRC)/findskycell.$(ARCH).o \
+$(SRC)/Shutdown.$(ARCH).o 
+
 $(ADDSTARC)   	  : $(INC)/addstar.h
 $(ADDSTARD)   	  : $(INC)/addstar.h
@@ -306,4 +312,5 @@
 $(SKYCELLS)    	  : $(INC)/addstar.h
 $(DUMPSKYCELLS)	  : $(INC)/addstar.h
+$(FINDSKYCELL)	  : $(INC)/addstar.h
 $(LOAD-2MASS)  	  : $(INC)/addstar.h $(INC)/2mass.h
 $(LOAD-WISE)   	  : $(INC)/addstar.h $(INC)/WISE.h
@@ -322,4 +329,5 @@
 $(BIN)/skycells.$(ARCH)       : $(SKYCELLS)
 $(BIN)/dumpskycells.$(ARCH)   : $(DUMPSKYCELLS)
+$(BIN)/findskycell.$(ARCH)    : $(FINDSKYCELL)
 $(BIN)/mkcmf.$(ARCH)          : $(MKCMF)
 
Index: /branches/eam_branches/ipp-20120627/Ohana/src/addstar/src/findskycell.c
===================================================================
--- /branches/eam_branches/ipp-20120627/Ohana/src/addstar/src/findskycell.c	(revision 34205)
+++ /branches/eam_branches/ipp-20120627/Ohana/src/addstar/src/findskycell.c	(revision 34205)
@@ -0,0 +1,225 @@
+# include "skycells.h"
+
+static double RA_offset_RINGS_V3[] = {360.000000,   40.000000,  24.000000,  17.142857,  13.333333,  10.909091,   9.230769,   8.000000,   7.200000,   
+				        6.545455,   6.000000,   5.625000,   5.294118,   5.000000,   4.736842,   4.556962,   4.390244,   
+				        4.285714,   4.186047,   4.090909,   4.044944,   4.044944,   4.000000,   4.000000,   4.044944,   
+				        4.044944,   4.090909,   4.186047,   4.285714,   4.390244,   4.556962,   4.736842,   5.000000,   
+				        5.294118,   5.625000,   6.000000,   6.545455,   7.200000,   8.000000,   9.230769,  10.909091,  
+				      13.333333,  17.142857,  24.000000,  40.000000,  360.000000};
+
+// in the general case, projection cell centers are arbitrary
+// in a more specific case, DEC[i] = DEC_origin + DEC_offset*zone
+// in an even more specific case, RA[i,zone] = RA_origin[zone] + RA_offset[zone]
+
+typedef struct {
+  int FixedGridDEC;	      // is the DEC sequence linear?
+  int FixedGridRA;	      // in the RA sequence in a zone linear?
+
+  double DEC_origin;
+  double DEC_offset;
+
+  int Nzone;
+  double *RA_origin;
+  double *RA_offset;
+
+  int *Nband;
+  int *NBAND;
+
+  int **cells;
+} BoundaryTree;
+
+enum {NONE, MAKE_TREE, USE_TREE};
+
+void usage (void) {
+  fprintf (stderr, "USAGE: findcell -mktree (tree) (catdir)\n");
+  fprintf (stderr, "USAGE: findcell -tree (tree) ra dec\n");
+  exit (2);
+}
+
+int mktree (char *treefile, char *catdir);
+int find_primary_cell_coords (BoundaryTree *tree, int *zone, int *band, double ra, double dec);
+
+int main (int argc, char **argv) {
+
+  int N;
+  char *treefile = NULL;
+
+  // what does this program do?
+
+  // 1) load the image table for a tessellation and generate the boundary tree
+
+  // 2) convert RA,DEC (or list?) to cell ID 
+
+  if (get_argument (argc, argv, "-help")) usage ();
+  if (get_argument (argc, argv, "-h")) usage ();
+
+  /* extra error messages */
+  MODE = NONE;
+  if ((N = get_argument (argc, argv, "-mktree"))) {
+    MODE = MAKE_TREE;
+    remove_argument (N, &argc, argv);
+    treefile = strcreate (argv[N]);
+    remove_argument (N, &argc, argv);
+  }
+  if ((N = get_argument (argc, argv, "-tree"))) {
+    MODE = USE_TREE;
+    remove_argument (N, &argc, argv);
+    treefile = strcreate (argv[N]);
+    remove_argument (N, &argc, argv);
+  }
+
+  // generate the boundary tree
+  if (argc < 2) usage();
+  if (!MODE) usage();
+
+  if (MODE == MAKE_TREE) {
+    mktree (treefile, argv[1]);
+    exit (0);
+  }
+
+  // mktree (treefile, argv[1]);
+  exit (0);
+}
+
+int mktree (char *treefile, char *catdir) {
+
+  int i, j, zone, band, status;
+  FITS_DB db;
+  Image *image;
+  off_t Nimage;
+  double x, y, ra, dec;
+
+  char imagefile[DVO_MAX_PATH];
+  snprintf (imagefile, DVO_MAX_PATH, "%s/Images.dat", catdir);
+
+  status = dvo_image_lock (&db, imagefile, 2.0, LCK_XCLD);
+  if (!status) Shutdown ("ERROR: failure to lock image catalog %s", db.filename);
+
+  /* load or create the image table */
+  if (db.dbstate == LCK_EMPTY) Shutdown ("can't read image catalog %s", db.filename);
+
+  if (!dvo_image_load (&db, TRUE, FALSE)) Shutdown ("can't read image catalog %s", db.filename);
+
+  // convert database table to internal structure (binary to Image)
+  // 'image' points to the same memory as db->ftable->buffer
+  image = gfits_table_get_Image (&db.ftable, &Nimage, &db.swapped);
+  if (!image) {
+      fprintf (stderr, "ERROR: failed to read images\n");
+      exit (2);
+  }
+  
+  // generate an empty BoundaryTree
+  BoundaryTree tree;
+  tree.FixedGridDEC = TRUE;
+  tree.FixedGridRA = TRUE;
+  
+  // for the moment, I'm going to hardwire the DEC bands to match RINGS.V3
+  tree.DEC_origin = -92.0;
+  tree.DEC_offset =   4.0;
+
+  tree.Nzone = 46;
+
+  ALLOCATE (tree.Nband, int, tree.Nzone);
+  ALLOCATE (tree.NBAND, int, tree.Nzone);
+
+  ALLOCATE (tree.RA_origin, double, tree.Nzone);
+  ALLOCATE (tree.RA_offset, double, tree.Nzone);
+  ALLOCATE (tree.cells,      int *, tree.Nzone);
+
+  // assign the bands for RINGS.V3
+  for (zone = 0; zone < tree.Nzone; zone++) {
+    tree.Nband[zone] =  0;
+    tree.NBAND[zone] = 10;
+    tree.RA_origin[zone] = -0.5*RA_offset_RINGS_V3[zone];
+    tree.RA_offset[zone] = RA_offset_RINGS_V3[zone];
+    ALLOCATE (tree.cells[zone], int, tree.NBAND[zone]);
+    for (band = 0; band < tree.NBAND[zone]; band++) {
+      tree.cells[zone][band] = -1;
+    }
+  }
+
+  // find the RA,DEC of the image centers & assign to cells
+  for (i = 0; i < Nimage; i++) {
+    x = 0.5*image[i].NX;
+    y = 0.5*image[i].NY;
+    XY_to_RD (&ra, &dec, x, y, &image[i].coords);
+
+    if (!find_primary_cell_coords (&tree, &zone, &band, ra, dec)) {
+      fprintf (stderr, "mismatch!\n");
+      continue;
+    }
+    // fprintf (stderr, "%d  %f %f  %f  %f %f  %d %d\n", i, x, y, tree.RA_offset[zone], ra, dec, zone, band);
+    
+    if (band >= tree.NBAND[zone]) {
+      int start = tree.NBAND[zone];
+      tree.Nband[zone] = band + 1;
+      tree.NBAND[zone] = band + 10;
+      REALLOCATE (tree.cells[zone], int, tree.NBAND[zone]);
+      for (j = start; j < tree.NBAND[zone]; j++) {
+	tree.cells[zone][j] = -1;
+      }
+    }
+    tree.cells[zone][band] = i;
+  }
+
+  // figure out the max band value for each zone?
+  for (zone = 0; zone < tree.Nzone; zone++) {
+    int found_last = FALSE;
+    int last_band = -1;
+    for (band = 0; band < tree.NBAND[zone]; band++) {
+      // all cells should be filled
+      if (tree.cells[zone][band] < 0) {
+	if (!found_last) {
+	  found_last = TRUE;
+	  last_band = band;
+	}
+      } else {
+	if (found_last) {
+	  fprintf (stderr, "error: empty cell (%d,%d) after last band (%d)\n", zone, band, last_band);
+	}
+      }
+    }
+    if (last_band == -1) {
+      last_band = tree.NBAND[zone];
+    }
+    tree.Nband[zone] = last_band;
+    // fprintf (stderr, "last_band : %d, Nband: %d\n", last_band, tree.Nband[zone]);
+  }
+
+  // test : find skycell for NN random points on the sky
+  long A = time(NULL);
+  long B = A + 10000;
+  srand48(B);
+  for (i = 0; i < 10000; i++) {
+    ra  = 360.0 * drand48();
+    dec = 180.0 * drand48() - 90.0;
+    if (!find_primary_cell_coords (&tree, &zone, &band, ra, dec)) {
+      fprintf (stderr, "failure for %f,%f\n", ra, dec);
+    }
+  }
+
+  return TRUE;
+}
+
+
+// the boundary tree...
+// given an (ra,dec) pair, find the containing projection cell.
+
+int find_primary_cell_coords (BoundaryTree *tree, int *zone, int *band, double ra, double dec) {
+
+  // first, find the containing zone
+
+  // if we know dDEC, we can get the bin instantly:
+  *zone = (dec - tree->DEC_origin) / tree->DEC_offset;
+  
+  if (*zone < 0) return FALSE;
+  if (*zone >= tree->Nzone) return FALSE;
+
+  // now select the RA bin for that zone
+  *band = (ra - tree->RA_origin[*zone]) / tree->RA_offset[*zone];
+  
+  if (*band < 0) return FALSE;
+  // if (*band >= tree->Nband[*zone]) return FALSE;
+
+  return TRUE;
+}