Index: trunk/Ohana/src/addstar/src/findskycell.c
===================================================================
--- trunk/Ohana/src/addstar/src/findskycell.c	(revision 40619)
+++ trunk/Ohana/src/addstar/src/findskycell.c	(revision 40621)
@@ -12,5 +12,5 @@
 // in an even more specific case, RA[i,zone] = RA_origin[zone] + RA_offset[zone]
 
-enum {TREE_NONE, TREE_MAKE, TREE_LOCAL, TREE_USE};
+enum {TREE_NONE, TREE_MAKE, TREE_CFIS, TREE_LOCAL, TREE_USE};
 enum {REGION_NONE, REGION_USER, REGION_MIN, REGION_MAX};
 
@@ -24,4 +24,5 @@
 
 int mktree (char *treefile, char *catdir);
+int mkcfis (char *treefile, char *catdir);
 int mklocal (char *treefile, char *catdir);
 int apply_tree (char *treefile, char *datafile);
@@ -131,4 +132,10 @@
     remove_argument (N, &argc, argv);
   }
+  if ((N = get_argument (argc, argv, "-mkcfis"))) {
+    MODE = TREE_CFIS;
+    remove_argument (N, &argc, argv);
+    treefile = strcreate (argv[N]);
+    remove_argument (N, &argc, argv);
+  }
   if ((N = get_argument (argc, argv, "-mklocal"))) {
     MODE = TREE_LOCAL;
@@ -153,4 +160,9 @@
   }
 
+  if (MODE == TREE_CFIS) {
+    mkcfis (treefile, argv[1]);
+    exit (0);
+  }
+
   if (MODE == TREE_LOCAL) {
     mklocal (treefile, argv[1]);
@@ -162,5 +174,5 @@
 }
 
-int mktree (char *treefile, char *catdir) {
+int mkcfis (char *treefile, char *catdir) {
 
   int i, j, zone, band, status;
@@ -195,12 +207,14 @@
   tree.FixedGridRA = TRUE;
   
+  float CELLSIZE = 0.5;
+
   // 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;
+  tree.DEC_origin = -90.0 - 0.5*CELLSIZE;
+  tree.DEC_offset = CELLSIZE;
+
+  tree.Nzone = (180 / CELLSIZE) + 1;
   tree.NX_SUB = NX_SUB;
   tree.NY_SUB = NY_SUB;
-  tree.dPix = SCALE/3600.0;
+  tree.dPix = SCALE/3600.0; // user-supplied (0.185768447409 for CFIS)
 
   ALLOCATE (tree.Nband, int, tree.Nzone);
@@ -238,6 +252,10 @@
     tree.Nband[zone] = 1000000;
     tree.NBAND[zone] = 10;
-    tree.RA_origin[zone] = -0.5*RA_offset_RINGS_V3[zone];
-    tree.RA_offset[zone] = RA_offset_RINGS_V3[zone];
+
+    float zoneDec = -90.0 + zone * CELLSIZE;
+    float dRA = CELLSIZE / cos(RAD_DEG*zoneDec); // dRA is a size in RA degrees == \alpha_n
+
+    tree.RA_origin[zone] = 0.0 - 0.5*dRA;
+    tree.RA_offset[zone] = dRA;
 
     // set the starting values here.  In fact, the MIN and MAX values need to be tweaked 
@@ -414,4 +432,287 @@
   MARKTIME("-- test %d pts: %f sec\n", Npts, dtime);
 
+  int Ntess = 0;
+  TessellationTable *tess = NULL;
+
+  if (APPEND) {
+    tess = TessellationTableLoad (treefile, &Ntess);
+    REALLOCATE (tess, TessellationTable, Ntess + 1);
+  } else {
+    ALLOCATE (tess, TessellationTable, Ntess + 1);
+  }    
+
+  TessellationTableInit (&tess[Ntess], 1);
+  tess[Ntess].tree = &tree;
+  tess[Ntess].type = TESS_RINGS;
+  tess[Ntess].Rmin =   0;
+  tess[Ntess].Rmax = 360;
+  tess[Ntess].Dmin = -90;
+  tess[Ntess].Dmax = +90;
+
+  if (BASENAME) {
+    tess[Ntess].Nbasename = strlen(BASENAME);
+    tess[Ntess].basename = strcreate(BASENAME);
+    tess[Ntess].projectIDoff = projectIDoff;
+    tess[Ntess].skycellIDoff = skycellIDoff;
+  }
+      
+  Ntess ++;
+  TessellationTableSave (treefile, tess, Ntess);
+
+  return TRUE;
+}
+
+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);
+
+  gfits_db_init (&db);
+  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.scaledValue, &db.nativeOrder);
+  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;
+  tree.NX_SUB = NX_SUB;
+  tree.NY_SUB = NY_SUB;
+  tree.dPix = SCALE/3600.0;
+
+  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.DEC_min,     double, tree.Nzone);
+  ALLOCATE (tree.DEC_max,     double, tree.Nzone);
+  ALLOCATE (tree.DEC_min_raw, double, tree.Nzone);
+  ALLOCATE (tree.DEC_max_raw, double, tree.Nzone);
+
+  ALLOCATE (tree.ra,   double *, tree.Nzone);
+  ALLOCATE (tree.dec,  double *, tree.Nzone);
+  ALLOCATE (tree.Xo,   double *, tree.Nzone);
+  ALLOCATE (tree.Yo,   double *, tree.Nzone);
+  ALLOCATE (tree.dX,    float *, tree.Nzone);
+  ALLOCATE (tree.dY,    float *, tree.Nzone);
+  ALLOCATE (tree.cell,    int *, tree.Nzone);
+  ALLOCATE (tree.name,  char **, tree.Nzone);
+
+  // NOTE 1: the RA_origin, RA_offsets for RINGS.V3 are defined so that the first projection cell 
+  // overlaps the RA = 0,360 boundary.  This make the split a bit of a hack.  We end up
+  // with Nbands, but the max boundary of the last band only goes to 360.0 -
+  // 0.5*RA_offset.  To get the right band number for the boundary region, if the
+  // calculation for the band number lands beyond the Nbands (ie, band >= Nbands), then 
+  // we need to loop back to the first band.
+
+  // NOTE 2: when we generate the zone & bands initially, we do not know the number of
+  // bands in the end.  we cannot use the test of band >= Nband unless we set an absurdly
+  // large default value.  Thus the value of 1000000 below.
+
+  // assign the bands for RINGS.V3
+  for (zone = 0; zone < tree.Nzone; zone++) {
+    tree.Nband[zone] = 1000000;
+    tree.NBAND[zone] = 10;
+    tree.RA_origin[zone] = -0.5*RA_offset_RINGS_V3[zone];
+    tree.RA_offset[zone] = RA_offset_RINGS_V3[zone];
+
+    // set the starting values here.  In fact, the MIN and MAX values need to be tweaked 
+    // on the basis of the real images...
+    tree.DEC_min_raw[zone] = +90.0;
+    tree.DEC_max_raw[zone] = -90.0;
+    tree.DEC_min[zone] = NAN;
+    tree.DEC_max[zone] = NAN;
+
+    ALLOCATE (tree.ra[zone],   double, tree.NBAND[zone]);
+    ALLOCATE (tree.dec[zone],  double, tree.NBAND[zone]);
+    ALLOCATE (tree.Xo[zone],   double, tree.NBAND[zone]);
+    ALLOCATE (tree.Yo[zone],   double, tree.NBAND[zone]);
+    ALLOCATE (tree.dX[zone],    float, tree.NBAND[zone]);
+    ALLOCATE (tree.dY[zone],    float, tree.NBAND[zone]);
+    ALLOCATE (tree.cell[zone], int,    tree.NBAND[zone]);
+    ALLOCATE (tree.name[zone], char *, tree.NBAND[zone]);
+    for (band = 0; band < tree.NBAND[zone]; band++) {
+      tree.ra[zone][band] = NAN;
+      tree.dec[zone][band] = NAN;
+      tree.Xo[zone][band] = NAN;
+      tree.Yo[zone][band] = NAN;
+      tree.dX[zone][band] = NAN;
+      tree.dY[zone][band] = NAN;
+      tree.cell[zone][band] = -1;
+      ALLOCATE (tree.name[zone][band], char, BOUNDARY_TREE_NAME_LENGTH);
+    }
+  }
+
+  // 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 (!BoundaryTreeCellCoords (&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 + 10;
+      REALLOCATE (tree.ra[zone],   double, tree.NBAND[zone]);
+      REALLOCATE (tree.dec[zone],  double, tree.NBAND[zone]);
+      REALLOCATE (tree.Xo[zone],   double, tree.NBAND[zone]);
+      REALLOCATE (tree.Yo[zone],   double, tree.NBAND[zone]);
+      REALLOCATE (tree.dX[zone],    float, tree.NBAND[zone]);
+      REALLOCATE (tree.dY[zone],    float, tree.NBAND[zone]);
+      REALLOCATE (tree.cell[zone],    int, tree.NBAND[zone]);
+      REALLOCATE (tree.name[zone], char *, tree.NBAND[zone]);
+      for (j = start; j < tree.NBAND[zone]; j++) {
+	tree.ra[zone][j] = NAN;
+	tree.dec[zone][j] = NAN;
+	tree.Xo[zone][band] = NAN;
+	tree.Yo[zone][band] = NAN;
+	tree.dX[zone][band] = NAN;
+	tree.dY[zone][band] = NAN;
+	tree.cell[zone][j] = -1;
+	ALLOCATE (tree.name[zone][j], char, BOUNDARY_TREE_NAME_LENGTH);
+      }
+    }
+    tree.ra[zone][band] = ra;
+    tree.dec[zone][band] = dec;
+    tree.Xo[zone][band] = x;
+    tree.Yo[zone][band] = y;
+    tree.dX[zone][band] = image[i].NX / NX_SUB;
+    tree.dY[zone][band] = image[i].NY / NY_SUB;
+    
+    tree.cell[zone][band] = i;
+
+    // what are the min and max DEC values for this zone? (test center and corners of top and bottom edge)  
+
+    if (dec > 0.0) {
+      // min DEC at bottom of the cell at the center
+      x = 0.5*image[i].NX;
+      y = 0.0*image[i].NY;
+      XY_to_RD (&ra, &dec, x, y, &image[i].coords);
+      tree.DEC_min_raw[zone] = MIN(tree.DEC_min_raw[zone], dec);
+  
+      // max DEC : find the intersection between the RA boundary and the top of the cell 
+      double ra_band_min = tree.RA_origin[zone] + tree.RA_offset[zone] * band;
+  
+      // does the parity matter?
+      x = 0.0*image[i].NX;
+      y = 1.0*image[i].NY;
+  
+      int Niter;
+      for (Niter = 0; Niter < 3; Niter ++) {
+        XY_to_RD (&ra, &dec, x, y, &image[i].coords);
+        RD_to_XY (&x, &y, ra_band_min, dec, &image[i].coords);
+        y = 1.0*image[i].NY;
+      }
+      tree.DEC_max_raw[zone] = MAX(tree.DEC_max_raw[zone], dec);
+    } else {
+      // max DEC at top of the cell at the center
+      x = 0.5*image[i].NX;
+      y = 1.0*image[i].NY;
+      XY_to_RD (&ra, &dec, x, y, &image[i].coords);
+      tree.DEC_max_raw[zone] = MAX(tree.DEC_max_raw[zone], dec);
+  
+      // max DEC : find the intersection between the RA boundary and the bottom of the cell 
+      double ra_band_min = tree.RA_origin[zone] + tree.RA_offset[zone] * band;
+  
+      // does the parity matter? (NO)
+      x = 0.0*image[i].NX;
+      y = 0.0*image[i].NY;
+  
+      int Niter;
+      for (Niter = 0; Niter < 3; Niter ++) {
+        XY_to_RD (&ra, &dec, x, y, &image[i].coords);
+        RD_to_XY (&x, &y, ra_band_min, dec, &image[i].coords);
+        y = 0.0*image[i].NY;
+      }
+      tree.DEC_min_raw[zone] = MIN(tree.DEC_min_raw[zone], dec);
+    }
+    memcpy (tree.name[zone][band], image[i].name, BOUNDARY_TREE_NAME_LENGTH);
+  }
+
+  // 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.cell[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]);
+  }
+
+  // figure out the max band value for each zone?
+  int Nm = 0;
+  tree.DEC_min[Nm] = -90.0;
+  tree.DEC_max[Nm] = 0.5*(tree.DEC_max_raw[Nm] + tree.DEC_min_raw[Nm + 1]);
+
+  int Np = tree.Nzone - 1;
+  tree.DEC_min[Np] = 0.5*(tree.DEC_min_raw[Np] + tree.DEC_max_raw[Np - 1]);
+  tree.DEC_max[Np] = +90.0;
+
+  for (zone = 1; zone < tree.Nzone - 1; zone++) {
+    tree.DEC_min[zone] = 0.5*(tree.DEC_min_raw[zone] + tree.DEC_max_raw[zone - 1]);
+    tree.DEC_max[zone] = 0.5*(tree.DEC_max_raw[zone] + tree.DEC_min_raw[zone + 1]);
+  }
+
+  INITTIME;
+
+  int Npts = 10000000;
+
+  // test : find skycell for NN random points on the sky
+  long A = time(NULL);
+  long B = A + 10000;
+  srand48(B);
+  for (i = 0; i < Npts; i++) {
+    ra  = 360.0 * drand48();
+    dec = 180.0 * drand48() - 90.0;
+    if (!BoundaryTreeCellCoords (&tree, &zone, &band, ra, dec)) {
+      fprintf (stderr, "failure for %f,%f\n", ra, dec);
+    }
+  }
+  MARKTIME("-- test %d pts: %f sec\n", Npts, dtime);
+
   if (APPEND) {
     int Ntess = 0;