IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 34210


Ignore:
Timestamp:
Jul 26, 2012, 8:10:41 AM (14 years ago)
Author:
eugene
Message:

adding tamas rings, tools for making boundary tree and using it in psphot, etc

Location:
branches/eam_branches/ipp-20120627/Ohana/src
Files:
2 added
12 edited

Legend:

Unmodified
Added
Removed

  • branches/eam_branches/ipp-20120627/Ohana/src/addstar/include/skycells.h

    r33719 r34210  
    1212# include <glob.h>
    1313
    14 enum {SQUARES, TRIANGLES, LOCAL, RINGS};
     14enum {SQUARES, TRIANGLES, LOCAL, RINGS, TAMAS};
    1515enum {TETRAHEDRON, CUBE, OCTOHEDRON, DODECAHEDRON, ICOSAHEDRON};
    1616
    … …  
    9595int          sky_tessellation_squares       PROTO((FITS_DB *db, int level, int Nmax));
    9696int          sky_tessellation_rings         PROTO((FITS_DB *db, int level, int Nmax));
     97int          sky_tessellation_tamas         PROTO((FITS_DB *db, int level, int Nmax));
    9798
    9899int          sky_triangle_to_image          PROTO((Image *image, SkyTriangle *triangle));
    … …  
    104105
    105106SkyRectangle *sky_rectangle_ring            PROTO((float dec, float dDEC, int *nring, char *format));
     107SkyRectangle *sky_rectangle_tamas           PROTO((double *Dec, double dm, double halfa, double halftheta, int *nring, char *format));
    106108
    107109SkyTriangle *sky_divide_triangles           PROTO((SkyTriangle *in, int *ntriangles));

  • branches/eam_branches/ipp-20120627/Ohana/src/addstar/src/BoundaryTreeIO.c

    r34209 r34210  
    11# include "addstar.h"
    22
    3 # define GET_COLUMN(OUT,NAME,TYPE)                                      \
     3# define GET_COLUMN_NEW(OUT,NAME,TYPE)                                  \
    44  TYPE *OUT = gfits_get_bintable_column_data (&theader, &ftable, NAME, type, &Nrow, &Ncol); \
    55  myAssert (!strcmp(type, #TYPE), "wrong column type");
    66
    7 // outstanding issues:
    8 //   create the ra,dec,name arrays in the original tree
    9 //   how do we deal with the name array in gfits?
     7# define GET_COLUMN_RAW(OUT,NAME,TYPE)                                  \
     8  OUT = gfits_get_bintable_column_data (&theader, &ftable, NAME, type, &Nrow, &Ncol); \
     9  myAssert (!strcmp(type, #TYPE), "wrong column type");
    1010
    1111BoundaryTree *BoundaryTreeLoad(char *filename) {
    1212
    13   int i, Ncol;
     13  int i, j, nz, nb, Ncol;
    1414  off_t Nrow;
    1515  char type[16];
    … …  
    4545  gfits_scan (&header, "DEC_ORI", "%lf", 1, &tree->DEC_origin);
    4646  gfits_scan (&header, "DEC_OFF", "%lf", 1, &tree->DEC_offset);
    47   // gfits_scan (&header, "NZONE",   "%d",  1, &tree->Nzone); XXX not needed (redundant with table length)
    4847
    4948  ftable.header = &theader;
    … …  
    5857 
    5958    // need to create and assign to flat-field correction
    60     GET_COLUMN(ZONE,            "ZONE",         int);
    61     GET_COLUMN(tree->RA_origin, "RA_ORIGIN",    double);
    62     GET_COLUMN(tree->RA_offset, "RA_OFFSET",    double);
    63     GET_COLUMN(tree->Nband,     "NBAND",        int);
     59    GET_COLUMN_RAW(tree->Nband,     "NBAND",     int);
     60    GET_COLUMN_RAW(tree->RA_origin, "RA_ORIGIN", double);
     61    GET_COLUMN_RAW(tree->RA_offset, "RA_OFFSET", double);
    6462    gfits_free_header (&theader);
    6563    gfits_free_table  (&ftable);
    … …  
    6967
    7068    // allocate the storage arrays
    71     ALLOCATE (tree->ra,  double *, tree->Nzone);
     69    ALLOCATE (tree->ra,   double *, tree->Nzone);
    7270    ALLOCATE (tree->dec,  double *, tree->Nzone);
    7371    ALLOCATE (tree->cell, int *, tree->Nzone);
    7472    ALLOCATE (tree->name, char **, tree->Nzone);
    7573    for (i = 0; i < tree->Nzone; i++) {
    76       ALLOCATE (tree->ra[i],   double *, tree->Nband[i]);
    77       ALLOCATE (tree->dec[i],  double *, tree->Nband[i]);
    78       ALLOCATE (tree->cell[i], int *,    tree->Nband[i]);
    79       ALLOCATE (tree->name[i], char **,  tree->Nband[i]);
     74      ALLOCATE (tree->ra[i],   double, tree->Nband[i]);
     75      ALLOCATE (tree->dec[i],  double, tree->Nband[i]);
     76      ALLOCATE (tree->cell[i], int,    tree->Nband[i]);
     77      ALLOCATE (tree->name[i], char *, tree->Nband[i]);
     78      for (j = 0; j < tree->Nband[i]; j++) {
     79        ALLOCATE (tree->name[i][j], char, BOUNDARY_TREE_NAME_LENGTH);
     80      }
    8081    }
    8182  }
    … …  
    8384  /*** cell information table ***/
    8485  {
    85 
    8686    // load data for this header
    8787    if (!gfits_load_header (f, &theader)) goto escape;
    … …  
    9191 
    9292    // need to create and assign to flat-field correction
    93     GET_COLUMN(R,     "RA",        double);
    94     GET_COLUMN(D,     "DEC",      double);
    95     GET_COLUMN(zone,  "NMEAS",       int);
    96     GET_COLUMN(band,  "MEASURE_OFF", int);
    97     GET_COLUMN(index, "FLAGS",       int);
    98     GET_COLUMN(name,  "CAT_ID",      char); // XXX how is this done?
     93    GET_COLUMN_NEW(R,     "RA",          double);
     94    GET_COLUMN_NEW(D,     "DEC",        double);
     95    GET_COLUMN_NEW(zone,  "ZONE",        int);
     96    GET_COLUMN_NEW(band,  "BAND",        int);
     97    GET_COLUMN_NEW(index, "INDEX",       int);
     98    GET_COLUMN_NEW(name,  "NAME",        char); // XXX how is this done?
    9999    gfits_free_header (&theader);
    100100    gfits_free_table  (&ftable);
    … …  
    107107      tree->ra[nz][nb] = R[i];
    108108      tree->dec[nz][nb] = D[i];
    109       tree->cells[nz][nb] = i; // XXX ?
    110       tree->name[nz][nb] = name[i];
     109      tree->cell[nz][nb] = i; // XXX ?
     110      memcpy(tree->name[nz][nb], &name[i*BOUNDARY_TREE_NAME_LENGTH], BOUNDARY_TREE_NAME_LENGTH);
    111111    }
    112112
    … …  
    139139int BoundaryTreeSave(char *filename, BoundaryTree *tree) {
    140140
    141   int i;
     141  int i, nz, nb;
    142142  Header header;
    143143  Header theader;
    … …  
    152152  FILE *f = fopen (filename, "w");
    153153  if (!f) {
    154     fprintf (stderr, "ERROR: cannot open image subset file for output %s\n", filename);
     154    fprintf (stderr, "ERROR: cannot open boundary tree file for output %s\n", filename);
    155155    return FALSE;
    156156  }
    … …  
    159159  gfits_modify (&header, "DEC_ORI", "%lf", 1, tree->DEC_origin);
    160160  gfits_modify (&header, "DEC_OFF", "%lf", 1, tree->DEC_offset);
    161   gfits_modify (&header, "NZONE",   "%d",  1, tree->Nzone);
    162161
    163162  gfits_fwrite_header  (f, &header);
    … …  
    170169    gfits_create_table_header (&theader, "BINTABLE", "ZONE_DATA");
    171170
    172     gfits_define_bintable_column (&theader, "J", "ZONE",      "zone sequence number", NULL, 1.0, 0.0);
     171    gfits_define_bintable_column (&theader, "J", "ZONE",      "zone sequence number", "none", 1.0, 0.0);
     172    gfits_define_bintable_column (&theader, "J", "NBAND",     "number of cells in each zone", "none", 1.0, 0.0);
    173173    gfits_define_bintable_column (&theader, "D", "RA_ORIGIN", "origin of ra cell sequence", "degree", 1.0, 0.0);
    174174    gfits_define_bintable_column (&theader, "D", "RA_OFFSET", "offset per cell of ra cell sequence", "degree/cell", 1.0, 0.0);
    175     gfits_define_bintable_column (&theader, "J", "NBAND",     "number of cells in each zone", NULL, 1.0, 0.0);
    176175
    177176    // generate the output array that carries the data
    … …  
    179178
    180179    // create intermediate storage arrays
    181     int   *zone      ; ALLOCATE (zone     ,  int  , tree->Nzone);
     180    int *zone = NULL; ALLOCATE (zone,  int, tree->Nzone);
    182181
    183182    // assign the storage arrays
    184183    for (i = 0; i < tree->Nzone; i++) {
    185         zone[i] = i;
     184      zone[i] = i;
    186185    }
    187186
    188187    // add the columns to the output array
    189188    gfits_set_bintable_column (&theader, &ftable, "ZONE",       zone,            tree->Nzone);
     189    gfits_set_bintable_column (&theader, &ftable, "NBAND",      tree->Nband,     tree->Nzone);
    190190    gfits_set_bintable_column (&theader, &ftable, "RA_ORIGIN",  tree->RA_origin, tree->Nzone);
    191191    gfits_set_bintable_column (&theader, &ftable, "RA_OFFSET",  tree->RA_offset, tree->Nzone);
    192     gfits_set_bintable_column (&theader, &ftable, "NBAND",      tree->Nband,     tree->Nzone);
    193 
    194192    free (zone);
    195193
    196194    gfits_fwrite_Theader (f, &theader);
    197     gfits_fwrite_table  (f, &ftable);
     195    gfits_fwrite_table (f, &ftable);
    198196    gfits_free_header (&theader);
    199197    gfits_free_table (&ftable);
    … …  
    204202    gfits_create_table_header (&theader, "BINTABLE", "CELL_DATA");
    205203
    206     gfits_define_bintable_column (&theader, "D", "RA",          "ra (J2000) of cell center", "degree", 1.0, 0.0);
    207     gfits_define_bintable_column (&theader, "D", "DEC",         "dec (J2000) of cell center", "degree", 1.0, 0.0);
    208     gfits_define_bintable_column (&theader, "J", "ZONE",       "zone sequence number",     NULL,    1.0, 0.0);
    209     gfits_define_bintable_column (&theader, "J", "BAND",       "band sequence number",  NULL,    1.0, 0.0);
    210     gfits_define_bintable_column (&theader, "J", "INDEX",      "cell index",  NULL,    1.0, 0.0);
    211     gfits_define_bintable_column (&theader, "J", "NAME",       "cell name",                  NULL,    1.0, 0.0);
     204    char fmt[16];
     205    snprintf (fmt, 16, "%dA", BOUNDARY_TREE_NAME_LENGTH);
     206    gfits_define_bintable_column (&theader, "D", "RA",   "ra (J2000) of cell center", "degree", 1.0, 0.0);
     207    gfits_define_bintable_column (&theader, "D", "DEC",  "dec (J2000) of cell center", "degree", 1.0, 0.0);
     208    gfits_define_bintable_column (&theader, "J", "ZONE", "zone sequence number", "none", 1.0, 0.0);
     209    gfits_define_bintable_column (&theader, "J", "BAND", "band sequence number", "none", 1.0, 0.0);
     210    gfits_define_bintable_column (&theader, "J", "INDEX","cell index", "none", 1.0, 0.0);
     211    gfits_define_bintable_column (&theader, fmt, "NAME", "cell name", "none", 1.0, 0.0);
    212212
    213213    // generate the output array that carries the data
    … …  
    216216    int Ncell = 0;
    217217    for (i = 0; i < tree->Nzone; i++) {
    218         Ncell += tree->Nband[i];
     218      Ncell += tree->Nband[i];
    219219    }
    220220
    … …  
    226226    int    *band          ; ALLOCATE (band,  int,    Ncell);
    227227    int    *index         ; ALLOCATE (index, int,    Ncell);
    228     char  **name          ; ALLOCATE (name,  char *, Ncell);
     228    char   *name          ; ALLOCATE (name,  char,   Ncell*BOUNDARY_TREE_NAME_LENGTH);
     229
     230    // NOTE: a table column of characters must be fixed width, and is passed as a
     231    // contiguous array of Nchar * Nrow values
    229232
    230233    // assign the storage arrays
    231234    i = 0;
    232235    for (nz = 0; nz < tree->Nzone; nz++) {
    233     for (nb = 0; nb < tree->Nband[nz]; nb++) {
    234       R[i]     = tree->ra[nz][nb];
    235       D[i]     = tree->dec[nz][nb];
    236       zone[i]  = nz;
    237       band[i]  = nb;
    238       index[i] = i; // or tree->cells[nz][nb] ?
    239       name[i]  = tree->name[nz][nb]; // need to memcopy this?
    240       i++;
     236      for (nb = 0; nb < tree->Nband[nz]; nb++) {
     237        R[i]     = tree->ra[nz][nb];
     238        D[i]     = tree->dec[nz][nb];
     239        zone[i]  = nz;
     240        band[i]  = nb;
     241        index[i] = i; // or tree->cells[nz][nb] ?
     242        memcpy(&name[i*BOUNDARY_TREE_NAME_LENGTH], tree->name[nz][nb], BOUNDARY_TREE_NAME_LENGTH);
     243        i++;
     244      }
    241245    }
    242246
    243247    // add the columns to the output array
    244     gfits_set_bintable_column (&theader, &ftable, "RA",          R,     Ncell);
    245     gfits_set_bintable_column (&theader, &ftable, "DEC",         D,     Ncell);
    246     gfits_set_bintable_column (&theader, &ftable, "ZONE",        zone,  Ncell);
    247     gfits_set_bintable_column (&theader, &ftable, "BAND",        band,  Ncell);
    248     gfits_set_bintable_column (&theader, &ftable, "INDEX",       index, Ncell);
    249     gfits_set_bintable_column (&theader, &ftable, "NAME",        name,  Ncell);
     248    gfits_set_bintable_column (&theader, &ftable, "RA",    R,     Ncell);
     249    gfits_set_bintable_column (&theader, &ftable, "DEC",   D,     Ncell);
     250    gfits_set_bintable_column (&theader, &ftable, "ZONE",  zone,  Ncell);
     251    gfits_set_bintable_column (&theader, &ftable, "BAND",  band,  Ncell);
     252    gfits_set_bintable_column (&theader, &ftable, "INDEX", index, Ncell);
     253    gfits_set_bintable_column (&theader, &ftable, "NAME",  name,  Ncell);
    250254
    251255    free (R     );
    … …  
    263267  return TRUE;
    264268}
     269
     270// the boundary tree...
     271// given an (ra,dec) pair, find the containing projection cell.
     272
     273int BoundaryTreeCellCoords (BoundaryTree *tree, int *zone, int *band, double ra, double dec) {
     274
     275  // first, find the containing zone
     276
     277  // if we know dDEC, we can get the bin instantly:
     278  *zone = (dec - tree->DEC_origin) / tree->DEC_offset;
     279 
     280  if (*zone < 0) return FALSE;
     281  if (*zone >= tree->Nzone) return FALSE;
     282
     283  // now select the RA bin for that zone
     284  *band = (ra - tree->RA_origin[*zone]) / tree->RA_offset[*zone];
     285 
     286  if (*band < 0) return FALSE;
     287  if (*band >= tree->Nband[*zone]) *band = 0;
     288
     289  return TRUE;
     290}
     291

  • branches/eam_branches/ipp-20120627/Ohana/src/addstar/src/args_skycells.c

    r31239 r34210  
    3737    if (!strcasecmp (argv[N], "rings")) {
    3838      MODE = RINGS;
     39    }
     40    if (!strcasecmp (argv[N], "tamas")) {
     41      MODE = TAMAS;
    3942    }
    4043    remove_argument (N, &argc, argv);
    … …  
    179182    } 
    180183    remove_argument (N, &argc, argv);
     184  }
     185  if (MODE == TAMAS) {
     186    CELLSIZE = 3.955;
     187    if ((N = get_argument (argc, argv, "-cellsize"))) {
     188      remove_argument (N, &argc, argv);
     189      CELLSIZE = strtod (argv[N], &ptr);
     190      if ((*ptr != 0) || (CELLSIZE < 0.0)) {
     191        fprintf (stderr, "-cellsize requires a floating-point argument\n");
     192        help ();
     193      } 
     194      remove_argument (N, &argc, argv);
     195    }
    181196  }
    182197

  • branches/eam_branches/ipp-20120627/Ohana/src/addstar/src/findskycell.c

    r34209 r34210  
    1 # include "skycells.h"
     1# include "addstar.h"
    22
    33static double RA_offset_RINGS_V3[] = {360.000000,   40.000000,  24.000000,  17.142857,  13.333333,  10.909091,   9.230769,   8.000000,   7.200000,   
    … …  
    1212// in an even more specific case, RA[i,zone] = RA_origin[zone] + RA_offset[zone]
    1313
    14 typedef struct {
    15   int FixedGridDEC;           // is the DEC sequence linear?
    16   int FixedGridRA;            // in the RA sequence in a zone linear?
    17 
    18   double DEC_origin;
    19   double DEC_offset;
    20 
    21   int Nzone;
    22   double *RA_origin;
    23   double *RA_offset;
    24 
    25   int *Nband;
    26   int *NBAND;
    27 
    28   int **cells;
    29 } BoundaryTree;
    30 
    31 enum {NONE, MAKE_TREE, USE_TREE};
     14enum {TREE_NONE, TREE_MAKE, TREE_USE};
    3215
    3316void usage (void) {
    3417  fprintf (stderr, "USAGE: findcell -mktree (tree) (catdir)\n");
    35   fprintf (stderr, "USAGE: findcell -tree (tree) ra dec\n");
     18  fprintf (stderr, "USAGE: findcell -tree (tree) (datafile)\n");
     19  fprintf (stderr, "   (datafile) should contain a list of RA,DEC pairs\n");
    3620  exit (2);
    3721}
    3822
    3923int mktree (char *treefile, char *catdir);
    40 int find_primary_cell_coords (BoundaryTree *tree, int *zone, int *band, double ra, double dec);
     24int apply_tree (char *treefile, char *datafile);
    4125
    4226int main (int argc, char **argv) {
    … …  
    5539
    5640  /* extra error messages */
    57   MODE = NONE;
     41  int MODE = TREE_NONE;
    5842  if ((N = get_argument (argc, argv, "-mktree"))) {
    59     MODE = MAKE_TREE;
     43    MODE = TREE_MAKE;
    6044    remove_argument (N, &argc, argv);
    6145    treefile = strcreate (argv[N]);
    … …  
    6347  }
    6448  if ((N = get_argument (argc, argv, "-tree"))) {
    65     MODE = USE_TREE;
     49    MODE = TREE_USE;
    6650    remove_argument (N, &argc, argv);
    6751    treefile = strcreate (argv[N]);
    … …  
    7357  if (!MODE) usage();
    7458
    75   if (MODE == MAKE_TREE) {
     59  if (MODE == TREE_MAKE) {
    7660    mktree (treefile, argv[1]);
    7761    exit (0);
    7862  }
    7963
    80   // mktree (treefile, argv[1]);
     64  apply_tree (treefile, argv[1]);
    8165  exit (0);
    8266}
    … …  
    131115  ALLOCATE (tree.RA_origin, double, tree.Nzone);
    132116  ALLOCATE (tree.RA_offset, double, tree.Nzone);
    133   ALLOCATE (tree.cells,      int *, tree.Nzone);
     117
     118  ALLOCATE (tree.ra,   double *, tree.Nzone);
     119  ALLOCATE (tree.dec,  double *, tree.Nzone);
     120  ALLOCATE (tree.cell,    int *, tree.Nzone);
     121  ALLOCATE (tree.name,  char **, tree.Nzone);
     122
     123  // NOTE 1: the RA_origin, RA_offsets for RINGS.V3 are defined so that the first projection cell
     124  // overlaps the RA = 0,360 boundary.  This make the split a bit of a hack.  We end up
     125  // with Nbands, but the max boundary of the last band only goes to 360.0 -
     126  // 0.5*RA_offset.  To get the right band number for the boundary region, if the
     127  // calculation for the band number lands beyond the Nbands (ie, band >= Nbands), then
     128  // we need to loop back to the first band.
     129
     130  // NOTE 2: when we generate the zone & bands initially, we do not know the number of
     131  // bands in the end.  we cannot use the test of band >= Nband unless we set an absurdly
     132  // large default value.  Thus the value of 1000000 below.
    134133
    135134  // assign the bands for RINGS.V3
    136135  for (zone = 0; zone < tree.Nzone; zone++) {
    137     tree.Nband[zone] =  0;
     136    tree.Nband[zone] = 1000000;
    138137    tree.NBAND[zone] = 10;
    139138    tree.RA_origin[zone] = -0.5*RA_offset_RINGS_V3[zone];
    140139    tree.RA_offset[zone] = RA_offset_RINGS_V3[zone];
    141     ALLOCATE (tree.cells[zone], int, tree.NBAND[zone]);
     140    ALLOCATE (tree.ra[zone],   double, tree.NBAND[zone]);
     141    ALLOCATE (tree.dec[zone],  double, tree.NBAND[zone]);
     142    ALLOCATE (tree.cell[zone], int,    tree.NBAND[zone]);
     143    ALLOCATE (tree.name[zone], char *, tree.NBAND[zone]);
    142144    for (band = 0; band < tree.NBAND[zone]; band++) {
    143       tree.cells[zone][band] = -1;
     145      tree.ra[zone][band] = NAN;
     146      tree.dec[zone][band] = NAN;
     147      tree.cell[zone][band] = -1;
     148      ALLOCATE (tree.name[zone][band], char, BOUNDARY_TREE_NAME_LENGTH);
    144149    }
    145150  }
    … …  
    151156    XY_to_RD (&ra, &dec, x, y, &image[i].coords);
    152157
    153     if (!find_primary_cell_coords (&tree, &zone, &band, ra, dec)) {
     158    if (!BoundaryTreeCellCoords (&tree, &zone, &band, ra, dec)) {
    154159      fprintf (stderr, "mismatch!\n");
    155160      continue;
    … …  
    159164    if (band >= tree.NBAND[zone]) {
    160165      int start = tree.NBAND[zone];
    161       tree.Nband[zone] = band + 1;
    162166      tree.NBAND[zone] = band + 10;
    163       REALLOCATE (tree.cells[zone], int, tree.NBAND[zone]);
     167      REALLOCATE (tree.ra[zone],   double, tree.NBAND[zone]);
     168      REALLOCATE (tree.dec[zone],  double, tree.NBAND[zone]);
     169      REALLOCATE (tree.cell[zone], int,    tree.NBAND[zone]);
     170      REALLOCATE (tree.name[zone], char *, tree.NBAND[zone]);
    164171      for (j = start; j < tree.NBAND[zone]; j++) {
    165         tree.cells[zone][j] = -1;
     172        tree.ra[zone][j] = NAN;
     173        tree.dec[zone][j] = NAN;
     174        tree.cell[zone][j] = -1;
     175        ALLOCATE (tree.name[zone][j], char, BOUNDARY_TREE_NAME_LENGTH);
    166176      }
    167177    }
    168     tree.cells[zone][band] = i;
     178    tree.ra[zone][band] = ra;
     179    tree.dec[zone][band] = dec;
     180    tree.cell[zone][band] = i;
     181    memcpy (tree.name[zone][band], image[i].name, BOUNDARY_TREE_NAME_LENGTH);
    169182  }
    170183
    … …  
    175188    for (band = 0; band < tree.NBAND[zone]; band++) {
    176189      // all cells should be filled
    177       if (tree.cells[zone][band] < 0) {
     190      if (tree.cell[zone][band] < 0) {
    178191        if (!found_last) {
    179192          found_last = TRUE;
    … …  
    205218    ra  = 360.0 * drand48();
    206219    dec = 180.0 * drand48() - 90.0;
    207     if (!find_primary_cell_coords (&tree, &zone, &band, ra, dec)) {
     220    if (!BoundaryTreeCellCoords (&tree, &zone, &band, ra, dec)) {
    208221      fprintf (stderr, "failure for %f,%f\n", ra, dec);
    209222    }
    … …  
    216229}
    217230
    218 // the boundary tree...
    219 // given an (ra,dec) pair, find the containing projection cell.
    220 
    221 int find_primary_cell_coords (BoundaryTree *tree, int *zone, int *band, double ra, double dec) {
    222 
    223   // first, find the containing zone
    224 
    225   // if we know dDEC, we can get the bin instantly:
    226   *zone = (dec - tree->DEC_origin) / tree->DEC_offset;
    227  
    228   if (*zone < 0) return FALSE;
    229   if (*zone >= tree->Nzone) return FALSE;
    230 
    231   // now select the RA bin for that zone
    232   *band = (ra - tree->RA_origin[*zone]) / tree->RA_offset[*zone];
    233  
    234   if (*band < 0) return FALSE;
    235   // if (*band >= tree->Nband[*zone]) return FALSE;
    236 
    237   return TRUE;
    238 }
     231int apply_tree (char *treefile, char *datafile) {
     232
     233  BoundaryTree *tree = BoundaryTreeLoad (treefile);
     234  if (!tree) {
     235    fprintf (stderr, "error loading boundary tree file %s\n", treefile);
     236    exit (2);
     237  }
     238
     239  FILE *f = fopen (datafile, "r");
     240  if (!f) {
     241    fprintf (stderr, "error opening data file %s\n", datafile);
     242    exit (3);
     243  }
     244
     245  double ra, dec;
     246  int Nvalue = 0;
     247  while ((Nvalue = fscanf (f, "%lf %lf", &ra, &dec)) != EOF) {
     248
     249    int zone, band;
     250    if (!BoundaryTreeCellCoords (tree, &zone, &band, ra, dec)) {
     251      fprintf (stderr, "error finding cell for %f,%f\n", ra, dec);
     252      continue;
     253    }
     254
     255    fprintf (stdout, "%10.6f %10.6f  %3d %3d  %s\n", ra, dec, zone, band, tree->name[zone][band]);
     256  }
     257
     258  exit (0);
     259}

  • branches/eam_branches/ipp-20120627/Ohana/src/addstar/src/sky_tessalation.c

    r34088 r34210  
    2323      sky_tessellation_rings (db, level, Nmax);
    2424      return TRUE;
     25    case TAMAS:
     26      sky_tessellation_tamas (db, level, Nmax);
     27      return TRUE;
    2528    default:
    2629      break;
    … …  
    284287    free (ring);
    285288    free (image);
     289  }   
     290  return (TRUE);
     291}
     292
     293// the RINGS tessellation uses the declination zones proposed by Tamas Budavari,
     294// based on code supplied by Tamas 2012.07.23
     295int sky_tessellation_tamas (FITS_DB *db, int level, int Nmax) {
     296
     297  int j, nDEC, Nimage, Nring, Ntotal, Ndigit;
     298  double dec, dDEC;
     299  SkyRectangle *ring;
     300  Image *image;
     301  char format[16];
     302
     303  // The tessellation has one input parameter: the approximate cell size.  Starting with
     304  // the cell size, determine the optimal projection cell height (dDEC) that results in an
     305  // integer number of dec zones between -90 and +90
     306
     307  // in fact, we place a single image on each pole, so the real range of dec is 180.0 - CELLSIZE:
     308
     309  nDEC = (180.0 - CELLSIZE) / CELLSIZE;
     310  dDEC = (180.0 - CELLSIZE) / nDEC;
     311  nDEC += 2;
     312
     313  // how many total projection cells for this realization?  divide sky area by cell area:
     314  // this is used to set the number of digits, so it does not need to be very accurate...
     315  Ntotal = 41254.2 / (dDEC*dDEC);
     316  Ndigit = (int)(log10(Ntotal)) + 1 ;
     317  snprintf (format, 16, "skycell.%%0%dd", Ndigit);
     318
     319  double d2r = M_PI / 180; // is RAD_DEG
     320
     321  // parameter 'a' is the cell size in degrees
     322  double adeg = 3.955;
     323
     324  // half of 'a' in radians and its atan
     325  double halfa = adeg / 2 * d2r;
     326  double halftheta = atan(halfa);
     327
     328  // loop init
     329  dec = 0; // starting Decl. - could change this...
     330 
     331  while (dec < M_PI / 2 - halftheta) {
     332        double dm = dec - halftheta; // eq.5
     333        if (dec == 0) dm = 0; // initial
     334
     335        // dec is modified by the call below
     336        ring = sky_rectangle_tamas (&dec, dm, halfa, halftheta, &Nring, format);
     337        if (!ring) continue;
     338
     339        // subdivide each image (Nx x Ny subcells)
     340        Nimage = NX_SUB*NY_SUB*Nring;
     341        ALLOCATE (image, Image, Nimage);
     342        for (j = 0; j < Nring; j++) {
     343          // convert the SkyRectangles to Images for output
     344          sky_subdivide_image (&image[j*NX_SUB*NY_SUB], &ring[j], NX_SUB, NY_SUB);
     345          // printf("%s %8.2f %8.2f\n", ring[j].name, ring[j].coords.crval1, ring[j].coords.crval2);
     346        }
     347
     348        /* add the new images and save */
     349        dvo_image_addrows (db, image, Nimage);
     350        SetProtect (TRUE);
     351        dvo_image_update (db, VERBOSE);
     352        SetProtect (FALSE);
     353        dvo_image_clear_vtable (db);
     354   
     355        free (ring);
     356        free (image);
    286357  }   
    287358  return (TRUE);
    … …  
    666737}
    667738
     739// define the parameters of a projection centers for this ring
     740// dec : ~ center of ring in Dec
     741// dDEC : approximate height
     742// nring : number of cells generated for this ring
     743// format : guide to generate the filenames (c-type string format)
     744SkyRectangle *sky_rectangle_tamas (double *Dec, double dm, double halfa, double halftheta, int *nring, char *format) {
     745
     746  static int Nname = 0;
     747  int i, j, NX, NY;
     748  SkyRectangle *ring;
     749
     750  double d2r = M_PI / 180; // is RAD_DEG
     751  double dec = *Dec;
     752
     753  int nRA = (int)ceil(M_PI * cos(dm) / halftheta);  // eq.6       
     754  double dRA = 2 * M_PI / nRA; // eq.7
     755  double dp = atan(tan(dec + halftheta) * cos(dRA / 2)); // eq.9
     756
     757  if (dec == 0.0) {
     758    ALLOCATE (ring, SkyRectangle, nRA);
     759  } else {
     760    ALLOCATE (ring, SkyRectangle, 2*nRA);
     761  }
     762
     763  for (i = 0; i < nRA; i++) {
     764    // R.A. can use different phase per ring
     765    double ra = i * dRA; // + phase (watch wraparound)
     766
     767    int npass = (dec == 0.0) ? 1 : 2;
     768    for (j = 0; j < npass; j++) {
     769
     770      int N = j*nRA + i;
     771
     772      memset (&ring[N], 0, sizeof(SkyRectangle));
     773      memset (&ring[N].coords, 0, sizeof(Coords));
     774
     775      ring[N].coords.crval1 = ra / d2r;
     776      ring[N].coords.crval2 = (j == 0) ? dec / d2r : -dec / d2r;
     777
     778      printf(" \t %d   %25.20f   %25.20f\n", i, ring[N].coords.crval2, ring[N].coords.crval1);
     779
     780      ring[N].coords.pc1_1 = +1.0 * X_PARITY;
     781      ring[N].coords.pc1_2 = +0.0;
     782      ring[N].coords.pc2_1 = -0.0;
     783      ring[N].coords.pc2_2 = +1.0;
     784 
     785      // range values are in projected degrees
     786      NX = cos(dec - halftheta) * dRA   * 3600.0 / SCALE / d2r;
     787      NY =    2 * halftheta * 3600.0 / SCALE / d2r;
     788
     789      // crpix1,crpix2 is the projection center
     790      ring[N].coords.crpix1 = 0.5*NX;
     791      ring[N].coords.crpix2 = 0.5*NY;
     792
     793      ring[N].coords.cdelt1 = SCALE / 3600.0;
     794      ring[N].coords.cdelt2 = SCALE / 3600.0;
     795
     796      strcpy (ring[N].coords.ctype, "DEC--TAN");
     797
     798      ring[N].NX = NX*(1.0 + PADDING);
     799      ring[N].NY = NY*(1.0 + PADDING);
     800      ring[N].photcode = 1; // this needs to be set more sensibly
     801
     802      snprintf (ring[N].name, DVO_IMAGE_NAME_LEN, format, Nname);
     803      Nname++;
     804    }
     805  }
     806
     807  // advance to next ring
     808  *Dec = halftheta + dp;
     809
     810  *nring = (dec == 0.0) ? nRA : 2*nRA;
     811  return ring;
     812}
     813
    668814// an allocated image set is supplied, we fill in the values
    669815int sky_subdivide_image (Image *output, SkyRectangle *input, int Nx, int Ny) {
    … …  
    682828  }
    683829
    684   Ndigit = (int)(log10(Nx*Ny)) + 1 ;
    685   snprintf (format, 24, "%s.%%0%dd", input[0].name, Ndigit);
     830  if (Nx * Ny > 1) {
     831    Ndigit = (int)(log10(Nx*Ny)) + 1 ;
     832    snprintf (format, 24, "%s.%%0%dd", input[0].name, Ndigit);
     833  } else {
     834    snprintf (format, 24, "%s", input[0].name);
     835  }
    686836
    687837  // if requested extend, the skycell boundaries so that skycells overlap
    … …  
    696846      memcpy (&output[N].coords, &input[0].coords, sizeof(Coords));
    697847
    698       snprintf (output[N].name, DVO_IMAGE_NAME_LEN, format, N);
     848      if (Nx + Ny > 1) {
     849        snprintf (output[N].name, DVO_IMAGE_NAME_LEN, format, N);
     850      } else {
     851        snprintf (output[N].name, DVO_IMAGE_NAME_LEN, "%s", format);
     852      }
     853
    699854      output[N].NX = NX + 2 * pad_x;
    700855      output[N].NY = NY + 2 * pad_y;

  • branches/eam_branches/ipp-20120627/Ohana/src/libdvo/Makefile

    r34138 r34210  
    102102$(SRC)/db_utils.$(ARCH).o               \
    103103$(SRC)/convert.$(ARCH).o                \
    104 $(SRC)/HostTable.$(ARCH).o
     104$(SRC)/HostTable.$(ARCH).o              \
     105$(SRC)/BoundaryTree.$(ARCH).o
    105106
    106107

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

    r34207 r34210  
    293293} DVOTinyValueMode;
    294294
     295# define BOUNDARY_TREE_NAME_LENGTH 128
     296
     297typedef struct {
     298  int FixedGridDEC;           // is the DEC sequence linear?
     299  int FixedGridRA;            // in the RA sequence in a zone linear?
     300
     301  double DEC_origin;
     302  double DEC_offset;
     303
     304  int Nzone;
     305  double *RA_origin;
     306  double *RA_offset;
     307
     308  int *Nband;
     309  int *NBAND;
     310
     311  double   **ra;
     312  double  **dec;
     313  int    **cell;
     314  char  ***name;
     315} BoundaryTree;
     316
    295317// XXX DROP? // a reduced-subset structure for relastro
    296318// XXX DROP? typedef struct {
    … …  
    662684int free_tiny_values (Catalog *catalog);
    663685
     686int BoundaryTreeCellCoords (BoundaryTree *tree, int *zone, int *band, double ra, double dec);
     687int BoundaryTreeSave(char *filename, BoundaryTree *tree);
     688BoundaryTree *BoundaryTreeLoad(char *filename);
     689
    664690# endif // DVO_H
    665 

  • branches/eam_branches/ipp-20120627/Ohana/src/relphot/Makefile

    r33963 r34210  
    4949$(SRC)/setExclusions.$(ARCH).o   \
    5050$(SRC)/setMrelFinal.$(ARCH).o    \
     51$(SRC)/BoundaryTreeOps.$(ARCH).o         \
    5152$(SRC)/write_coords.$(ARCH).o
    5253
    … …  
    7778$(SRC)/setExclusions.$(ARCH).o   \
    7879$(SRC)/setMrelFinal.$(ARCH).o    \
     80$(SRC)/BoundaryTreeOps.$(ARCH).o         \
    7981$(SRC)/write_coords.$(ARCH).o
    8082

  • branches/eam_branches/ipp-20120627/Ohana/src/relphot/include/relphot.h

    r34207 r34210  
    111111char        *BCATALOG;
    112112ModeType     MODE;
     113
     114char        *BOUNDARY_TREE;
    113115
    114116double MAG_LIM;
    … …  
    341343int client_logger_message (char *format,...);
    342344
     345int MatchImageName (off_t meas, int cat, char *name);
     346
     347int load_tree (char *treefile);
     348int BoundaryTreePrimaryCell (char *primaryCellName, double ra, double dec);

  • branches/eam_branches/ipp-20120627/Ohana/src/relphot/src/ImageOps.c

    r34207 r34210  
    385385  distance = hypot (measure[0].Xccd - Xcenter, measure[0].Yccd - Ycenter);
    386386  return (distance);
     387}
     388
     389// returns image.Mcal - ff(x,y)
     390int MatchImageName (off_t meas, int cat, char *name) {
     391
     392  off_t i;
     393
     394  if (!name) return FALSE;
     395  if (!name[0]) return FALSE;
     396
     397  if (!MeasureToImage) return FALSE;
     398
     399  i = MeasureToImage[cat][meas];
     400  if (i == -1) return FALSE;
     401
     402  // this is a bit crude: stack image names are of the form
     403  // RINGS.V3.skycell.1495.027.sky.191211.stk.988232.cmf
     404  // the primaryCell has a name of the form RINGS.V3.skycell.1495
     405
     406  if (!strncmp(image[i].name, name, strlen(name))) return TRUE;
     407  return FALSE;
    387408}
    388409

  • branches/eam_branches/ipp-20120627/Ohana/src/relphot/src/StarOps.c

    r34207 r34210  
    306306int setMrel_catalog (Catalog *catalog, int Nc, int pass, FlatCorrectionTable *flatcorr, SetMrelInfo *results, int Nsecfilt) {
    307307
    308   off_t j, k, m;
     308  off_t j, k, m, ID;
    309309  int N;
    310310  float Msys, Mcal, Mmos, Mgrid;
    … …  
    326326  int isSetMrelFinal = (pass >= 0);
    327327
     328  char *primaryCell = NULL;
     329  ALLOCATE (primaryCell, char, DVO_MAX_PATH);
     330
    328331  for (j = 0; j < catalog[Nc].Naverage; j++) {
    329332    // XXX accumulate all secfilt values in a single pass?
    … …  
    334337      print_measure_set (&catalog[Nc].average[j], &catalog[Nc].secfilt[j*Nsecfilt], catalog[Nc].measure);
    335338    }
     339
     340    BoundaryTreePrimaryCell(primaryCell, catalog[Nc].average[j].R, catalog[Nc].average[j].D);
    336341
    337342    int GoodPS1 = FALSE;
    … …  
    361366      int haveStack = FALSE;
    362367
    363       float stackFluxPSF = NAN;
    364       float stackdFluxPSF = NAN;
    365       float stackFluxKron = NAN;
    366       float stackdFluxKron = NAN;
     368      // need to find the measurement closest to the center of its skycell, as well as the
     369      // closest for the subset of primary projection cells
     370
    367371      float stackCenterOffsetMin = 1e9;
    368       unsigned int stackImageIDmin;
     372      int stackCenterIDmin = -1;
     373      off_t stackCenterMeasureMin = -1;
     374
     375      float stackPrimaryOffsetMin = 1e9;
     376      int stackPrimaryIDmin = -1;
     377      off_t stackPrimaryMeasureMin = -1;
    369378
    370379      int forceSynth = FALSE;
    … …  
    436445
    437446            unsigned int stackImageID;
     447
     448            // which stack image should we use for the mean value?
     449            // if we request the primary (USE_TREE_FOR_PRIMARY), then find the min distances for data from the primary cell
     450            if (MatchImageName (m, Nc, primaryCell)) {
     451              float stackPrimaryOffset = getCenterOffset (m, Nc, &catalog[Nc].measure[m], &stackImageID);
     452              if (stackPrimaryOffset < stackPrimaryOffsetMin) {
     453                stackPrimaryIDmin = stackImageID;
     454                stackPrimaryMeasureMin = m;
     455              }
     456            }
     457
     458            // get the center distance for the generic case:
    438459            float stackCenterOffset = getCenterOffset (m, Nc, &catalog[Nc].measure[m], &stackImageID);
    439460            if (stackCenterOffset < stackCenterOffsetMin) {
    440               stackImageIDmin = stackImageID;
    441 
    442               float zp = Mcal + Mmos + Mgrid + PhotZeroPoint (&catalog[Nc].measure[m], &catalog[Nc].average[j], &catalog[Nc].secfilt[j*Nsecfilt]);
    443 
    444               // flux_cgs : erg sec^1 cm^-2 Hz^-1
    445               // mag_inst : -2.5 log (cts/sec)
    446               // mag_inst : -2.5 log (flux_inst)
    447               // flux_inst = ten(-0.4*mag_inst)
    448 
    449               // mag_AB = -2.5 log (flux_cgs) - 48.6 (~by definition) [~Vega flux in V-band]
    450               // flux_cgs = ten(-0.4*(mag_AB + 48.6))
    451 
    452               // flux_AB : ten(-0.4*mag_AB)
    453 
    454               // flux_cgs = ten(-0.4*48.6) * flux_AB
    455               // flux_AB  = ten(+0.4*48.6) * flux_cgs
    456            
    457               // flux_Jy : flux_cgs * 10^23
    458 
    459               // flux_AB = ten(+0.4*48.6) * ten(-23) * flux_Jy
    460 
    461               // mag_AB = mag_inst + ZP
    462 
    463               // flux_inst = ten(-0.4*(mag_AB - ZP)) = ten(0.4*ZP) * flux_AB
    464 
    465               // flux_AB = flux_inst * ten(-0.4*ZP)
    466 
    467               // flux_inst * ten(-0.4*ZP) = ten(+0.4*48.6 - 23) * flux_Jy
    468 
    469               // flux_inst = flux_Jy * ten(0.4*ZP + 0.4*48.6 - 23)
    470               // flux_inst = flux_Jy * ten(0.4*ZP - 3.56)
    471               // flux_Jy = flux_inst * ten(-0.4*ZP + 3.56)
    472 
    473               // zpFactor to go from instrumental flux to Janskies
    474               float zpFactor = pow(10.0, -0.4*zp + 3.56);
    475 
    476               // need to put in AB mag factor to get to Janskies (or uJy?)
    477               stackFluxPSF   = zpFactor * catalog[Nc].measure[m].FluxPSF;
    478               stackdFluxPSF  = zpFactor * catalog[Nc].measure[m].dFluxPSF;
    479               stackFluxKron  = zpFactor * catalog[Nc].measure[m].FluxKron;
    480               stackdFluxKron = zpFactor * catalog[Nc].measure[m].dFluxKron;
     461              stackCenterIDmin = stackImageID;
     462              stackCenterMeasureMin = m;
    481463            }
    482464          }
    … …  
    619601
    620602        if (haveStack) {
    621           catalog[Nc].secfilt[Nsecfilt*j+Nsec].FluxPSF = stackFluxPSF;
    622           catalog[Nc].secfilt[Nsecfilt*j+Nsec].dFluxPSF = stackdFluxPSF;
    623           catalog[Nc].secfilt[Nsecfilt*j+Nsec].FluxKron = stackFluxKron;
    624           catalog[Nc].secfilt[Nsecfilt*j+Nsec].dFluxKron = stackdFluxKron;
    625           catalog[Nc].secfilt[Nsecfilt*j+Nsec].stackID = stackImageIDmin;
     603          m  = (stackPrimaryMeasureMin >= 0) ? stackPrimaryMeasureMin : stackCenterMeasureMin;
     604          ID = (stackPrimaryMeasureMin >= 0) ? stackPrimaryIDmin      : stackCenterIDmin;
     605
     606          // get the zero point for the selected image
     607          float zp = Mcal + Mmos + Mgrid + PhotZeroPoint (&catalog[Nc].measure[m], &catalog[Nc].average[j], &catalog[Nc].secfilt[j*Nsecfilt]);
     608
     609          // flux_cgs : erg sec^1 cm^-2 Hz^-1
     610          // mag_inst : -2.5 log (cts/sec)
     611          // mag_inst : -2.5 log (flux_inst)
     612          // flux_inst = ten(-0.4*mag_inst)
     613
     614          // mag_AB = -2.5 log (flux_cgs) - 48.6 (~by definition) [~Vega flux in V-band]
     615          // flux_cgs = ten(-0.4*(mag_AB + 48.6))
     616
     617          // flux_AB : ten(-0.4*mag_AB)
     618
     619          // flux_cgs = ten(-0.4*48.6) * flux_AB
     620          // flux_AB  = ten(+0.4*48.6) * flux_cgs
     621           
     622          // flux_Jy : flux_cgs * 10^23
     623
     624          // flux_AB = ten(+0.4*48.6) * ten(-23) * flux_Jy
     625
     626          // mag_AB = mag_inst + ZP
     627
     628          // flux_inst = ten(-0.4*(mag_AB - ZP)) = ten(0.4*ZP) * flux_AB
     629
     630          // flux_AB = flux_inst * ten(-0.4*ZP)
     631
     632          // flux_inst * ten(-0.4*ZP) = ten(+0.4*48.6 - 23) * flux_Jy
     633
     634          // flux_inst = flux_Jy * ten(0.4*ZP + 0.4*48.6 - 23)
     635          // flux_inst = flux_Jy * ten(0.4*ZP - 3.56)
     636          // flux_Jy = flux_inst * ten(-0.4*ZP + 3.56)
     637
     638          // zpFactor to go from instrumental flux to Janskies
     639          float zpFactor = pow(10.0, -0.4*zp + 3.56);
     640
     641          // need to put in AB mag factor to get to Janskies (or uJy?)
     642          catalog[Nc].secfilt[Nsecfilt*j+Nsec].FluxPSF   = zpFactor * catalog[Nc].measure[m].FluxPSF; 
     643          catalog[Nc].secfilt[Nsecfilt*j+Nsec].dFluxPSF  = zpFactor * catalog[Nc].measure[m].dFluxPSF;
     644          catalog[Nc].secfilt[Nsecfilt*j+Nsec].FluxKron  = zpFactor * catalog[Nc].measure[m].FluxKron;
     645          catalog[Nc].secfilt[Nsecfilt*j+Nsec].dFluxKron = zpFactor * catalog[Nc].measure[m].dFluxKron;
     646
     647          catalog[Nc].secfilt[Nsecfilt*j+Nsec].stackID   = ID;
    626648        }
    627649
    … …  
    680702    }
    681703  }
     704  if (primaryCell) free (primaryCell);
    682705  return (TRUE);
    683706}

  • branches/eam_branches/ipp-20120627/Ohana/src/relphot/src/args.c

    r33963 r34210  
    197197  }
    198198
     199  if ((N = get_argument (argc, argv, "-boundary-tree"))) {
     200    remove_argument (N, &argc, argv);
     201    load_tree (argv[N]);
     202    remove_argument (N, &argc, argv);
     203  }
     204
    199205  SHOW_PARAMS = FALSE;
    200206  if ((N = get_argument (argc, argv, "-params"))) {
Note: See TracChangeset for help on using the changeset viewer.