Changeset 36474

Timestamp:

Feb 5, 2014, 6:09:47 AM (12 years ago)

Author:

eugene

Message:

finished the tess table modifications; code seems ready to support 3pi and MD

Location:

branches/eam_branches/ipp-20131211/Ohana/src

Files:

• addstar/src/findskycell.c (modified) (9 diffs)
• libdvo/include/dvo.h (modified) (1 diff)
• libdvo/src/TessellationTable.c (modified) (14 diffs)
• relphot/include/relphot.h (modified) (2 diffs)
• relphot/src/BoundaryTreeOps.c (modified) (4 diffs)
• relphot/src/ImageOps.c (modified) (8 diffs)
• relphot/src/ImageSubset.c (modified) (9 diffs)
• relphot/src/StarOps.c (modified) (2 diffs)
• relphot/src/args.c (modified) (2 diffs)
• relphot/src/setMrelCatalog.c (modified) (2 diffs)

branches/eam_branches/ipp-20131211/Ohana/src/addstar/src/findskycell.c

@@ -41 +41 @@
 int skycellIDoff = -1;
 
+int APPEND = FALSE;
+
 int main (int argc, char **argv) {
 
@@ -55 +57 @@
   if (get_argument (argc, argv, "-h")) usage ();
 
-  NX_SUB = NY_SUB = 1;
   if ((N = get_argument (argc, argv, "-nx"))) {
     remove_argument (N, &argc, argv);
@@ -68 +69 @@
 
   /* pixel scale (arcsec/pixel) */
-  SCALE = 1.0;
   if ((N = get_argument (argc, argv, "-scale"))) {
     remove_argument (N, &argc, argv);
     SCALE = atof (argv[N]);
     remove_argument (N, &argc, argv);
+  }
+
+  /* pixel scale (arcsec/pixel) */
+  if ((N = get_argument (argc, argv, "-append"))) {
+    remove_argument (N, &argc, argv);
+    APPEND = TRUE;
   }
 
@@ -410 +416 @@
   MARKTIME("-- test %d pts: %f sec\n", Npts, dtime);
 
-  BoundaryTreeSave (treefile, &tree);
+  if (APPEND) {
+    int Ntess = 0;
+    TessellationTable *tess = TessellationTableLoad (treefile, &Ntess);
+    REALLOCATE (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;
+    }
+      
+    // add basename an related here...
+    Ntess ++;
+    TessellationTableSave (treefile, tess, Ntess);
+  } else {
+    BoundaryTreeSave (treefile, &tree);
+  }
 
   return TRUE;
@@ -416 +446 @@
 
 int apply_tree (char *treefile, char *datafile) {
+
+  int Ntess = 0;
+  TessellationTable *tess = TessellationTableLoad (treefile, &Ntess);
+  if (!tess) {
+    fprintf (stderr, "error loading tessellation table file %s\n", treefile);
+    exit (2);
+  }
+
+  FILE *f = fopen (datafile, "r");
+  if (!f) {
+    fprintf (stderr, "error opening data file %s\n", datafile);
+    exit (3);
+  }
+
+  double ra, dec;
+  int Nvalue = 0;
+  while ((Nvalue = fscanf (f, "%lf %lf", &ra, &dec)) != EOF) {
+
+    int tessID, projID, skycellID;
+
+    if (!TessellationPrimaryCellIDs (tess, Ntess, &tessID, &projID, &skycellID, ra, dec)) {
+      fprintf (stderr, "error finding cell for %f,%f\n", ra, dec);
+      continue;
+    }
+
+    fprintf (stdout, "%10.6f %10.6f : %2d  %04d %03d : %s\n", ra, dec, tessID, projID, skycellID, tess[tessID].basename);
+  }
+
+  exit (0);
+}
+
+int apply_tree_old (char *treefile, char *datafile) {
 
   BoundaryTree *tree = BoundaryTreeLoad (treefile);
@@ -479 +541 @@
 int mklocal (char *treefile, char *catdir) {
 
-  int i, status;
+  int i, j, status;
   FITS_DB db;
   Image *image;
@@ -515 +577 @@
   
   // generate an empty BoundaryTree
-  int Ntess = Nimage;
   TessellationTable *tess = NULL;
-  ALLOCATE (tess, TessellationTable, Ntess);
-  
+
+  int NtessDisk = 0;
+  if (APPEND) {
+    tess = TessellationTableLoad (treefile, &NtessDisk);
+    REALLOCATE (tess, TessellationTable, Nimage + NtessDisk);
+  } else {
+    ALLOCATE (tess, TessellationTable, Nimage + NtessDisk);
+  }
+
+  TessellationTableInit (&tess[NtessDisk], Nimage);
+
   // find the RA,DEC of the image centers & assign to cells
-  for (i = 0; i < Ntess; i++) {
+  int Ntess = NtessDisk;
+  for (i = 0; i < Nimage; i++) {
     // user supplied values, do not try to derive from Image.dat
-    tess[i].NX_SUB = NX_SUB;
-    tess[i].NY_SUB = NY_SUB;
-    tess[i].dPix = SCALE/3600.0;
+    tess[Ntess].NX_SUB = NX_SUB;
+    tess[Ntess].NY_SUB = NY_SUB;
+    tess[Ntess].dPix = SCALE/3600.0;
 
     x = 0.5*image[i].NX;
@@ -530 +601 @@
     XY_to_RD (&ra, &dec, x, y, &image[i].coords);
 
-    tess[i].Ro  = ra;
-    tess[i].Do  = dec;
-    tess[i].Xo  = x;
-    tess[i].Yo  = y;
-    tess[i].dX  = image[i].NX / NX_SUB;
-    tess[i].dY  = image[i].NY / NY_SUB;
+    tess[Ntess].Ro  = ra;
+    tess[Ntess].Do  = dec;
+    tess[Ntess].Xo  = x;
+    tess[Ntess].Yo  = y;
+    tess[Ntess].dX  = image[i].NX / NX_SUB;
+    tess[Ntess].dY  = image[i].NY / NY_SUB;
 
     // find the minimum or maximum containing region
@@ -547 +618 @@
 
       if (REGION_OPTION == REGION_MIN) { 
-        for (i = 0; i < 4; i++) {
-          R_MIN = (R[i] < ra)  ? (isfinite(R_MIN) ? MAX(R_MIN, R[i]) : R[i]) : R_MIN;
-          R_MAX = (R[i] > ra)  ? (isfinite(R_MAX) ? MIN(R_MAX, R[i]) : R[i]) : R_MAX;
-          D_MIN = (D[i] < dec) ? (isfinite(D_MIN) ? MAX(D_MIN, D[i]) : D[i]) : D_MIN;
-          D_MAX = (D[i] > dec) ? (isfinite(D_MAX) ? MIN(D_MAX, D[i]) : D[i]) : D_MAX;
+        for (j = 0; j < 4; j++) {
+          R_MIN = (R[j] < ra)  ? (isfinite(R_MIN) ? MAX(R_MIN, R[j]) : R[j]) : R_MIN;
+          R_MAX = (R[j] > ra)  ? (isfinite(R_MAX) ? MIN(R_MAX, R[j]) : R[j]) : R_MAX;
+          D_MIN = (D[j] < dec) ? (isfinite(D_MIN) ? MAX(D_MIN, D[j]) : D[j]) : D_MIN;
+          D_MAX = (D[j] > dec) ? (isfinite(D_MAX) ? MIN(D_MAX, D[j]) : D[j]) : D_MAX;
         }
       } else {
-        for (i = 0; i < 4; i++) {
-          R_MIN = (R[i] < ra)  ? (isfinite(R_MIN) ? MIN(R_MIN, R[i]) : R[i]) : R_MIN;
-          R_MAX = (R[i] > ra)  ? (isfinite(R_MAX) ? MAX(R_MAX, R[i]) : R[i]) : R_MAX;
-          D_MIN = (D[i] < dec) ? (isfinite(D_MIN) ? MIN(D_MIN, D[i]) : D[i]) : D_MIN;
-          D_MAX = (D[i] > dec) ? (isfinite(D_MAX) ? MAX(D_MAX, D[i]) : D[i]) : D_MAX;
+        for (j = 0; j < 4; j++) {
+          R_MIN = (R[j] < ra)  ? (isfinite(R_MIN) ? MIN(R_MIN, R[j]) : R[j]) : R_MIN;
+          R_MAX = (R[j] > ra)  ? (isfinite(R_MAX) ? MAX(R_MAX, R[j]) : R[j]) : R_MAX;
+          D_MIN = (D[j] < dec) ? (isfinite(D_MIN) ? MIN(D_MIN, D[j]) : D[j]) : D_MIN;
+          D_MAX = (D[j] > dec) ? (isfinite(D_MAX) ? MAX(D_MAX, D[j]) : D[j]) : D_MAX;
         }
       }
     }
-    tess[i].Rmin = R_MIN;
-    tess[i].Rmax = R_MAX;
-    tess[i].Dmin = D_MIN;
-    tess[i].Dmax = D_MAX;
-
-    tess[i].type = TESS_LOCAL;
+    tess[Ntess].Rmin = R_MIN;
+    tess[Ntess].Rmax = R_MAX;
+    tess[Ntess].Dmin = D_MIN;
+    tess[Ntess].Dmax = D_MAX;
+
+    tess[Ntess].type = TESS_LOCAL;
 
     // XXX I don't really want to do the work of discovering the rule...
-    tess[i].Nbasename = strlen(BASENAME);
-    tess[i].basename = strcreate(BASENAME);
-
-    tess[i].projectIDoff = projectIDoff;
-    tess[i].skycellIDoff = skycellIDoff;
+    tess[Ntess].Nbasename = strlen(BASENAME);
+    tess[Ntess].basename = strcreate(BASENAME);
+
+    tess[Ntess].projectIDoff = projectIDoff;
+    tess[Ntess].skycellIDoff = skycellIDoff;
+
+    Ntess ++;
   }
 

branches/eam_branches/ipp-20131211/Ohana/src/libdvo/include/dvo.h

@@ -810 +810 @@
 TessellationTable *TessellationTableLoad(char *filename, int *Ntess);
 int TessellationTableSave(char *filename, TessellationTable *tess, int Ntess);
+int TessellationPrimaryCellIDs (TessellationTable *tess, int Ntess, int *tessID, int *projID, int *skycellID, double ra, double dec);
+void TessellationTableInit (TessellationTable *tess, int Ntess);
 
 void dvo_average_init (Average *average);

branches/eam_branches/ipp-20131211/Ohana/src/libdvo/src/TessellationTable.c

@@ -10 +10 @@
 
 # define DEBUG 0
+
+void TessellationTableInit (TessellationTable *tess, int Ntess) {
+
+  int i;
+
+  for (i = 0; i < Ntess; i++) {
+    tess[i].Rmin = NAN;
+    tess[i].Rmax = NAN;
+    tess[i].Dmin = NAN;
+    tess[i].Dmax = NAN;
+
+    tess[i].Xo = NAN;
+    tess[i].Yo = NAN;
+    tess[i].Ro = NAN;
+    tess[i].Do = NAN;
+
+    tess[i].dPix = NAN;
+
+    tess[i].dX = 0;
+    tess[i].dY = 0;
+
+    tess[i].NX_SUB = 0;
+    tess[i].NY_SUB = 0;
+
+    tess[i].basename = NULL;
+    tess[i].Nbasename = 0;
+
+    tess[i].projectIDoff = 0;
+    tess[i].skycellIDoff = 0;
+
+    tess[i].type = TESS_NONE;
+    tess[i].tree = NULL;
+  }
+  return;
+}
 
 // backwards compatible load : if the first table is a BoundaryTree, not a TessellationTable, just load that.
@@ -64 +99 @@
     
     ALLOCATE (tess, TessellationTable, 1);
+    TessellationTableInit (tess, 1);
     tess[0].tree = tree;
     tess[0].type = TESS_RINGS;
+    tess[0].Rmin =   0;
+    tess[0].Rmax = 360;
+    tess[0].Dmin = -90;
+    tess[0].Dmax = +90;
     *Ntess = 1;
     return tess;
@@ -86 +126 @@
   GET_COLUMN_NEW(dY,     "Y_GRID", int);
   GET_COLUMN_NEW(Rmin,   "R_MIN",  double);
-  GET_COLUMN_NEW(Rmax,   "R_MIN",  double);
+  GET_COLUMN_NEW(Rmax,   "R_MAX",  double);
   GET_COLUMN_NEW(Dmin,   "D_MIN",  double);
   GET_COLUMN_NEW(Dmax,   "D_MAX",  double);
@@ -92 +132 @@
   GET_COLUMN_NEW(NX_SUB, "NX_SUB", int);
   GET_COLUMN_NEW(NY_SUB, "NY_SUB", int);
+  GET_COLUMN_NEW(TYPE,   "TYPE",   int);
+
+  GET_COLUMN_NEW(PROJECT_ID_OFF, "PROJECT_ID_OFF", int);
+  GET_COLUMN_NEW(SKYCELL_ID_OFF, "SKYCELL_ID_OFF", int);
+
+  GET_COLUMN_NEW(BASENAME, "BASENAME", char);
 
   *Ntess = Nrow;
@@ -109 +155 @@
     tess[i].NX_SUB = NX_SUB[i];
     tess[i].NY_SUB = NY_SUB[i];
-    tess[i].type   = TESS_LOCAL;
+    tess[i].type   = TYPE[i];
+
+    tess[i].projectIDoff = PROJECT_ID_OFF[i];
+    tess[i].skycellIDoff = SKYCELL_ID_OFF[i];
+
+    ALLOCATE (tess[i].basename, char, BOUNDARY_TREE_NAME_LENGTH);
+    memcpy (tess[i].basename, &BASENAME[i*BOUNDARY_TREE_NAME_LENGTH], BOUNDARY_TREE_NAME_LENGTH);
+    tess[i].Nbasename = strlen (tess[i].basename);
+
     tess[i].tree   = NULL;
   }
@@ -126 +180 @@
   free(NX_SUB);
   free(NY_SUB);
+  free(TYPE);
+
+  free(PROJECT_ID_OFF);
+  free(SKYCELL_ID_OFF);
 
   gfits_free_header (&theader);
@@ -134 +192 @@
   /*** check for a RINGS entry ***/
 
+  // find the RINGS tessellation (if present)
+  TessellationTable *rings = NULL;
+  for (i = 0; i < *Ntess; i++) {
+    if (tess[i].type == TESS_LOCAL) continue;
+    if (tess[i].type != TESS_RINGS) continue;
+    rings = &tess[i];
+  }
+
+  if (!rings) {
+    fclose (f);
+    gfits_free_header (&header);
+    gfits_free_matrix (&matrix);
+    return (tess);
+  }
+
   // load data for this header, or exit
   if (!gfits_load_header (f, &theader))  {
+    fprintf (stderr, "warning : RINGS tessellation listed in table but not present\n");
+    gfits_free_header (&header);
+    gfits_free_matrix (&matrix);
     return tess;
   }
@@ -142 +218 @@
   gfits_scan (&theader, "EXTNAME", "%s", 1, extname);
 
-  // is this a BoundaryTree?  if so, it must be stand-alone: load just that and generate a
-  // containing TessellationTable.
-  if (!strcmp(extname, "ZONE_DATA")) {
-    BoundaryTree *tree = BoundaryTreeRead (&header, &theader, f);
+  if (strcmp(extname, "ZONE_DATA")) {
+    fprintf (stderr, "warning : RINGS tessellation listed in table but not present\n");
     gfits_free_header (&theader);
     gfits_free_header (&header);
     gfits_free_matrix (&matrix);
-    fclose (f);
-    
-    (*Ntess) ++;
-    REALLOCATE (tess, TessellationTable, *Ntess);
-    tess[(*Ntess - 1)].tree = tree;
-    tess[(*Ntess - 1)].type = TESS_RINGS;
     return tess;
   }
-  
+
+  // is this a BoundaryTree?  if so, it must be stand-alone: load just that and generate a
+  // containing TessellationTable.
+  BoundaryTree *tree = BoundaryTreeRead (&header, &theader, f);
   gfits_free_header (&theader);
   gfits_free_header (&header);
   gfits_free_matrix (&matrix);
   fclose (f);
-
+    
+  rings->tree = tree;
   return tess;
 
@@ -237 +309 @@
   gfits_define_bintable_column (&theader, "J", "NX_SUB", "skycell subdivision in x", "none", 1.0, 0.0);
   gfits_define_bintable_column (&theader, "J", "NY_SUB", "skycell subdivision in y", "none", 1.0, 0.0);
+  gfits_define_bintable_column (&theader, "J", "TYPE",   "type of tessellation", "none", 1.0, 0.0);
+
+  gfits_define_bintable_column (&theader, "J", "PROJECT_ID_OFF", "offset in image name to projection cell ID", "none", 1.0, 0.0);
+  gfits_define_bintable_column (&theader, "J", "SKYCELL_ID_OFF", "offset in image name to skycell ID", "none", 1.0, 0.0);
+
+  char fmt[16];
+  snprintf (fmt, 16, "%dA", BOUNDARY_TREE_NAME_LENGTH);
+  gfits_define_bintable_column (&theader, fmt, "BASENAME", "image base name", "none", 1.0, 0.0);
 
   // generate the output array that carries the data
   gfits_create_table (&theader, &ftable);
   
-  int Nout = (rings == NULL) ? Ntess : Ntess - 1;
+  // int Nout = (rings == NULL) ? Ntess : Ntess - 1;
+  int Nout = Ntess;
 
   // create intermediate storage arrays
@@ -258 +339 @@
   int    *NX_SUB        ; ALLOCATE (NX_SUB,    int,    Nout);
   int    *NY_SUB        ; ALLOCATE (NY_SUB,    int,    Nout);
+  int    *TYPE          ; ALLOCATE (TYPE,      int,    Nout);
+
+  int    *PROJECT_ID_OFF; ALLOCATE (PROJECT_ID_OFF, int, Nout);
+  int    *SKYCELL_ID_OFF; ALLOCATE (SKYCELL_ID_OFF, int, Nout);
+
+  char   *BASENAME      ; ALLOCATE (BASENAME,  char,   Nout*BOUNDARY_TREE_NAME_LENGTH);
+  memset (BASENAME, 0, Nout*BOUNDARY_TREE_NAME_LENGTH);
 
   // assign the storage arrays
+  // XXX I was not writing out the RINGS tess entry, but that probably is not useful (drop
+  // N vs i in future?)
   int N = 0;
   for (i = 0; i < Ntess; i++) {
-    if (tess->type != TESS_LOCAL) continue;
+    // if (tess->type != TESS_LOCAL) continue;
     Ro[N]     = tess[i].Ro;
     Do[N]     = tess[i].Do;
@@ -276 +366 @@
     NX_SUB[N] = tess[i].NX_SUB;
     NY_SUB[N] = tess[i].NY_SUB;
+    TYPE[N]   = tess[i].type;
+
+    PROJECT_ID_OFF[N] = tess[i].projectIDoff;
+    SKYCELL_ID_OFF[N] = tess[i].skycellIDoff;
+
+    if (tess[i].basename) {
+      memcpy(&BASENAME[N*BOUNDARY_TREE_NAME_LENGTH], tess[i].basename, BOUNDARY_TREE_NAME_LENGTH);
+    }
+
     N++; 
   }
@@ -293 +392 @@
   gfits_set_bintable_column (&theader, &ftable, "NX_SUB", NX_SUB, Nout);
   gfits_set_bintable_column (&theader, &ftable, "NY_SUB", NY_SUB, Nout);
+  gfits_set_bintable_column (&theader, &ftable, "TYPE",   TYPE,   Nout);
+
+  gfits_set_bintable_column (&theader, &ftable, "PROJECT_ID_OFF", PROJECT_ID_OFF, Nout);
+  gfits_set_bintable_column (&theader, &ftable, "SKYCELL_ID_OFF", SKYCELL_ID_OFF, Nout);
+  gfits_set_bintable_column (&theader, &ftable, "BASENAME", BASENAME, Nout);
 
   free (Ro    );
@@ -307 +411 @@
   free (NX_SUB);
   free (NY_SUB);
+  free (TYPE);
+
+  free (PROJECT_ID_OFF);
+  free (SKYCELL_ID_OFF);
+  free (BASENAME);
 
   gfits_fwrite_Theader (f, &theader);
@@ -321 +430 @@
   return TRUE;
 }
+
+int TessellationLocalProjection (double *x, double *y, double r, double d, TessellationTable *tess) {
+
+    double Xo = tess->Xo;
+    double Yo = tess->Yo;
+    double Ro = tess->Ro;
+    double Do = tess->Do;
+    double dPix = tess->dPix;
+
+    // this block only depends on Ro, Do
+
+    double sdp  = sin(RAD_DEG*Do);
+    double cdp  = cos(RAD_DEG*Do);
+    double salp = sin(RAD_DEG*(r - Ro));
+    double calp = cos(RAD_DEG*(r - Ro));
+    double sdel = sin(RAD_DEG*d);
+    double cdel = cos(RAD_DEG*d);
+    
+    double stht = sdel*sdp + cdel*cdp*calp;    /* sin(theta) */
+    double sphi = cdel*salp;                   /* = cos(theta)*sin(phi) */
+    double cphi = cdel*sdp*calp - sdel*cdp;    /* = cos(theta)*cos(phi) */
+
+    // defines the TAN projection (one of zenithal projections available, libdvo/src/coordops.c
+    // R = cot (theta) = cos(theta) / sin(theta)
+    double L, M;
+    if (stht == 0) {
+        double Rc = hypot(sphi, cphi);
+        L = 180.0 * sphi / Rc;
+        M = 180.0 * cphi / Rc;
+    } else {
+        L = +DEG_RAD * sphi / stht;
+        M = -DEG_RAD * cphi / stht;
+    }
+
+    // scale, rotation, parity:
+    // rotation == 0.0 (pc1_1 == pc2_2 == 1.0, pc1_2 = pc2_1 = 0.0)
+
+    // if there were rotation or parity:
+    // Ro = (coords[0].pc1_1*coords[0].pc2_2 - coords[0].pc1_2*coords[0].pc2_1);
+    // Xo = (coords[0].pc2_2*L - coords[0].pc1_2*M) / Ro;
+    // Yo = (coords[0].pc1_1*M - coords[0].pc2_1*L) / Ro;
+
+    double X = L;
+    double Y = M;
+
+    // scale is dPix
+
+    *x = Xo - X / dPix;
+    *y = Yo + Y / dPix;
+    
+    return TRUE;
+}
+
+static int Nfail = 0;
+
+// for the given ra,dec : find the valid tessellation and the containing projection / skycell IDs
+int TessellationPrimaryCellIDs (TessellationTable *tess, int Ntess, int *tessID, int *projID, int *skycellID, double ra, double dec) {
+
+  int i, zone, band;
+
+  *tessID = -1;
+  *projID = -1;
+  *skycellID = -1;
+
+  if (!tess) return FALSE;
+
+  // find the tessellation which includes this location
+  int myTess = -1;
+  for (i = 0; i < Ntess; i++) {
+    if (ra  <  tess[i].Rmin) continue;
+    if (ra  >= tess[i].Rmax) continue;
+    if (dec <  tess[i].Dmin) continue;
+    if (dec >= tess[i].Dmax) continue;
+    myTess = i;
+    break;
+  }
+  if (myTess < 0) {
+    if (Nfail < 100) {
+      fprintf (stderr, "no matching tessellation @ %f,%f\n", ra, dec);
+    }
+    Nfail ++;
+    return FALSE;
+  }
+
+  if (tess[myTess].type == TESS_NONE) {
+    if (Nfail < 100) {
+      fprintf (stderr, "tess has invalid type\n");
+      Nfail ++;
+    }
+    return FALSE;
+  }
+
+  if (tess[myTess].type == TESS_LOCAL) {
+    // I have ra, dec, and the primary projection cell.  In order to choose the primary skycell,
+    // I just need to project to ra,dec to X,Y based on the center of the cell and then get the subdivision right.
+    
+    double x = 0.0;
+    double y = 0.0;
+    TessellationLocalProjection (&x, &y, ra, dec, &tess[myTess]);
+  
+    int xi = x / tess[myTess].dX;
+    int yi = y / tess[myTess].dY;
+    int N = xi + tess[myTess].NX_SUB * yi;
+  
+    *tessID = myTess;
+    *projID = 0;
+    *skycellID = N;
+
+    // *projID = tess[myTess].projID :: NOTE : The projection cell may be non-zero like eg
+    // STS.  However, I'm not certain we are going to implement it like this.
+
+    return TRUE;
+  }
+
+  if (tess[myTess].type == TESS_RINGS) {
+    if (!BoundaryTreeCellCoords (tess[myTess].tree, &zone, &band, ra, dec)) {
+      fprintf (stderr, "mismatch!\n");
+      return FALSE;
+    }
+  
+    // I have ra, dec, and the primary projection cell.  In order to choose the primary skycell,
+    // I just need to project to ra,dec to X,Y based on the center of the cell and then get the subdivision right.
+    
+    double x = 0.0;
+    double y = 0.0;
+    BoundaryTreeProjection (&x, &y, ra, dec, tess[myTess].tree, zone, band);
+  
+    int xi = x / tess[myTess].tree->dX[zone][band];
+    int yi = y / tess[myTess].tree->dY[zone][band];
+    int N = xi + tess[myTess].tree->NX_SUB * yi;
+  
+    *tessID = myTess;
+    *projID = tess[myTess].tree->projID[zone][band];
+    *skycellID = N;
+
+    return TRUE;
+  }
+
+  fprintf (stderr, "invalid tessellation type!\n");
+  return FALSE;
+}
+

branches/eam_branches/ipp-20131211/Ohana/src/relphot/include/relphot.h

@@ -116 +116 @@
   unsigned int photom_map_id;
   unsigned int flags;
+  int tessID;
   int projID;
   int skycellID;
@@ -401 +402 @@
 
 int MatchImageName (off_t meas, int cat, char *name);
-int MatchImageSkycellID (off_t meas, int cat, int myProjectionID, int mySkycellID);
-
-int load_tree (char *treefile);
-int BoundaryTreePrimaryCell (char *primaryCellName, double ra, double dec);
-int BoundaryTreePrimaryCellIDs (int *projID, int *skycellID, double ra, double dec);
+int MatchImageSkycellID (off_t meas, int cat, int myTessID, int myProjectionID, int mySkycellID);
+
+int load_tess (char *treefile);
+int get_tess_ids (int *tessID, int *projID, int *skycellID, double ra, double dec);
+int TessellationIDsByImageName (int *tessID, int *projID, int *skycellID, char *name);
+
+// int BoundaryTreePrimaryCell (char *primaryCellName, double ra, double dec);
+// int BoundaryTreePrimaryCellIDs (int *projID, int *skycellID, double ra, double dec);
 
 int print_measure_set_alt (Average *average, SecFilt *secfilt, Measure *measure);

branches/eam_branches/ipp-20131211/Ohana/src/relphot/src/BoundaryTreeOps.c

@@ -1 +1 @@
 # include "relphot.h"
 
-static int Nfail = 0;
 static int Ntess = 0;
-static TessllationBoundary *tess = NULL;
+static TessellationTable *tess = NULL;
 
-// for the given ra,dec : find the valid tessallation and the containing projection / skycell IDs
-int TessallationPrimaryCellIDs (int *tessID, int *projID, int *skycellID, double ra, double dec) {
+int TessellationIDsByImageName (int *tessID, int *projID, int *skycellID, char *name) {
 
-  int zone, band;
-
-  *tessID = -1;
-  *projID = -1;
-  *skycellID = -1;
-
-  if (!tess) return FALSE;
-
-  // find the tessellation which includes this location
-  int myTess = -1;
-  for (i = 0; i < Ntess; i++) {
-    if (ra  <  tess[i].Rmin) continue;
-    if (ra  >= tess[i].Rmax) continue;
-    if (dec <  tess[i].Dmin) continue;
-    if (dec >= tess[i].Dmax) continue;
-    myTess = i;
-    break;
-  }
-  if (myTess < 0) {
-    if (Nfail < 100) {
-      fprintf (stderr, "no matching tessallation @ %f,%f\n", ra, dec);
-    }
-    Nfail ++;
-    return FALSE;
-  }
-
-  if (tess[myTess].type == TESS_NONE) {
-    if (Nfail < 100) {
-      fprintf (stderr, "tess has invalid type\n");
-      Nfail ++;
-    }
-    return FALSE;
-  }
-
-  if (tess[myTess].type == TESS_LOCAL) {
-    // I have ra, dec, and the primary projection cell.  In order to choose the primary skycell,
-    // I just need to project to ra,dec to X,Y based on the center of the cell and then get the subdivision right.
-    
-    double x = 0.0;
-    double y = 0.0;
-    TessallationLocalProjection (&x, &y, ra, dec, &tess[myTess]);
-  
-    int xi = x / tess[myTess].dX;
-    int yi = y / tess[myTess].dY;
-    int N = xi + tess[myTess].NX_SUB * yi;
-  
-    *tessID = myTess;
-    *projID = tess[myTess].projID;
-    *skycellID = N;
-
-    return TRUE;
-  }
-
-  if (tess[myTess].type == TESS_RINGS) {
-    if (!BoundaryTreeCellCoords (tess[myTess].tree, &zone, &band, ra, dec)) {
-      fprintf (stderr, "mismatch!\n");
-      return FALSE;
-    }
-  
-    // I have ra, dec, and the primary projection cell.  In order to choose the primary skycell,
-    // I just need to project to ra,dec to X,Y based on the center of the cell and then get the subdivision right.
-    
-    double x = 0.0;
-    double y = 0.0;
-    BoundaryTreeProjection (&x, &y, ra, dec, tess[myTess].tree, zone, band);
-  
-    int xi = x / tess[myTess].tree->dX[zone][band];
-    int yi = y / tess[myTess].tree->dY[zone][band];
-    int N = xi + tess[myTess].tree->NX_SUB * yi;
-  
-    *tessID = myTess;
-    *projID = tess[myTess].tree->projID[zone][band];
-    *skycellID = N;
-
-    return TRUE;
-  }
-}
-
-int TessallationIDsByImageName (int *tessID, int *projID, int *skycellID, char *name) {
+  int i;
 
   *tessID = -1;
@@ -98 +18 @@
     if (!strncmp (name, tess[i].basename, tess[i].Nbasename)) {
       *tessID    = i;
-      *projID    = atoi(&name[tess[i].projectIDoff]);
-      *skycellID = atoi(&name[tess[i].skycellIDoff]);
+      if (tess[i].projectIDoff >= 0) {
+        *projID    = atoi(&name[tess[i].projectIDoff]);
+      } else {
+        *projID    = 0;
+      }
+      if (tess[i].skycellIDoff >= 0) {
+        *skycellID = atoi(&name[tess[i].skycellIDoff]);
+      } else {
+        *skycellID = 0;
+      }
       return TRUE;
     }
@@ -108 +36 @@
 int load_tess (char *tessfile) {
 
-  tess = TessBoundaryLoad (tessfile, &Ntess);
+  tess = TessellationTableLoad (tessfile, &Ntess);
   if (!tess) {
-    fprintf (stderr, "failed to load tessallation boundary file %s\n", tessfile);
+    fprintf (stderr, "failed to load tessellation boundary file %s\n", tessfile);
     exit (2);
   }
@@ -116 +44 @@
   return TRUE;
 }
+
+int get_tess_ids (int *tessID, int *projID, int *skycellID, double ra, double dec) {
+
+  int status;
+
+  if (!tess) return FALSE;
+
+  status = TessellationPrimaryCellIDs(tess, Ntess, tessID, projID, skycellID, ra, dec);
+  return status;
+}

branches/eam_branches/ipp-20131211/Ohana/src/relphot/src/ImageOps.c

@@ -48 +48 @@
 // stack image.  for now, we generate these ID arrays based on the image names when we load in the image table (initImages).
 // When we pass data to the remote clients via the ImageSubset, the projID/skycellID values are carried directly in the table.
+int *tessID    = NULL;
 int *projectID = NULL;
 int *skycellID = NULL;
@@ -81 +82 @@
 
   // for stack images, assign projection cell ID and skycell ID based on filenames
+  ALLOCATE (tessID,      int, Nimage);
   ALLOCATE (projectID,   int, Nimage);
   ALLOCATE (skycellID,   int, Nimage);
@@ -88 +90 @@
     imageIDs[i] = image[i].imageID;
 
-    tessID[i] = -1;
+    tessID[i]    = -1;
     projectID[i] = -1;
     skycellID[i] = -1;
 
-    TessallationIDsByImageName (&tessID[i], &projectID[i], &skycellID[i], image[i].name);
+    TessellationIDsByImageName (&tessID[i], &projectID[i], &skycellID[i], image[i].name);
   }
 
@@ -107 +109 @@
   // create full a Image array and save the needed values
   ALLOCATE (image, Image, N);
+
+  ALLOCATE (tessID,    int, N);
   ALLOCATE (projectID, int, N);
   ALLOCATE (skycellID, int, N);
@@ -119 +123 @@
     image[i].trate         = input[i].trate        ;
     image[i].ubercalDist   = input[i].ubercalDist  ;
+    tessID[i]              = input[i].tessID       ;
     projectID[i]           = input[i].projID       ;
     skycellID[i]           = input[i].skycellID    ;
@@ -154 +159 @@
     subset[i].trate         = image[i].trate        ;
     subset[i].ubercalDist   = image[i].ubercalDist  ;
+    subset[i].tessID        = tessID[i];
     subset[i].projID        = projectID[i];
     subset[i].skycellID     = skycellID[i];
@@ -462 +468 @@
 
 // returns image.Mcal - ff(x,y)
-int MatchImageSkycellID (off_t meas, int cat, int myProjectionID, int mySkycellID) {
+int MatchImageSkycellID (off_t meas, int cat, int myTessID, int myProjectionID, int mySkycellID) {
 
   off_t i;
@@ -471 +477 @@
   if (i == -1) return FALSE;
 
+  if (tessID[i]    == -1) return FALSE;
   if (projectID[i] == -1) return FALSE;
   if (skycellID[i] == -1) return FALSE;
 
+  if (tessID[i]    != myTessID) return FALSE;
   if (projectID[i] != myProjectionID) return FALSE;
   if (skycellID[i] != mySkycellID) return FALSE;

branches/eam_branches/ipp-20131211/Ohana/src/relphot/src/ImageSubset.c

@@ -63 +63 @@
   GET_COLUMN (map,       "PHOTOM_MAP",   int);
   GET_COLUMN (flags,     "FLAGS",        int);
+  GET_COLUMN (tessID,    "TESS_ID",      int);
   GET_COLUMN (projID,    "PROJ_ID",      int);
   GET_COLUMN (skycellID, "SKYCELL_ID",   int);
@@ -78 +79 @@
     image[i].photom_map_id = map[i];
     image[i].flags         = flags[i];
+    image[i].tessID        = tessID[i];
     image[i].projID        = projID[i];
     image[i].skycellID     = skycellID[i];
@@ -91 +93 @@
   free (map);
   free (flags);
+  free (tessID);
   free (projID);
   free (skycellID);
@@ -130 +133 @@
   gfits_define_bintable_column (&theader, "J", "FLAGS", "flags", NULL, 1.0, 1.0*0x8000);
 
+  gfits_define_bintable_column (&theader, "J", "TESS_ID", "ID", NULL, 1.0, 0.0);
   gfits_define_bintable_column (&theader, "J", "PROJ_ID", "ID", NULL, 1.0, 0.0);
   gfits_define_bintable_column (&theader, "J", "SKYCELL_ID", "ID", NULL, 1.0, 0.0);
@@ -143 +147 @@
   float *Mcal, *dMcal;
   unsigned int *imageID, *map, *flags, *tzero;
-  int *projID, *skycellID;
+  int *tessID, *projID, *skycellID;
   unsigned short *trate;
   short *ucdist;
@@ -153 +157 @@
   ALLOCATE (map,       unsigned int,   Nimage);
   ALLOCATE (flags,     unsigned int,   Nimage);
+  ALLOCATE (tessID,    int,            Nimage);
   ALLOCATE (projID,    int,            Nimage);
   ALLOCATE (skycellID, int,            Nimage);
@@ -169 +174 @@
     trate[i]     = image[i].trate;
     ucdist[i]    = image[i].ubercalDist;
+    tessID[i]    = image[i].tessID;
     projID[i]    = image[i].projID;
     skycellID[i] = image[i].skycellID;
@@ -179 +185 @@
   gfits_set_bintable_column (&theader, &ftable, "PHOTOM_MAP",   map,       Nimage);
   gfits_set_bintable_column (&theader, &ftable, "FLAGS",        flags,     Nimage);
+  gfits_set_bintable_column (&theader, &ftable, "TESS_ID",      tessID,    Nimage);
   gfits_set_bintable_column (&theader, &ftable, "PROJ_ID",      projID,    Nimage);
   gfits_set_bintable_column (&theader, &ftable, "SKYCELL_ID",   skycellID, Nimage);
@@ -190 +197 @@
   free (map);
   free (flags);
+  free (tessID);
   free (projID);
   free (skycellID);

branches/eam_branches/ipp-20131211/Ohana/src/relphot/src/StarOps.c

@@ -328 +328 @@
 int setMrel_catalog (Catalog *catalog, int Nc, int pass, FlatCorrectionTable *flatcorr, SetMrelInfo *results, int Nsecfilt) {
 
+  fprintf (stderr, "??? should you still be using the old version of setMrel_catalog??\n");
+  exit (3);
+
   off_t j, k, m, ID;
   int N;
@@ -366 +369 @@
     if (isSetMrelFinal) {
       // set the name of the primary skycell (this is used in a strcmp to match the skycells in stack detections)
-      BoundaryTreePrimaryCell(primaryCell, catalog[Nc].average[j].R, catalog[Nc].average[j].D);
+      // XXX : this whole function is deprecated
+      // BoundaryTreePrimaryCell(primaryCell, catalog[Nc].average[j].R, catalog[Nc].average[j].D);
     }
 

branches/eam_branches/ipp-20131211/Ohana/src/relphot/src/args.c

@@ -223 +223 @@
     remove_argument (N, &argc, argv);
     BOUNDARY_TREE = strcreate(argv[N]);
-    load_tree (BOUNDARY_TREE);
+    load_tess (BOUNDARY_TREE);
     remove_argument (N, &argc, argv);
   }
@@ -426 +426 @@
     remove_argument (N, &argc, argv);
     BOUNDARY_TREE = strcreate(argv[N]);
-    load_tree (BOUNDARY_TREE);
+    load_tess (BOUNDARY_TREE);
     remove_argument (N, &argc, argv);
   }

branches/eam_branches/ipp-20131211/Ohana/src/relphot/src/setMrelCatalog.c

@@ -409 +409 @@
   int projectID = -1;
   int skycellID = -1;
-  TessallationTreePrimaryCellIDs(&tessID, &projectID, &skycellID, average[0].R, average[0].D);
+  get_tess_ids(&tessID, &projectID, &skycellID, average[0].R, average[0].D);
 
   int NstackGood = 0;
@@ -477 +477 @@
       // which stack image should we use for the mean value?
       // if we request the primary (USE_TREE_FOR_PRIMARY), then find the min distances for data from the primary cell
-      if (MatchImageSkycellID (meas, cat, projectID, skycellID)) {
+      if (MatchImageSkycellID (meas, cat, tessID, projectID, skycellID)) {
         float stackPrimaryOffset = getCenterOffset (meas, cat, &measure[k], &stackImageID);
         if (stackPrimaryOffset < stackPrimaryOffsetMin) {