Index: /branches/eam_branches/ipp-20131211/Ohana/src/addstar/src/findskycell.c
===================================================================
--- /branches/eam_branches/ipp-20131211/Ohana/src/addstar/src/findskycell.c	(revision 36467)
+++ /branches/eam_branches/ipp-20131211/Ohana/src/addstar/src/findskycell.c	(revision 36468)
@@ -405,2 +405,198 @@
 }
 
+// a given tess is defined by a catdir (more than one per catdir?)
+// this function takes a catdir and generates an extension for a tess file for a local tess.
+
+int mklocal (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
+  TessallationBoundary tess;
+  
+  // user supplied values, do not try to derive from Image.dat
+  tree.NX_SUB = NX_SUB;
+  tree.NY_SUB = NY_SUB;
+  tree.dPix = SCALE/3600.0;
+
+  // derive these or assign these?
+  tree.Dmin = NAN;
+  tree.Dmax = NAN;
+  tree.Rmin = NAN;
+  tree.Rmax = NAN;
+
+  // 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],      int, tree.NBAND[zone]);
+      REALLOCATE (tree.dY[zone],      int, 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] = -1;
+	tree.dY[zone][band] = -1;
+	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]);
+  }
+
+  struct timeval start, stop;
+  gettimeofday (&start, (void *) NULL);
+
+  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);
+
+  BoundaryTreeSave (treefile, &tree);
+
+  return TRUE;
+}
+