Index: /branches/eam_branches/ipp-20120627/Ohana/src/addstar/src/BoundaryTreeIO.c
===================================================================
--- /branches/eam_branches/ipp-20120627/Ohana/src/addstar/src/BoundaryTreeIO.c	(revision 34209)
+++ /branches/eam_branches/ipp-20120627/Ohana/src/addstar/src/BoundaryTreeIO.c	(revision 34209)
@@ -0,0 +1,264 @@
+# include "addstar.h"
+
+# define GET_COLUMN(OUT,NAME,TYPE)					\
+  TYPE *OUT = gfits_get_bintable_column_data (&theader, &ftable, NAME, type, &Nrow, &Ncol); \
+  myAssert (!strcmp(type, #TYPE), "wrong column type");
+
+// outstanding issues: 
+//   create the ra,dec,name arrays in the original tree
+//   how do we deal with the name array in gfits?
+
+BoundaryTree *BoundaryTreeLoad(char *filename) {
+
+  int i, Ncol;
+  off_t Nrow;
+  char type[16];
+  Header header;
+  Header theader;
+  Matrix matrix;
+  FTable ftable;
+
+  header.buffer = NULL;
+  matrix.buffer = NULL;
+  ftable.buffer = NULL;
+  theader.buffer = NULL;
+  BoundaryTree *tree = NULL;
+
+  FILE *f = fopen (filename, "r");
+  if (!f) {
+    fprintf (stderr, "ERROR: cannot open image subset file %s\n", filename);
+    return NULL;
+  }
+
+  /* load in PHU segment (ignore) */
+  if (!gfits_fread_header (f, &header)) {
+    if (VERBOSE) fprintf (stderr, "can't read image subset header\n");
+    goto escape;
+  }
+  if (!gfits_fread_matrix (f, &matrix, &header)) {
+    if (VERBOSE) fprintf (stderr, "can't read image subset matrix\n");
+    goto escape;
+  }
+
+  ALLOCATE (tree, BoundaryTree, 1);
+
+  gfits_scan (&header, "DEC_ORI", "%lf", 1, &tree->DEC_origin);
+  gfits_scan (&header, "DEC_OFF", "%lf", 1, &tree->DEC_offset);
+  // gfits_scan (&header, "NZONE",   "%d",  1, &tree->Nzone); XXX not needed (redundant with table length)
+
+  ftable.header = &theader;
+
+  /*** zone information table ***/
+  { 
+    // load data for this header 
+    if (!gfits_load_header (f, &theader)) goto escape;
+
+    // read the fits table bytes
+    if (!gfits_fread_ftable_data (f, &ftable, FALSE)) goto escape;
+ 
+    // need to create and assign to flat-field correction
+    GET_COLUMN(ZONE,            "ZONE",   	int);
+    GET_COLUMN(tree->RA_origin, "RA_ORIGIN",  	double);
+    GET_COLUMN(tree->RA_offset, "RA_OFFSET",  	double);
+    GET_COLUMN(tree->Nband,     "NBAND",  	int);
+    gfits_free_header (&theader);
+    gfits_free_table  (&ftable);
+
+    fprintf (stderr, "loaded data for %lld zones\n", (long long) Nrow);
+    tree->Nzone = Nrow;
+
+    // allocate the storage arrays
+    ALLOCATE (tree->ra,  double *, tree->Nzone);
+    ALLOCATE (tree->dec,  double *, tree->Nzone);
+    ALLOCATE (tree->cell, int *, tree->Nzone);
+    ALLOCATE (tree->name, char **, tree->Nzone);
+    for (i = 0; i < tree->Nzone; i++) {
+      ALLOCATE (tree->ra[i],   double *, tree->Nband[i]);
+      ALLOCATE (tree->dec[i],  double *, tree->Nband[i]);
+      ALLOCATE (tree->cell[i], int *,    tree->Nband[i]);
+      ALLOCATE (tree->name[i], char **,  tree->Nband[i]);
+    }
+  }
+
+  /*** cell information table ***/
+  { 
+
+    // load data for this header 
+    if (!gfits_load_header (f, &theader)) goto escape;
+
+    // read the fits table bytes
+    if (!gfits_fread_ftable_data (f, &ftable, FALSE)) goto escape;
+ 
+    // need to create and assign to flat-field correction
+    GET_COLUMN(R,     "RA",   	   double);
+    GET_COLUMN(D,     "DEC",  	   double);
+    GET_COLUMN(zone,  "NMEAS",       int);
+    GET_COLUMN(band,  "MEASURE_OFF", int);
+    GET_COLUMN(index, "FLAGS",       int);
+    GET_COLUMN(name,  "CAT_ID",      char); // XXX how is this done?
+    gfits_free_header (&theader);
+    gfits_free_table  (&ftable);
+    fprintf (stderr, "loaded data for %lld cells\n", (long long) Nrow);
+
+    // assign the storage arrays
+    for (i = 0; i < Nrow; i++) {
+      nz = zone[i];
+      nb = band[i];
+      tree->ra[nz][nb] = R[i];
+      tree->dec[nz][nb] = D[i];
+      tree->cells[nz][nb] = i; // XXX ?
+      tree->name[nz][nb] = name[i];
+    }
+
+    free (R     );
+    free (D     );
+    free (zone  );
+    free (band  );
+    free (index );
+    free (name  );
+  }
+
+  gfits_free_header (&header);
+  gfits_free_matrix (&matrix);
+  fclose (f);
+
+  return tree;
+
+escape:
+  gfits_free_header (&header);
+  gfits_free_matrix (&matrix);
+  gfits_free_header (&theader);
+  gfits_free_table  (&ftable);
+  if (tree) free (tree);
+
+  fclose (f);
+  return NULL;
+}
+
+// we are passed a BoundaryTree structure, write it to a FITS table (3 ext)
+int BoundaryTreeSave(char *filename, BoundaryTree *tree) {
+
+  int i;
+  Header header;
+  Header theader;
+  Matrix matrix;
+  FTable ftable;
+
+  gfits_init_header (&header);
+  header.extend = TRUE;
+  gfits_create_header (&header);
+  gfits_create_matrix (&header, &matrix);
+
+  FILE *f = fopen (filename, "w");
+  if (!f) {
+    fprintf (stderr, "ERROR: cannot open image subset file for output %s\n", filename);
+    return FALSE;
+  }
+
+  // we need some information in the header to define the layout
+  gfits_modify (&header, "DEC_ORI", "%lf", 1, tree->DEC_origin);
+  gfits_modify (&header, "DEC_OFF", "%lf", 1, tree->DEC_offset);
+  gfits_modify (&header, "NZONE",   "%d",  1, tree->Nzone);
+
+  gfits_fwrite_header  (f, &header);
+  gfits_fwrite_matrix  (f, &matrix);
+  gfits_free_header (&header);
+  gfits_free_matrix (&matrix);
+
+  /*** zone information table ***/
+  {
+    gfits_create_table_header (&theader, "BINTABLE", "ZONE_DATA");
+
+    gfits_define_bintable_column (&theader, "J", "ZONE",      "zone sequence number", NULL, 1.0, 0.0);
+    gfits_define_bintable_column (&theader, "D", "RA_ORIGIN", "origin of ra cell sequence", "degree", 1.0, 0.0);
+    gfits_define_bintable_column (&theader, "D", "RA_OFFSET", "offset per cell of ra cell sequence", "degree/cell", 1.0, 0.0);
+    gfits_define_bintable_column (&theader, "J", "NBAND",     "number of cells in each zone", NULL, 1.0, 0.0);
+
+    // generate the output array that carries the data
+    gfits_create_table (&theader, &ftable);
+
+    // create intermediate storage arrays
+    int   *zone      ; ALLOCATE (zone     ,  int  , tree->Nzone);
+
+    // assign the storage arrays
+    for (i = 0; i < tree->Nzone; i++) {
+	zone[i] = i;
+    }
+
+    // add the columns to the output array
+    gfits_set_bintable_column (&theader, &ftable, "ZONE",   	zone,            tree->Nzone);
+    gfits_set_bintable_column (&theader, &ftable, "RA_ORIGIN", 	tree->RA_origin, tree->Nzone);
+    gfits_set_bintable_column (&theader, &ftable, "RA_OFFSET", 	tree->RA_offset, tree->Nzone);
+    gfits_set_bintable_column (&theader, &ftable, "NBAND",   	tree->Nband,     tree->Nzone);
+
+    free (zone);
+
+    gfits_fwrite_Theader (f, &theader);
+    gfits_fwrite_table  (f, &ftable);
+    gfits_free_header (&theader);
+    gfits_free_table (&ftable);
+  }
+
+  /*** cell information table ***/
+  {
+    gfits_create_table_header (&theader, "BINTABLE", "CELL_DATA");
+
+    gfits_define_bintable_column (&theader, "D", "RA",   	"ra (J2000) of cell center", "degree", 1.0, 0.0);
+    gfits_define_bintable_column (&theader, "D", "DEC",  	"dec (J2000) of cell center", "degree", 1.0, 0.0);
+    gfits_define_bintable_column (&theader, "J", "ZONE",       "zone sequence number",     NULL,    1.0, 0.0);
+    gfits_define_bintable_column (&theader, "J", "BAND",       "band sequence number",  NULL,    1.0, 0.0);
+    gfits_define_bintable_column (&theader, "J", "INDEX",      "cell index",  NULL,    1.0, 0.0);
+    gfits_define_bintable_column (&theader, "J", "NAME",       "cell name",                  NULL,    1.0, 0.0);
+
+    // generate the output array that carries the data
+    gfits_create_table (&theader, &ftable);
+
+    int Ncell = 0;
+    for (i = 0; i < tree->Nzone; i++) {
+	Ncell += tree->Nband[i];
+    }
+
+    // create intermediate storage arrays
+    // NOTE: we have to unroll the 2D arrays in tree into 1D arrays
+    double *R             ; ALLOCATE (R,     double, Ncell);
+    double *D             ; ALLOCATE (D,     double, Ncell);
+    int    *zone          ; ALLOCATE (zone,  int,    Ncell);
+    int    *band          ; ALLOCATE (band,  int,    Ncell);
+    int    *index         ; ALLOCATE (index, int,    Ncell);
+    char  **name          ; ALLOCATE (name,  char *, Ncell);
+
+    // assign the storage arrays
+    i = 0;
+    for (nz = 0; nz < tree->Nzone; nz++) {
+    for (nb = 0; nb < tree->Nband[nz]; nb++) {
+      R[i]     = tree->ra[nz][nb];
+      D[i]     = tree->dec[nz][nb];
+      zone[i]  = nz;
+      band[i]  = nb;
+      index[i] = i; // or tree->cells[nz][nb] ?
+      name[i]  = tree->name[nz][nb]; // need to memcopy this?
+      i++; 
+    }
+
+    // add the columns to the output array
+    gfits_set_bintable_column (&theader, &ftable, "RA",          R,     Ncell);
+    gfits_set_bintable_column (&theader, &ftable, "DEC",         D,     Ncell);
+    gfits_set_bintable_column (&theader, &ftable, "ZONE",        zone,  Ncell);
+    gfits_set_bintable_column (&theader, &ftable, "BAND",        band,  Ncell);
+    gfits_set_bintable_column (&theader, &ftable, "INDEX",       index, Ncell);
+    gfits_set_bintable_column (&theader, &ftable, "NAME",        name,  Ncell);
+
+    free (R     );
+    free (D     );
+    free (zone  );
+    free (band  );
+    free (index );
+    free (name  );
+
+    gfits_fwrite_Theader (f, &theader);
+    gfits_fwrite_table  (f, &ftable);
+    gfits_free_header (&theader);
+    gfits_free_table (&ftable);
+  }
+  return TRUE;
+}
Index: /branches/eam_branches/ipp-20120627/Ohana/src/addstar/src/findskycell.c
===================================================================
--- /branches/eam_branches/ipp-20120627/Ohana/src/addstar/src/findskycell.c	(revision 34208)
+++ /branches/eam_branches/ipp-20120627/Ohana/src/addstar/src/findskycell.c	(revision 34209)
@@ -211,99 +211,8 @@
   MARKTIME("-- test %d pts: %f sec\n", Npts, dtime);
 
-  // save the tree
-  Header header, theader;
-  Matrix matrix;
-  FTable ftable;
-
-  gfits_init_header (&header);
-  header.extend = TRUE;
-  gfits_create_header (&header);
-  gfits_create_matrix (&header, &matrix);
-
-  ftable.header = &theader;
-  
-  FILE *f = fopen (filename, "w");
-  if (!f) {
-    fprintf (stderr, "ERROR: cannot open image subset file for output %s\n", filename);
-    return FALSE;
-  }
-
-  gfits_fwrite_header  (f, &header);
-  gfits_fwrite_matrix  (f, &matrix);
-  gfits_free_header (&header);
-  gfits_free_matrix (&matrix);
-
-  /*** MeasureTiny ***/
-  {
-    gfits_create_table_header (&theader, "BINTABLE", "MEASURE_TINY");
-
-    gfits_define_bintable_column (&theader, "E", "RA_OFF",   "ra offset",                  "arcsec", 1.0, 0.0);
-    gfits_define_bintable_column (&theader, "E", "DEC_OFF",  "dec offset",                 "arcsec", 1.0, 0.0);
-    gfits_define_bintable_column (&theader, "E", "MAG_SYS",  "magnitude (sys)",             NULL,    1.0, 0.0);
-    gfits_define_bintable_column (&theader, "E", "MAG_CAL",  "magnitude (cal)",             NULL,    1.0, 0.0);
-    gfits_define_bintable_column (&theader, "E", "MAG_ERR",  "magnitude (err)",             NULL,    1.0, 0.0);
-    gfits_define_bintable_column (&theader, "E", "AIRMASS",  "airmass",                	    NULL,    1.0, 0.0);
-    gfits_define_bintable_column (&theader, "E", "X_CCD",    "ccd x coord",            	   "pix",    1.0, 0.0);
-    gfits_define_bintable_column (&theader, "E", "Y_CCD",    "ccd y coord",            	   "pix",    1.0, 0.0);
-    gfits_define_bintable_column (&theader, "E", "EXPTIME",  "-2.5 * log (exposure time)", "sec",    1.0, 0.0);
-    gfits_define_bintable_column (&theader, "J", "TIME",     "time of exp",                "sec",    1.0, 1.0*0x8000);
-    gfits_define_bintable_column (&theader, "J", "AVE_REF",  "pointer to average table",    NULL,    1.0, 1.0*0x8000);
-    gfits_define_bintable_column (&theader, "J", "IMAGE_ID", "image",                       NULL,    1.0, 1.0*0x8000);
-    gfits_define_bintable_column (&theader, "J", "DB_FLAGS", "flags",                       NULL,    1.0, 1.0*0x8000);
-    gfits_define_bintable_column (&theader, "J", "PHOT_FLAGS", "photflags",                 NULL,    1.0, 1.0*0x8000);
-    gfits_define_bintable_column (&theader, "J", "CAT_ID",   "catalog",                     NULL,    1.0, 1.0*0x8000);
-    gfits_define_bintable_column (&theader, "I", "PHOTCODE", "photcode",                    NULL,    1.0, 1.0*0x80);
-
-    // generate the output array that carries the data
-    gfits_create_table (&theader, &ftable);
-
-    // create intermediate storage arrays
-    float *dR        ; ALLOCATE (dR       ,  float, catalog->Nmeasure);
-    float *dD        ; ALLOCATE (dD       ,  float, catalog->Nmeasure);
-    float *M         ; ALLOCATE (M        ,  float, catalog->Nmeasure);
-    float *Mcal      ; ALLOCATE (Mcal     ,  float, catalog->Nmeasure);
-    float *dM        ; ALLOCATE (dM       ,  float, catalog->Nmeasure);
-    float *airmass   ; ALLOCATE (airmass  ,  float, catalog->Nmeasure);
-    float *Xccd      ; ALLOCATE (Xccd     ,  float, catalog->Nmeasure);
-    float *Yccd      ; ALLOCATE (Yccd     ,  float, catalog->Nmeasure);
-    float *dt        ; ALLOCATE (dt       ,  float, catalog->Nmeasure);
-    int   *t         ; ALLOCATE (t        ,  int  , catalog->Nmeasure);
-    int   *averef    ; ALLOCATE (averef   ,  int  , catalog->Nmeasure);
-    int   *imageID   ; ALLOCATE (imageID  ,  int  , catalog->Nmeasure);
-    int   *dbFlags   ; ALLOCATE (dbFlags  ,  int  , catalog->Nmeasure);
-    int   *photFlags ; ALLOCATE (photFlags,  int  , catalog->Nmeasure);
-    int   *catID     ; ALLOCATE (catID    ,  int  , catalog->Nmeasure);
-    short *photcode  ; ALLOCATE (photcode ,  short, catalog->Nmeasure);
-
-    // assign the storage arrays
-    MeasureTiny *measure = catalog->measure;
-    for (i = 0; i < catalog->Nmeasure; i++) {
-      dR[i]       = measure[i].dR       ;
-      dD[i]       = measure[i].dD       ;
-      M[i]  	  = measure[i].M        ;
-      Mcal[i]     = measure[i].Mcal     ;
-      dM[i]       = measure[i].dM       ;
-      airmass[i]  = measure[i].airmass  ;
-      Xccd[i]     = measure[i].Xccd     ;
-      Yccd[i]     = measure[i].Yccd     ;
-      dt[i]       = measure[i].dt       ;
-      t[i] 	  = measure[i].t        ;
-      averef[i]   = measure[i].averef   ;
-      catID[i]    = measure[i].catID    ;
-      imageID[i]  = measure[i].imageID  ;
-      dbFlags[i]  = measure[i].dbFlags  ;
-      photFlags[i]= measure[i].photFlags;
-      photcode[i] = measure[i].photcode ;
-    }
-
-  header.buffer = NULL;
-  matrix.buffer = NULL;
-  ftable.buffer = NULL;
-  theader.buffer = NULL;
-
+  BoundaryTreeSave (treefile, &tree);
 
   return TRUE;
 }
-
 
 // the boundary tree...