Index: /branches/eam_branches/ipp-20120627/Ohana/src/addstar/include/skycells.h =================================================================== --- /branches/eam_branches/ipp-20120627/Ohana/src/addstar/include/skycells.h (revision 34209) +++ /branches/eam_branches/ipp-20120627/Ohana/src/addstar/include/skycells.h (revision 34210) @@ -12,5 +12,5 @@ # include -enum {SQUARES, TRIANGLES, LOCAL, RINGS}; +enum {SQUARES, TRIANGLES, LOCAL, RINGS, TAMAS}; enum {TETRAHEDRON, CUBE, OCTOHEDRON, DODECAHEDRON, ICOSAHEDRON}; @@ -95,4 +95,5 @@ int sky_tessellation_squares PROTO((FITS_DB *db, int level, int Nmax)); int sky_tessellation_rings PROTO((FITS_DB *db, int level, int Nmax)); +int sky_tessellation_tamas PROTO((FITS_DB *db, int level, int Nmax)); int sky_triangle_to_image PROTO((Image *image, SkyTriangle *triangle)); @@ -104,4 +105,5 @@ SkyRectangle *sky_rectangle_ring PROTO((float dec, float dDEC, int *nring, char *format)); +SkyRectangle *sky_rectangle_tamas PROTO((double *Dec, double dm, double halfa, double halftheta, int *nring, char *format)); SkyTriangle *sky_divide_triangles PROTO((SkyTriangle *in, int *ntriangles)); 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 34210) @@ -1,15 +1,15 @@ # include "addstar.h" -# define GET_COLUMN(OUT,NAME,TYPE) \ +# 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"); -// outstanding issues: -// create the ra,dec,name arrays in the original tree -// how do we deal with the name array in gfits? +# 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"); BoundaryTree *BoundaryTreeLoad(char *filename) { - int i, Ncol; + int i, j, nz, nb, Ncol; off_t Nrow; char type[16]; @@ -45,5 +45,4 @@ 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; @@ -58,8 +57,7 @@ // 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); + GET_COLUMN_RAW(tree->Nband, "NBAND", int); + GET_COLUMN_RAW(tree->RA_origin, "RA_ORIGIN", double); + GET_COLUMN_RAW(tree->RA_offset, "RA_OFFSET", double); gfits_free_header (&theader); gfits_free_table (&ftable); @@ -69,13 +67,16 @@ // allocate the storage arrays - ALLOCATE (tree->ra, double *, tree->Nzone); + 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]); + 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]); + for (j = 0; j < tree->Nband[i]; j++) { + ALLOCATE (tree->name[i][j], char, BOUNDARY_TREE_NAME_LENGTH); + } } } @@ -83,5 +84,4 @@ /*** cell information table ***/ { - // load data for this header if (!gfits_load_header (f, &theader)) goto escape; @@ -91,10 +91,10 @@ // 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? + 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(name, "NAME", char); // XXX how is this done? gfits_free_header (&theader); gfits_free_table (&ftable); @@ -107,6 +107,6 @@ tree->ra[nz][nb] = R[i]; tree->dec[nz][nb] = D[i]; - tree->cells[nz][nb] = i; // XXX ? - tree->name[nz][nb] = name[i]; + tree->cell[nz][nb] = i; // XXX ? + memcpy(tree->name[nz][nb], &name[i*BOUNDARY_TREE_NAME_LENGTH], BOUNDARY_TREE_NAME_LENGTH); } @@ -139,5 +139,5 @@ int BoundaryTreeSave(char *filename, BoundaryTree *tree) { - int i; + int i, nz, nb; Header header; Header theader; @@ -152,5 +152,5 @@ FILE *f = fopen (filename, "w"); if (!f) { - fprintf (stderr, "ERROR: cannot open image subset file for output %s\n", filename); + fprintf (stderr, "ERROR: cannot open boundary tree file for output %s\n", filename); return FALSE; } @@ -159,5 +159,4 @@ 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); @@ -170,8 +169,8 @@ 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, "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, "J", "NBAND", "number of cells in each zone", NULL, 1.0, 0.0); // generate the output array that carries the data @@ -179,21 +178,20 @@ // create intermediate storage arrays - int *zone ; ALLOCATE (zone , int , tree->Nzone); + int *zone = NULL; ALLOCATE (zone, int, tree->Nzone); // assign the storage arrays for (i = 0; i < tree->Nzone; i++) { - zone[i] = 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, "NBAND", tree->Nband, 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_fwrite_table (f, &ftable); gfits_free_header (&theader); gfits_free_table (&ftable); @@ -204,10 +202,12 @@ 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); + char fmt[16]; + snprintf (fmt, 16, "%dA", BOUNDARY_TREE_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, fmt, "NAME", "cell name", "none", 1.0, 0.0); // generate the output array that carries the data @@ -216,5 +216,5 @@ int Ncell = 0; for (i = 0; i < tree->Nzone; i++) { - Ncell += tree->Nband[i]; + Ncell += tree->Nband[i]; } @@ -226,26 +226,30 @@ int *band ; ALLOCATE (band, int, Ncell); int *index ; ALLOCATE (index, int, Ncell); - char **name ; ALLOCATE (name, char *, Ncell); + char *name ; ALLOCATE (name, char, Ncell*BOUNDARY_TREE_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 < 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++; + 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] ? + memcpy(&name[i*BOUNDARY_TREE_NAME_LENGTH], tree->name[nz][nb], BOUNDARY_TREE_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, "NAME", name, Ncell); + 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 ); @@ -263,2 +267,25 @@ return TRUE; } + +// the boundary tree... +// given an (ra,dec) pair, find the containing projection cell. + +int BoundaryTreeCellCoords (BoundaryTree *tree, int *zone, int *band, double ra, double dec) { + + // first, find the containing zone + + // if we know dDEC, we can get the bin instantly: + *zone = (dec - tree->DEC_origin) / tree->DEC_offset; + + if (*zone < 0) return FALSE; + if (*zone >= tree->Nzone) return FALSE; + + // now select the RA bin for that zone + *band = (ra - tree->RA_origin[*zone]) / tree->RA_offset[*zone]; + + if (*band < 0) return FALSE; + if (*band >= tree->Nband[*zone]) *band = 0; + + return TRUE; +} + Index: /branches/eam_branches/ipp-20120627/Ohana/src/addstar/src/args_skycells.c =================================================================== --- /branches/eam_branches/ipp-20120627/Ohana/src/addstar/src/args_skycells.c (revision 34209) +++ /branches/eam_branches/ipp-20120627/Ohana/src/addstar/src/args_skycells.c (revision 34210) @@ -37,4 +37,7 @@ if (!strcasecmp (argv[N], "rings")) { MODE = RINGS; + } + if (!strcasecmp (argv[N], "tamas")) { + MODE = TAMAS; } remove_argument (N, &argc, argv); @@ -179,4 +182,16 @@ } remove_argument (N, &argc, argv); + } + if (MODE == TAMAS) { + CELLSIZE = 3.955; + if ((N = get_argument (argc, argv, "-cellsize"))) { + remove_argument (N, &argc, argv); + CELLSIZE = strtod (argv[N], &ptr); + if ((*ptr != 0) || (CELLSIZE < 0.0)) { + fprintf (stderr, "-cellsize requires a floating-point argument\n"); + help (); + } + remove_argument (N, &argc, argv); + } } Index: /branches/eam_branches/ipp-20120627/Ohana/src/addstar/src/findskycell.c =================================================================== --- /branches/eam_branches/ipp-20120627/Ohana/src/addstar/src/findskycell.c (revision 34209) +++ /branches/eam_branches/ipp-20120627/Ohana/src/addstar/src/findskycell.c (revision 34210) @@ -1,3 +1,3 @@ -# include "skycells.h" +# include "addstar.h" static double RA_offset_RINGS_V3[] = {360.000000, 40.000000, 24.000000, 17.142857, 13.333333, 10.909091, 9.230769, 8.000000, 7.200000, @@ -12,31 +12,15 @@ // in an even more specific case, RA[i,zone] = RA_origin[zone] + RA_offset[zone] -typedef struct { - int FixedGridDEC; // is the DEC sequence linear? - int FixedGridRA; // in the RA sequence in a zone linear? - - double DEC_origin; - double DEC_offset; - - int Nzone; - double *RA_origin; - double *RA_offset; - - int *Nband; - int *NBAND; - - int **cells; -} BoundaryTree; - -enum {NONE, MAKE_TREE, USE_TREE}; +enum {TREE_NONE, TREE_MAKE, TREE_USE}; void usage (void) { fprintf (stderr, "USAGE: findcell -mktree (tree) (catdir)\n"); - fprintf (stderr, "USAGE: findcell -tree (tree) ra dec\n"); + fprintf (stderr, "USAGE: findcell -tree (tree) (datafile)\n"); + fprintf (stderr, " (datafile) should contain a list of RA,DEC pairs\n"); exit (2); } int mktree (char *treefile, char *catdir); -int find_primary_cell_coords (BoundaryTree *tree, int *zone, int *band, double ra, double dec); +int apply_tree (char *treefile, char *datafile); int main (int argc, char **argv) { @@ -55,7 +39,7 @@ /* extra error messages */ - MODE = NONE; + int MODE = TREE_NONE; if ((N = get_argument (argc, argv, "-mktree"))) { - MODE = MAKE_TREE; + MODE = TREE_MAKE; remove_argument (N, &argc, argv); treefile = strcreate (argv[N]); @@ -63,5 +47,5 @@ } if ((N = get_argument (argc, argv, "-tree"))) { - MODE = USE_TREE; + MODE = TREE_USE; remove_argument (N, &argc, argv); treefile = strcreate (argv[N]); @@ -73,10 +57,10 @@ if (!MODE) usage(); - if (MODE == MAKE_TREE) { + if (MODE == TREE_MAKE) { mktree (treefile, argv[1]); exit (0); } - // mktree (treefile, argv[1]); + apply_tree (treefile, argv[1]); exit (0); } @@ -131,15 +115,36 @@ ALLOCATE (tree.RA_origin, double, tree.Nzone); ALLOCATE (tree.RA_offset, double, tree.Nzone); - ALLOCATE (tree.cells, int *, tree.Nzone); + + ALLOCATE (tree.ra, double *, tree.Nzone); + ALLOCATE (tree.dec, double *, 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] = 0; + 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]; - ALLOCATE (tree.cells[zone], int, tree.NBAND[zone]); + ALLOCATE (tree.ra[zone], double, tree.NBAND[zone]); + ALLOCATE (tree.dec[zone], double, 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.cells[zone][band] = -1; + tree.ra[zone][band] = NAN; + tree.dec[zone][band] = NAN; + tree.cell[zone][band] = -1; + ALLOCATE (tree.name[zone][band], char, BOUNDARY_TREE_NAME_LENGTH); } } @@ -151,5 +156,5 @@ XY_to_RD (&ra, &dec, x, y, &image[i].coords); - if (!find_primary_cell_coords (&tree, &zone, &band, ra, dec)) { + if (!BoundaryTreeCellCoords (&tree, &zone, &band, ra, dec)) { fprintf (stderr, "mismatch!\n"); continue; @@ -159,12 +164,20 @@ if (band >= tree.NBAND[zone]) { int start = tree.NBAND[zone]; - tree.Nband[zone] = band + 1; tree.NBAND[zone] = band + 10; - REALLOCATE (tree.cells[zone], int, tree.NBAND[zone]); + REALLOCATE (tree.ra[zone], double, tree.NBAND[zone]); + REALLOCATE (tree.dec[zone], double, 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.cells[zone][j] = -1; + tree.ra[zone][j] = NAN; + tree.dec[zone][j] = NAN; + tree.cell[zone][j] = -1; + ALLOCATE (tree.name[zone][j], char, BOUNDARY_TREE_NAME_LENGTH); } } - tree.cells[zone][band] = i; + tree.ra[zone][band] = ra; + tree.dec[zone][band] = dec; + tree.cell[zone][band] = i; + memcpy (tree.name[zone][band], image[i].name, BOUNDARY_TREE_NAME_LENGTH); } @@ -175,5 +188,5 @@ for (band = 0; band < tree.NBAND[zone]; band++) { // all cells should be filled - if (tree.cells[zone][band] < 0) { + if (tree.cell[zone][band] < 0) { if (!found_last) { found_last = TRUE; @@ -205,5 +218,5 @@ ra = 360.0 * drand48(); dec = 180.0 * drand48() - 90.0; - if (!find_primary_cell_coords (&tree, &zone, &band, ra, dec)) { + if (!BoundaryTreeCellCoords (&tree, &zone, &band, ra, dec)) { fprintf (stderr, "failure for %f,%f\n", ra, dec); } @@ -216,23 +229,31 @@ } -// the boundary tree... -// given an (ra,dec) pair, find the containing projection cell. - -int find_primary_cell_coords (BoundaryTree *tree, int *zone, int *band, double ra, double dec) { - - // first, find the containing zone - - // if we know dDEC, we can get the bin instantly: - *zone = (dec - tree->DEC_origin) / tree->DEC_offset; - - if (*zone < 0) return FALSE; - if (*zone >= tree->Nzone) return FALSE; - - // now select the RA bin for that zone - *band = (ra - tree->RA_origin[*zone]) / tree->RA_offset[*zone]; - - if (*band < 0) return FALSE; - // if (*band >= tree->Nband[*zone]) return FALSE; - - return TRUE; -} +int apply_tree (char *treefile, char *datafile) { + + BoundaryTree *tree = BoundaryTreeLoad (treefile); + if (!tree) { + fprintf (stderr, "error loading boundary tree 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 zone, band; + if (!BoundaryTreeCellCoords (tree, &zone, &band, ra, dec)) { + fprintf (stderr, "error finding cell for %f,%f\n", ra, dec); + continue; + } + + fprintf (stdout, "%10.6f %10.6f %3d %3d %s\n", ra, dec, zone, band, tree->name[zone][band]); + } + + exit (0); +} Index: /branches/eam_branches/ipp-20120627/Ohana/src/addstar/src/sky_tessalation.c =================================================================== --- /branches/eam_branches/ipp-20120627/Ohana/src/addstar/src/sky_tessalation.c (revision 34209) +++ /branches/eam_branches/ipp-20120627/Ohana/src/addstar/src/sky_tessalation.c (revision 34210) @@ -23,4 +23,7 @@ sky_tessellation_rings (db, level, Nmax); return TRUE; + case TAMAS: + sky_tessellation_tamas (db, level, Nmax); + return TRUE; default: break; @@ -284,4 +287,72 @@ free (ring); free (image); + } + return (TRUE); +} + +// the RINGS tessellation uses the declination zones proposed by Tamas Budavari, +// based on code supplied by Tamas 2012.07.23 +int sky_tessellation_tamas (FITS_DB *db, int level, int Nmax) { + + int j, nDEC, Nimage, Nring, Ntotal, Ndigit; + double dec, dDEC; + SkyRectangle *ring; + Image *image; + char format[16]; + + // The tessellation has one input parameter: the approximate cell size. Starting with + // the cell size, determine the optimal projection cell height (dDEC) that results in an + // integer number of dec zones between -90 and +90 + + // in fact, we place a single image on each pole, so the real range of dec is 180.0 - CELLSIZE: + + nDEC = (180.0 - CELLSIZE) / CELLSIZE; + dDEC = (180.0 - CELLSIZE) / nDEC; + nDEC += 2; + + // how many total projection cells for this realization? divide sky area by cell area: + // this is used to set the number of digits, so it does not need to be very accurate... + Ntotal = 41254.2 / (dDEC*dDEC); + Ndigit = (int)(log10(Ntotal)) + 1 ; + snprintf (format, 16, "skycell.%%0%dd", Ndigit); + + double d2r = M_PI / 180; // is RAD_DEG + + // parameter 'a' is the cell size in degrees + double adeg = 3.955; + + // half of 'a' in radians and its atan + double halfa = adeg / 2 * d2r; + double halftheta = atan(halfa); + + // loop init + dec = 0; // starting Decl. - could change this... + + while (dec < M_PI / 2 - halftheta) { + double dm = dec - halftheta; // eq.5 + if (dec == 0) dm = 0; // initial + + // dec is modified by the call below + ring = sky_rectangle_tamas (&dec, dm, halfa, halftheta, &Nring, format); + if (!ring) continue; + + // subdivide each image (Nx x Ny subcells) + Nimage = NX_SUB*NY_SUB*Nring; + ALLOCATE (image, Image, Nimage); + for (j = 0; j < Nring; j++) { + // convert the SkyRectangles to Images for output + sky_subdivide_image (&image[j*NX_SUB*NY_SUB], &ring[j], NX_SUB, NY_SUB); + // printf("%s %8.2f %8.2f\n", ring[j].name, ring[j].coords.crval1, ring[j].coords.crval2); + } + + /* add the new images and save */ + dvo_image_addrows (db, image, Nimage); + SetProtect (TRUE); + dvo_image_update (db, VERBOSE); + SetProtect (FALSE); + dvo_image_clear_vtable (db); + + free (ring); + free (image); } return (TRUE); @@ -666,4 +737,79 @@ } +// define the parameters of a projection centers for this ring +// dec : ~ center of ring in Dec +// dDEC : approximate height +// nring : number of cells generated for this ring +// format : guide to generate the filenames (c-type string format) +SkyRectangle *sky_rectangle_tamas (double *Dec, double dm, double halfa, double halftheta, int *nring, char *format) { + + static int Nname = 0; + int i, j, NX, NY; + SkyRectangle *ring; + + double d2r = M_PI / 180; // is RAD_DEG + double dec = *Dec; + + int nRA = (int)ceil(M_PI * cos(dm) / halftheta); // eq.6 + double dRA = 2 * M_PI / nRA; // eq.7 + double dp = atan(tan(dec + halftheta) * cos(dRA / 2)); // eq.9 + + if (dec == 0.0) { + ALLOCATE (ring, SkyRectangle, nRA); + } else { + ALLOCATE (ring, SkyRectangle, 2*nRA); + } + + for (i = 0; i < nRA; i++) { + // R.A. can use different phase per ring + double ra = i * dRA; // + phase (watch wraparound) + + int npass = (dec == 0.0) ? 1 : 2; + for (j = 0; j < npass; j++) { + + int N = j*nRA + i; + + memset (&ring[N], 0, sizeof(SkyRectangle)); + memset (&ring[N].coords, 0, sizeof(Coords)); + + ring[N].coords.crval1 = ra / d2r; + ring[N].coords.crval2 = (j == 0) ? dec / d2r : -dec / d2r; + + printf(" \t %d %25.20f %25.20f\n", i, ring[N].coords.crval2, ring[N].coords.crval1); + + ring[N].coords.pc1_1 = +1.0 * X_PARITY; + ring[N].coords.pc1_2 = +0.0; + ring[N].coords.pc2_1 = -0.0; + ring[N].coords.pc2_2 = +1.0; + + // range values are in projected degrees + NX = cos(dec - halftheta) * dRA * 3600.0 / SCALE / d2r; + NY = 2 * halftheta * 3600.0 / SCALE / d2r; + + // crpix1,crpix2 is the projection center + ring[N].coords.crpix1 = 0.5*NX; + ring[N].coords.crpix2 = 0.5*NY; + + ring[N].coords.cdelt1 = SCALE / 3600.0; + ring[N].coords.cdelt2 = SCALE / 3600.0; + + strcpy (ring[N].coords.ctype, "DEC--TAN"); + + ring[N].NX = NX*(1.0 + PADDING); + ring[N].NY = NY*(1.0 + PADDING); + ring[N].photcode = 1; // this needs to be set more sensibly + + snprintf (ring[N].name, DVO_IMAGE_NAME_LEN, format, Nname); + Nname++; + } + } + + // advance to next ring + *Dec = halftheta + dp; + + *nring = (dec == 0.0) ? nRA : 2*nRA; + return ring; +} + // an allocated image set is supplied, we fill in the values int sky_subdivide_image (Image *output, SkyRectangle *input, int Nx, int Ny) { @@ -682,6 +828,10 @@ } - Ndigit = (int)(log10(Nx*Ny)) + 1 ; - snprintf (format, 24, "%s.%%0%dd", input[0].name, Ndigit); + if (Nx * Ny > 1) { + Ndigit = (int)(log10(Nx*Ny)) + 1 ; + snprintf (format, 24, "%s.%%0%dd", input[0].name, Ndigit); + } else { + snprintf (format, 24, "%s", input[0].name); + } // if requested extend, the skycell boundaries so that skycells overlap @@ -696,5 +846,10 @@ memcpy (&output[N].coords, &input[0].coords, sizeof(Coords)); - snprintf (output[N].name, DVO_IMAGE_NAME_LEN, format, N); + if (Nx + Ny > 1) { + snprintf (output[N].name, DVO_IMAGE_NAME_LEN, format, N); + } else { + snprintf (output[N].name, DVO_IMAGE_NAME_LEN, "%s", format); + } + output[N].NX = NX + 2 * pad_x; output[N].NY = NY + 2 * pad_y; Index: /branches/eam_branches/ipp-20120627/Ohana/src/libdvo/Makefile =================================================================== --- /branches/eam_branches/ipp-20120627/Ohana/src/libdvo/Makefile (revision 34209) +++ /branches/eam_branches/ipp-20120627/Ohana/src/libdvo/Makefile (revision 34210) @@ -102,5 +102,6 @@ $(SRC)/db_utils.$(ARCH).o \ $(SRC)/convert.$(ARCH).o \ -$(SRC)/HostTable.$(ARCH).o +$(SRC)/HostTable.$(ARCH).o \ +$(SRC)/BoundaryTree.$(ARCH).o Index: /branches/eam_branches/ipp-20120627/Ohana/src/libdvo/include/dvo.h =================================================================== --- /branches/eam_branches/ipp-20120627/Ohana/src/libdvo/include/dvo.h (revision 34209) +++ /branches/eam_branches/ipp-20120627/Ohana/src/libdvo/include/dvo.h (revision 34210) @@ -293,4 +293,26 @@ } DVOTinyValueMode; +# define BOUNDARY_TREE_NAME_LENGTH 128 + +typedef struct { + int FixedGridDEC; // is the DEC sequence linear? + int FixedGridRA; // in the RA sequence in a zone linear? + + double DEC_origin; + double DEC_offset; + + int Nzone; + double *RA_origin; + double *RA_offset; + + int *Nband; + int *NBAND; + + double **ra; + double **dec; + int **cell; + char ***name; +} BoundaryTree; + // XXX DROP? // a reduced-subset structure for relastro // XXX DROP? typedef struct { @@ -662,4 +684,7 @@ int free_tiny_values (Catalog *catalog); +int BoundaryTreeCellCoords (BoundaryTree *tree, int *zone, int *band, double ra, double dec); +int BoundaryTreeSave(char *filename, BoundaryTree *tree); +BoundaryTree *BoundaryTreeLoad(char *filename); + # endif // DVO_H - Index: /branches/eam_branches/ipp-20120627/Ohana/src/libdvo/src/BoundaryTree.c =================================================================== --- /branches/eam_branches/ipp-20120627/Ohana/src/libdvo/src/BoundaryTree.c (revision 34210) +++ /branches/eam_branches/ipp-20120627/Ohana/src/libdvo/src/BoundaryTree.c (revision 34210) @@ -0,0 +1,293 @@ +# 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 + +BoundaryTree *BoundaryTreeLoad(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; + 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 (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 (tree, BoundaryTree, 1); + + gfits_scan (&header, "DEC_ORI", "%lf", 1, &tree->DEC_origin); + gfits_scan (&header, "DEC_OFF", "%lf", 1, &tree->DEC_offset); + + 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(tree->Nband, "NBAND", int); + GET_COLUMN_RAW(tree->RA_origin, "RA_ORIGIN", double); + GET_COLUMN_RAW(tree->RA_offset, "RA_OFFSET", double); + 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]); + for (j = 0; j < tree->Nband[i]; j++) { + ALLOCATE (tree->name[i][j], char, BOUNDARY_TREE_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(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]; + tree->ra[nz][nb] = R[i]; + tree->dec[nz][nb] = D[i]; + tree->cell[nz][nb] = i; // XXX ? + memcpy(tree->name[nz][nb], &name[i*BOUNDARY_TREE_NAME_LENGTH], BOUNDARY_TREE_NAME_LENGTH); + } + + 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, 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 tree 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_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); + + // generate the output array that carries the data + gfits_create_table (&theader, &ftable); + + // create intermediate storage arrays + int *zone = NULL; 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, "NBAND", tree->Nband, 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); + 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_TREE_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, 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 < 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*BOUNDARY_TREE_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 < 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] ? + memcpy(&name[i*BOUNDARY_TREE_NAME_LENGTH], tree->name[nz][nb], BOUNDARY_TREE_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, "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; +} + +// the boundary tree... +// given an (ra,dec) pair, find the containing projection cell. + +int BoundaryTreeCellCoords (BoundaryTree *tree, int *zone, int *band, double ra, double dec) { + + // first, find the containing zone + + // if we know dDEC, we can get the bin instantly: + *zone = (dec - tree->DEC_origin) / tree->DEC_offset; + + if (*zone < 0) return FALSE; + if (*zone >= tree->Nzone) return FALSE; + + // now select the RA bin for that zone + *band = (ra - tree->RA_origin[*zone]) / tree->RA_offset[*zone]; + + if (*band < 0) return FALSE; + if (*band >= tree->Nband[*zone]) *band = 0; + + return TRUE; +} + Index: /branches/eam_branches/ipp-20120627/Ohana/src/relphot/Makefile =================================================================== --- /branches/eam_branches/ipp-20120627/Ohana/src/relphot/Makefile (revision 34209) +++ /branches/eam_branches/ipp-20120627/Ohana/src/relphot/Makefile (revision 34210) @@ -49,4 +49,5 @@ $(SRC)/setExclusions.$(ARCH).o \ $(SRC)/setMrelFinal.$(ARCH).o \ +$(SRC)/BoundaryTreeOps.$(ARCH).o \ $(SRC)/write_coords.$(ARCH).o @@ -77,4 +78,5 @@ $(SRC)/setExclusions.$(ARCH).o \ $(SRC)/setMrelFinal.$(ARCH).o \ +$(SRC)/BoundaryTreeOps.$(ARCH).o \ $(SRC)/write_coords.$(ARCH).o Index: /branches/eam_branches/ipp-20120627/Ohana/src/relphot/include/relphot.h =================================================================== --- /branches/eam_branches/ipp-20120627/Ohana/src/relphot/include/relphot.h (revision 34209) +++ /branches/eam_branches/ipp-20120627/Ohana/src/relphot/include/relphot.h (revision 34210) @@ -111,4 +111,6 @@ char *BCATALOG; ModeType MODE; + +char *BOUNDARY_TREE; double MAG_LIM; @@ -341,2 +343,6 @@ int client_logger_message (char *format,...); +int MatchImageName (off_t meas, int cat, char *name); + +int load_tree (char *treefile); +int BoundaryTreePrimaryCell (char *primaryCellName, double ra, double dec); Index: /branches/eam_branches/ipp-20120627/Ohana/src/relphot/src/BoundaryTreeOps.c =================================================================== --- /branches/eam_branches/ipp-20120627/Ohana/src/relphot/src/BoundaryTreeOps.c (revision 34210) +++ /branches/eam_branches/ipp-20120627/Ohana/src/relphot/src/BoundaryTreeOps.c (revision 34210) @@ -0,0 +1,36 @@ +# include "relphot.h" + +// XXX for the moment, only load one boundary tree at a time +// XXX in fact, only allow RINGS.V3... + +static BoundaryTree *tree = NULL; + +int BoundaryTreePrimaryCell (char *primaryCellName, double ra, double dec) { + + int zone, band; + + if (!primaryCellName) return FALSE; + + primaryCellName[0] = 0; + + if (!tree) return FALSE; + + if (!BoundaryTreeCellCoords (tree, &zone, &band, ra, dec)) { + fprintf (stderr, "mismatch!\n"); + return FALSE; + } + + snprintf (primaryCellName, DVO_MAX_PATH, "RINGS.V3.%s", tree->name[zone][band]); + return TRUE; +} + +int load_tree (char *treefile) { + + tree = BoundaryTreeLoad (treefile); + if (!tree) { + fprintf (stderr, "failed to load boundary tree %s\n", treefile); + exit (2); + } + + return TRUE; +} Index: /branches/eam_branches/ipp-20120627/Ohana/src/relphot/src/ImageOps.c =================================================================== --- /branches/eam_branches/ipp-20120627/Ohana/src/relphot/src/ImageOps.c (revision 34209) +++ /branches/eam_branches/ipp-20120627/Ohana/src/relphot/src/ImageOps.c (revision 34210) @@ -385,4 +385,25 @@ distance = hypot (measure[0].Xccd - Xcenter, measure[0].Yccd - Ycenter); return (distance); +} + +// returns image.Mcal - ff(x,y) +int MatchImageName (off_t meas, int cat, char *name) { + + off_t i; + + if (!name) return FALSE; + if (!name[0]) return FALSE; + + if (!MeasureToImage) return FALSE; + + i = MeasureToImage[cat][meas]; + if (i == -1) return FALSE; + + // this is a bit crude: stack image names are of the form + // RINGS.V3.skycell.1495.027.sky.191211.stk.988232.cmf + // the primaryCell has a name of the form RINGS.V3.skycell.1495 + + if (!strncmp(image[i].name, name, strlen(name))) return TRUE; + return FALSE; } Index: /branches/eam_branches/ipp-20120627/Ohana/src/relphot/src/StarOps.c =================================================================== --- /branches/eam_branches/ipp-20120627/Ohana/src/relphot/src/StarOps.c (revision 34209) +++ /branches/eam_branches/ipp-20120627/Ohana/src/relphot/src/StarOps.c (revision 34210) @@ -306,5 +306,5 @@ int setMrel_catalog (Catalog *catalog, int Nc, int pass, FlatCorrectionTable *flatcorr, SetMrelInfo *results, int Nsecfilt) { - off_t j, k, m; + off_t j, k, m, ID; int N; float Msys, Mcal, Mmos, Mgrid; @@ -326,4 +326,7 @@ int isSetMrelFinal = (pass >= 0); + char *primaryCell = NULL; + ALLOCATE (primaryCell, char, DVO_MAX_PATH); + for (j = 0; j < catalog[Nc].Naverage; j++) { // XXX accumulate all secfilt values in a single pass? @@ -334,4 +337,6 @@ print_measure_set (&catalog[Nc].average[j], &catalog[Nc].secfilt[j*Nsecfilt], catalog[Nc].measure); } + + BoundaryTreePrimaryCell(primaryCell, catalog[Nc].average[j].R, catalog[Nc].average[j].D); int GoodPS1 = FALSE; @@ -361,10 +366,14 @@ int haveStack = FALSE; - float stackFluxPSF = NAN; - float stackdFluxPSF = NAN; - float stackFluxKron = NAN; - float stackdFluxKron = NAN; + // need to find the measurement closest to the center of its skycell, as well as the + // closest for the subset of primary projection cells + float stackCenterOffsetMin = 1e9; - unsigned int stackImageIDmin; + int stackCenterIDmin = -1; + off_t stackCenterMeasureMin = -1; + + float stackPrimaryOffsetMin = 1e9; + int stackPrimaryIDmin = -1; + off_t stackPrimaryMeasureMin = -1; int forceSynth = FALSE; @@ -436,47 +445,20 @@ unsigned int stackImageID; + + // 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 (MatchImageName (m, Nc, primaryCell)) { + float stackPrimaryOffset = getCenterOffset (m, Nc, &catalog[Nc].measure[m], &stackImageID); + if (stackPrimaryOffset < stackPrimaryOffsetMin) { + stackPrimaryIDmin = stackImageID; + stackPrimaryMeasureMin = m; + } + } + + // get the center distance for the generic case: float stackCenterOffset = getCenterOffset (m, Nc, &catalog[Nc].measure[m], &stackImageID); if (stackCenterOffset < stackCenterOffsetMin) { - stackImageIDmin = stackImageID; - - float zp = Mcal + Mmos + Mgrid + PhotZeroPoint (&catalog[Nc].measure[m], &catalog[Nc].average[j], &catalog[Nc].secfilt[j*Nsecfilt]); - - // flux_cgs : erg sec^1 cm^-2 Hz^-1 - // mag_inst : -2.5 log (cts/sec) - // mag_inst : -2.5 log (flux_inst) - // flux_inst = ten(-0.4*mag_inst) - - // mag_AB = -2.5 log (flux_cgs) - 48.6 (~by definition) [~Vega flux in V-band] - // flux_cgs = ten(-0.4*(mag_AB + 48.6)) - - // flux_AB : ten(-0.4*mag_AB) - - // flux_cgs = ten(-0.4*48.6) * flux_AB - // flux_AB = ten(+0.4*48.6) * flux_cgs - - // flux_Jy : flux_cgs * 10^23 - - // flux_AB = ten(+0.4*48.6) * ten(-23) * flux_Jy - - // mag_AB = mag_inst + ZP - - // flux_inst = ten(-0.4*(mag_AB - ZP)) = ten(0.4*ZP) * flux_AB - - // flux_AB = flux_inst * ten(-0.4*ZP) - - // flux_inst * ten(-0.4*ZP) = ten(+0.4*48.6 - 23) * flux_Jy - - // flux_inst = flux_Jy * ten(0.4*ZP + 0.4*48.6 - 23) - // flux_inst = flux_Jy * ten(0.4*ZP - 3.56) - // flux_Jy = flux_inst * ten(-0.4*ZP + 3.56) - - // zpFactor to go from instrumental flux to Janskies - float zpFactor = pow(10.0, -0.4*zp + 3.56); - - // need to put in AB mag factor to get to Janskies (or uJy?) - stackFluxPSF = zpFactor * catalog[Nc].measure[m].FluxPSF; - stackdFluxPSF = zpFactor * catalog[Nc].measure[m].dFluxPSF; - stackFluxKron = zpFactor * catalog[Nc].measure[m].FluxKron; - stackdFluxKron = zpFactor * catalog[Nc].measure[m].dFluxKron; + stackCenterIDmin = stackImageID; + stackCenterMeasureMin = m; } } @@ -619,9 +601,49 @@ if (haveStack) { - catalog[Nc].secfilt[Nsecfilt*j+Nsec].FluxPSF = stackFluxPSF; - catalog[Nc].secfilt[Nsecfilt*j+Nsec].dFluxPSF = stackdFluxPSF; - catalog[Nc].secfilt[Nsecfilt*j+Nsec].FluxKron = stackFluxKron; - catalog[Nc].secfilt[Nsecfilt*j+Nsec].dFluxKron = stackdFluxKron; - catalog[Nc].secfilt[Nsecfilt*j+Nsec].stackID = stackImageIDmin; + m = (stackPrimaryMeasureMin >= 0) ? stackPrimaryMeasureMin : stackCenterMeasureMin; + ID = (stackPrimaryMeasureMin >= 0) ? stackPrimaryIDmin : stackCenterIDmin; + + // get the zero point for the selected image + float zp = Mcal + Mmos + Mgrid + PhotZeroPoint (&catalog[Nc].measure[m], &catalog[Nc].average[j], &catalog[Nc].secfilt[j*Nsecfilt]); + + // flux_cgs : erg sec^1 cm^-2 Hz^-1 + // mag_inst : -2.5 log (cts/sec) + // mag_inst : -2.5 log (flux_inst) + // flux_inst = ten(-0.4*mag_inst) + + // mag_AB = -2.5 log (flux_cgs) - 48.6 (~by definition) [~Vega flux in V-band] + // flux_cgs = ten(-0.4*(mag_AB + 48.6)) + + // flux_AB : ten(-0.4*mag_AB) + + // flux_cgs = ten(-0.4*48.6) * flux_AB + // flux_AB = ten(+0.4*48.6) * flux_cgs + + // flux_Jy : flux_cgs * 10^23 + + // flux_AB = ten(+0.4*48.6) * ten(-23) * flux_Jy + + // mag_AB = mag_inst + ZP + + // flux_inst = ten(-0.4*(mag_AB - ZP)) = ten(0.4*ZP) * flux_AB + + // flux_AB = flux_inst * ten(-0.4*ZP) + + // flux_inst * ten(-0.4*ZP) = ten(+0.4*48.6 - 23) * flux_Jy + + // flux_inst = flux_Jy * ten(0.4*ZP + 0.4*48.6 - 23) + // flux_inst = flux_Jy * ten(0.4*ZP - 3.56) + // flux_Jy = flux_inst * ten(-0.4*ZP + 3.56) + + // zpFactor to go from instrumental flux to Janskies + float zpFactor = pow(10.0, -0.4*zp + 3.56); + + // need to put in AB mag factor to get to Janskies (or uJy?) + catalog[Nc].secfilt[Nsecfilt*j+Nsec].FluxPSF = zpFactor * catalog[Nc].measure[m].FluxPSF; + catalog[Nc].secfilt[Nsecfilt*j+Nsec].dFluxPSF = zpFactor * catalog[Nc].measure[m].dFluxPSF; + catalog[Nc].secfilt[Nsecfilt*j+Nsec].FluxKron = zpFactor * catalog[Nc].measure[m].FluxKron; + catalog[Nc].secfilt[Nsecfilt*j+Nsec].dFluxKron = zpFactor * catalog[Nc].measure[m].dFluxKron; + + catalog[Nc].secfilt[Nsecfilt*j+Nsec].stackID = ID; } @@ -680,4 +702,5 @@ } } + if (primaryCell) free (primaryCell); return (TRUE); } Index: /branches/eam_branches/ipp-20120627/Ohana/src/relphot/src/args.c =================================================================== --- /branches/eam_branches/ipp-20120627/Ohana/src/relphot/src/args.c (revision 34209) +++ /branches/eam_branches/ipp-20120627/Ohana/src/relphot/src/args.c (revision 34210) @@ -197,4 +197,10 @@ } + if ((N = get_argument (argc, argv, "-boundary-tree"))) { + remove_argument (N, &argc, argv); + load_tree (argv[N]); + remove_argument (N, &argc, argv); + } + SHOW_PARAMS = FALSE; if ((N = get_argument (argc, argv, "-params"))) {