Index: /branches/eam_branches/ipp-20131211/Ohana/src/addstar/src/findskycell.c =================================================================== --- /branches/eam_branches/ipp-20131211/Ohana/src/addstar/src/findskycell.c (revision 36470) +++ /branches/eam_branches/ipp-20131211/Ohana/src/addstar/src/findskycell.c (revision 36471) @@ -12,8 +12,10 @@ // in an even more specific case, RA[i,zone] = RA_origin[zone] + RA_offset[zone] -enum {TREE_NONE, TREE_MAKE, TREE_USE}; +enum {TREE_NONE, TREE_MAKE, TREE_LOCAL, TREE_USE}; +enum {REGION_NONE, REGION_USER, REGION_MIN, REGION_MAX}; void usage (void) { fprintf (stderr, "USAGE: findcell -mktree (tree) (catdir) [-nx Nx] [-ny Ny]\n"); + fprintf (stderr, "USAGE: findcell -mklocal (file) (catdir) [-nx Nx] [-ny Ny]\n"); fprintf (stderr, "USAGE: findcell -tree (tree) (datafile)\n"); fprintf (stderr, " (datafile) should contain a list of RA,DEC pairs\n"); @@ -28,4 +30,14 @@ int NY_SUB = 1; +double R_MIN = NAN; +double R_MAX = NAN; +double D_MIN = NAN; +double D_MAX = NAN; +int REGION_OPTION = REGION_NONE; + +char *BASENAME = NULL; +int projectIDoff = -1; +int skycellIDoff = -1; + int main (int argc, char **argv) { @@ -59,4 +71,42 @@ remove_argument (N, &argc, argv); SCALE = atof (argv[N]); + remove_argument (N, &argc, argv); + } + + // user-specified region + if ((N = get_argument (argc, argv, "-region"))) { + remove_argument (N, &argc, argv); + if (N > argc - 4) { + fprintf (stderr, "USAGE: -region requires 4 arguments (Rmin Rmax Dmin Dmax)\n"); + exit (1); + } + R_MIN = atof (argv[N]); remove_argument (N, &argc, argv); + R_MAX = atof (argv[N]); remove_argument (N, &argc, argv); + D_MIN = atof (argv[N]); remove_argument (N, &argc, argv); + D_MAX = atof (argv[N]); remove_argument (N, &argc, argv); + REGION_OPTION = REGION_USER; + } + + if ((N = get_argument (argc, argv, "-region-min"))) { + remove_argument (N, &argc, argv); + REGION_OPTION = REGION_MIN; + } + + if ((N = get_argument (argc, argv, "-region-max"))) { + remove_argument (N, &argc, argv); + REGION_OPTION = REGION_MAX; + } + + if ((N = get_argument (argc, argv, "-basename"))) { + remove_argument (N, &argc, argv); + if (N > argc - 3) { + fprintf (stderr, "USAGE: -basename (name) (proj_offset) (skycell_offset)\n"); + exit (1); + } + BASENAME = strcreate (argv[N]); + remove_argument (N, &argc, argv); + projectIDoff = atoi(argv[N]); + remove_argument (N, &argc, argv); + skycellIDoff = atoi(argv[N]); remove_argument (N, &argc, argv); } @@ -70,4 +120,10 @@ remove_argument (N, &argc, argv); } + if ((N = get_argument (argc, argv, "-mklocal"))) { + MODE = TREE_LOCAL; + remove_argument (N, &argc, argv); + treefile = strcreate (argv[N]); + remove_argument (N, &argc, argv); + } if ((N = get_argument (argc, argv, "-tree"))) { MODE = TREE_USE; @@ -83,4 +139,9 @@ if (MODE == TREE_MAKE) { mktree (treefile, argv[1]); + exit (0); + } + + if (MODE == TREE_LOCAL) { + mklocal (treefile, argv[1]); exit (0); } @@ -408,4 +469,11 @@ // this function takes a catdir and generates an extension for a tess file for a local tess. +/* + + LOCAL projections are defined by single projection cells. should I be supplying the info on + the cmd line or figure out the bounds from the catdir? + + */ + int mklocal (char *treefile, char *catdir) { @@ -416,4 +484,9 @@ double x, y, ra, dec; + if (REGION_OPTION == REGION_NONE) { + fprintf (stderr, "ERROR: need to define bounding region (-region Rmin Rmax Dmin Dmax | -region-min | -region-max)\n"); + exit (2); + } + char imagefile[DVO_MAX_PATH]; snprintf (imagefile, DVO_MAX_PATH, "%s/Images.dat", catdir); @@ -431,168 +504,69 @@ image = gfits_table_get_Image (&db.ftable, &Nimage, &db.swapped); if (!image) { - fprintf (stderr, "ERROR: failed to read images\n"); - exit (2); + 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; - + TessallationBoundary *tess = NULL; + ALLOCATE (tess, TessallationBoundary, Nimage); + // find the RA,DEC of the image centers & assign to cells 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; + 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; + tess[i].ra = ra; + tess[i].dec = 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; + + // find the minimum or maximum containing region + if ((REGION_OPTION == REGION_MIN) || (REGION_OPTION == REGION_MAX)) { + // XXX short cut for now : if the projection cell bounds the equator or 0,360 boundary, this will fail: + double R[4], D[4]; + XY_to_RD (&R[0], &D[0], 0, 0, &image[i].coords); + XY_to_RD (&R[1], &D[1], image[i].NX, 0, &image[i].coords); + XY_to_RD (&R[2], &D[2], 0, image[i].NY, &image[i].coords); + XY_to_RD (&R[3], &D[3], image[i].NX, image[i].NY, &image[i].coords); + + 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; } } else { - if (found_last) { - fprintf (stderr, "error: empty cell (%d,%d) after last band (%d)\n", zone, band, last_band); + 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; } } } - 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); + 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; + + // XXX I don't really want to do the work of discovering the rule... + memcpy (tess[i].basename, BASENAME); + tess[i].Nbasename = strlen(BASENAME); + + tess[i].projectIDoff = projectIDoff; + tess[i].skycellIDoff = skycellIDoff; + } BoundaryTreeSave (treefile, &tree); Index: /branches/eam_branches/ipp-20131211/Ohana/src/libdvo/include/dvo.h =================================================================== --- /branches/eam_branches/ipp-20131211/Ohana/src/libdvo/include/dvo.h (revision 36470) +++ /branches/eam_branches/ipp-20131211/Ohana/src/libdvo/include/dvo.h (revision 36471) @@ -346,6 +346,9 @@ typedef enum { TESS_NONE, TESS_LOCAL, TESS_RINGS } TessType; -// TessallationBoundaries is a structure to describe the LOCAL skycell boundaries in terms of lines of constant (RA,DEC) -// the structure describes the boundaries of a SINGLE projection cell with Nx * Ny skycells +// TessallationTable is a structure to describe the parameters of a set of "tessellations" +// (these are not strictly tessellations but projection sets as only the non-local +// versions can cover the full sky). For LOCAL projection cells, the structure describes +// the boundaries of a SINGLE projection cell with Nx * Ny skycells and includes some +// basic parameters (not used by the fullsky, eg RINGS, tessellations) typedef struct { double Rmin; // this tessellation is valid only for RA >= Rmin @@ -365,7 +368,7 @@ int NY_SUB; - TessType type; + TessType type; // BoundaryTree *tree; -} TessallationBoundary; +} TessallationTable; // a reduced-subset structure for relphot Index: /branches/eam_branches/ipp-20131211/Ohana/src/libdvo/src/TessellationTable.c =================================================================== --- /branches/eam_branches/ipp-20131211/Ohana/src/libdvo/src/TessellationTable.c (revision 36471) +++ /branches/eam_branches/ipp-20131211/Ohana/src/libdvo/src/TessellationTable.c (revision 36471) @@ -0,0 +1,329 @@ +# include "dvo.h" + +# define GET_COLUMN_NEW(OUT,NAME,TYPE) \ + TYPE *OUT = gfits_get_bintable_column_data (&theader, &ftable, NAME, type, &Nrow, &Ncol); \ + myAssert (!strcmp(type, #TYPE), "wrong column type"); + +# define GET_COLUMN_RAW(OUT,NAME,TYPE) \ + OUT = gfits_get_bintable_column_data (&theader, &ftable, NAME, type, &Nrow, &Ncol); \ + myAssert (!strcmp(type, #TYPE), "wrong column type"); + +# define DEBUG 0 + +TessellationTable *TessellationTableLoad(char *filename) { + + int i, j, nz, nb, 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; + TessellationTable *tess = 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 (DEBUG) fprintf (stderr, "can't read image subset header\n"); + goto escape; + } + if (!gfits_fread_matrix (f, &matrix, &header)) { + if (DEBUG) fprintf (stderr, "can't read image subset matrix\n"); + goto escape; + } + + ALLOCATE (tess, TessellationTable, 1); + + 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_RAW(tess->Nband, "NBAND", int); + GET_COLUMN_RAW(tess->RA_origin, "RA_ORIGIN", double); + GET_COLUMN_RAW(tess->RA_offset, "RA_OFFSET", double); + GET_COLUMN_RAW(tess->DEC_min , "DEC_MIN", double); + GET_COLUMN_RAW(tess->DEC_max , "DEC_MAX", double); + GET_COLUMN_RAW(tess->DEC_min_raw, "DEC_MIN_RAW", double); + GET_COLUMN_RAW(tess->DEC_max_raw, "DEC_MAX_RAW", double); + gfits_free_header (&theader); + gfits_free_table (&ftable); + + fprintf (stderr, "loaded data for %lld zones\n", (long long) Nrow); + tess->Nzone = Nrow; + + // allocate the storage arrays + ALLOCATE (tess->ra, double *, tess->Nzone); + ALLOCATE (tess->dec, double *, tess->Nzone); + ALLOCATE (tess->Xo, double *, tess->Nzone); + ALLOCATE (tess->Yo, double *, tess->Nzone); + ALLOCATE (tess->dX, int *, tess->Nzone); + ALLOCATE (tess->dY, int *, tess->Nzone); + ALLOCATE (tess->cell, int *, tess->Nzone); + ALLOCATE (tess->projID, int *, tess->Nzone); + ALLOCATE (tess->name, char **, tess->Nzone); + for (i = 0; i < tess->Nzone; i++) { + ALLOCATE (tess->ra[i], double, tess->Nband[i]); + ALLOCATE (tess->dec[i], double, tess->Nband[i]); + ALLOCATE (tess->Xo[i], double, tess->Nband[i]); + ALLOCATE (tess->Yo[i], double, tess->Nband[i]); + ALLOCATE (tess->dX[i], int, tess->Nband[i]); + ALLOCATE (tess->dY[i], int, tess->Nband[i]); + ALLOCATE (tess->cell[i], int, tess->Nband[i]); + ALLOCATE (tess->projID[i], int, tess->Nband[i]); + ALLOCATE (tess->name[i], char *, tess->Nband[i]); + for (j = 0; j < tess->Nband[i]; j++) { + ALLOCATE (tess->name[i][j], char, BOUNDARY_TESS_NAME_LENGTH); + } + } + } + + /*** 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_NEW(R, "RA", double); + GET_COLUMN_NEW(D, "DEC", double); + GET_COLUMN_NEW(zone, "ZONE", int); + GET_COLUMN_NEW(band, "BAND", int); + GET_COLUMN_NEW(index, "INDEX", int); + GET_COLUMN_NEW(Xo, "X_CENT", double); + GET_COLUMN_NEW(Yo, "Y_CENT", double); + GET_COLUMN_NEW(dX, "X_GRID", int); + GET_COLUMN_NEW(dY, "Y_GRID", int); + GET_COLUMN_NEW(name, "NAME", 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]; + tess->ra[nz][nb] = R[i]; + tess->dec[nz][nb] = D[i]; + tess->Xo[nz][nb] = Xo[i]; + tess->Yo[nz][nb] = Yo[i]; + tess->dX[nz][nb] = dX[i]; + tess->dY[nz][nb] = dY[i]; + tess->cell[nz][nb] = i; // XXX ? + memcpy(tess->name[nz][nb], &name[i*BOUNDARY_TESS_NAME_LENGTH], BOUNDARY_TESS_NAME_LENGTH); + // XXX parse out the ID from the name (skycell.NNNN) + tess->projID[nz][nb] = atoi(&tess->name[nz][nb][8]); + } + + free (R ); + free (D ); + free (zone ); + free (band ); + free (Xo ); + free (Yo ); + free (dX ); + free (dY ); + free (index ); + free (name ); + } + + gfits_free_header (&header); + gfits_free_matrix (&matrix); + fclose (f); + + return tess; + +escape: + gfits_free_header (&header); + gfits_free_matrix (&matrix); + gfits_free_header (&theader); + gfits_free_table (&ftable); + if (tess) free (tess); + + fclose (f); + return NULL; +} + +// we are passed a TessellationTable structure, write it to a FITS table (3 ext) +int TessellationTableSave(char *filename, TessellationTable *tess) { + + int i, nz, nb; + 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 boundary tess 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, tess->DEC_origin); + gfits_modify (&header, "DEC_OFF", "%lf", 1, tess->DEC_offset); + + gfits_modify (&header, "NX_SUB", "%d", 1, tess->NX_SUB); + gfits_modify (&header, "NY_SUB", "%d", 1, tess->NY_SUB); + gfits_modify (&header, "PIXSCALE", "%lf", 1, tess->dPix); + + 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", "none", 1.0, 0.0); + gfits_define_bintable_column (&theader, "J", "NBAND", "number of cells in each zone", "none", 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, "D", "DEC_MIN", "min dec for zone", "degree", 1.0, 0.0); + gfits_define_bintable_column (&theader, "D", "DEC_MAX", "max dec for zone", "degree", 1.0, 0.0); + gfits_define_bintable_column (&theader, "D", "DEC_MIN_RAW", "min dec for zone", "degree", 1.0, 0.0); + gfits_define_bintable_column (&theader, "D", "DEC_MAX_RAW", "max dec for zone", "degree", 1.0, 0.0); + + // generate the output array that carries the data + gfits_create_table (&theader, &ftable); + + // create intermediate storage arrays + int *zone = NULL; ALLOCATE (zone, int, tess->Nzone); + + // assign the storage arrays + for (i = 0; i < tess->Nzone; i++) { + zone[i] = i; + } + + // add the columns to the output array + gfits_set_bintable_column (&theader, &ftable, "ZONE", zone, tess->Nzone); + gfits_set_bintable_column (&theader, &ftable, "NBAND", tess->Nband, tess->Nzone); + gfits_set_bintable_column (&theader, &ftable, "RA_ORIGIN", tess->RA_origin, tess->Nzone); + gfits_set_bintable_column (&theader, &ftable, "RA_OFFSET", tess->RA_offset, tess->Nzone); + gfits_set_bintable_column (&theader, &ftable, "DEC_MIN", tess->DEC_min, tess->Nzone); + gfits_set_bintable_column (&theader, &ftable, "DEC_MAX", tess->DEC_max, tess->Nzone); + gfits_set_bintable_column (&theader, &ftable, "DEC_MIN_RAW", tess->DEC_min_raw, tess->Nzone); + gfits_set_bintable_column (&theader, &ftable, "DEC_MAX_RAW", tess->DEC_max_raw, tess->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"); + + char fmt[16]; + snprintf (fmt, 16, "%dA", BOUNDARY_TESS_NAME_LENGTH); + 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", "none", 1.0, 0.0); + gfits_define_bintable_column (&theader, "J", "BAND", "band sequence number", "none", 1.0, 0.0); + gfits_define_bintable_column (&theader, "J", "INDEX","cell index", "none", 1.0, 0.0); + gfits_define_bintable_column (&theader, "D", "X_CENT", "projection cell center pixel", "none", 1.0, 0.0); + gfits_define_bintable_column (&theader, "D", "Y_CENT", "projection cell center pixel", "none", 1.0, 0.0); + gfits_define_bintable_column (&theader, "J", "X_GRID", "skycell grid spacing", "none", 1.0, 0.0); + gfits_define_bintable_column (&theader, "J", "Y_GRID", "skycell grid spacing", "none", 1.0, 0.0); + gfits_define_bintable_column (&theader, fmt, "NAME", "cell name", "none", 1.0, 0.0); + + // generate the output array that carries the data + gfits_create_table (&theader, &ftable); + + int Ncell = 0; + for (i = 0; i < tess->Nzone; i++) { + Ncell += tess->Nband[i]; + } + + // create intermediate storage arrays + // NOTE: we have to unroll the 2D arrays in tess 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); + double *Xo ; ALLOCATE (Xo, double, Ncell); + double *Yo ; ALLOCATE (Yo, double, Ncell); + int *dX ; ALLOCATE (dX, int, Ncell); + int *dY ; ALLOCATE (dY, int, Ncell); + char *name ; ALLOCATE (name, char, Ncell*BOUNDARY_TESS_NAME_LENGTH); + + // NOTE: a table column of characters must be fixed width, and is passed as a + // contiguous array of Nchar * Nrow values + + // assign the storage arrays + i = 0; + for (nz = 0; nz < tess->Nzone; nz++) { + for (nb = 0; nb < tess->Nband[nz]; nb++) { + R[i] = tess->ra[nz][nb]; + D[i] = tess->dec[nz][nb]; + Xo[i] = tess->Xo[nz][nb]; + Yo[i] = tess->Yo[nz][nb]; + dX[i] = tess->dX[nz][nb]; + dY[i] = tess->dY[nz][nb]; + zone[i] = nz; + band[i] = nb; + index[i] = i; // or tess->cells[nz][nb] ? + memcpy(&name[i*BOUNDARY_TESS_NAME_LENGTH], tess->name[nz][nb], BOUNDARY_TESS_NAME_LENGTH); + 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, "X_CENT", Xo, Ncell); + gfits_set_bintable_column (&theader, &ftable, "Y_CENT", Yo, Ncell); + gfits_set_bintable_column (&theader, &ftable, "X_GRID", dX, Ncell); + gfits_set_bintable_column (&theader, &ftable, "Y_GRID", dY, Ncell); + gfits_set_bintable_column (&theader, &ftable, "NAME", name, Ncell); + + free (R ); + free (D ); + free (zone ); + free (band ); + free (index ); + free (Xo ); + free (Yo ); + free (dX ); + free (dY ); + free (name ); + + gfits_fwrite_Theader (f, &theader); + gfits_fwrite_table (f, &ftable); + gfits_free_header (&theader); + gfits_free_table (&ftable); + } + return TRUE; +} + +