IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 38179


Ignore:
Timestamp:
Apr 23, 2015, 8:22:11 PM (11 years ago)
Author:
eugene
Message:

adding ukirt

Location:
branches/eam_branches/ipp-20150419/Ohana/src/addstar/src
Files:
2 edited

Legend:

Unmodified
Added
Removed

  • branches/eam_branches/ipp-20150419/Ohana/src/addstar/src/LoadDataUKIRT.c

    r38169 r38179  
    22
    33// examine the header sets and set the Image entries for the the valid images
    4 // there should only be a single data set (phu + table) in this file
    5 // each SDSS data set corresponds to 5 images (ugriz)
    64int LoadDataUKIRT (FILE *f, char *file, Image **images, off_t *nimages, Stars **stars, unsigned int *Nstars, Header **headers, off_t *extsize, HeaderSet *headerSets, off_t NheaderSets) {
    75
    … …  
    1513  if (images[0] == NULL) {
    1614    Nvalid = 0;
    17     NVALID = 4;
     15    NVALID = NheaderSets;
    1816    ALLOCATE (images[0], Image, NVALID);
    1917  } else {
    2018    Nvalid = *nimages;
    21     NVALID = Nvalid + 4;
     19    NVALID = Nvalid + NheaderSets;
    2220    REALLOCATE (images[0], Image, NVALID);
    2321  }   

  • branches/eam_branches/ipp-20150419/Ohana/src/addstar/src/ReadStarsUKIRT.c

    r38169 r38179  
    11# include "addstar.h"
    22
    3 int SetSDSSFlags (Stars *star, unsigned int flags1, unsigned int flags2);
    4 
    5 # define NFILTER 5
    6 
    7 # define GET_COLUMN_5(NAME,TYPE) \
    8   TYPE *NAME; \
    9   NAME = (TYPE *) gfits_get_bintable_column_data (table.header, &table, #NAME, type, &Nrow, &Ncol); \
    10   assert (NAME); assert (!strcmp (type, #TYPE)); assert (Nrow == Nstars); assert (Ncol == NFILTER);
    11 
    12 # define GET_COLUMN_5_NOASSERT(NAME,TYPE) \
    13   TYPE *NAME; \
    14   NAME = (TYPE *) gfits_get_bintable_column_data (table.header, &table, #NAME, type, &Nrow, &Ncol);
    15 
    16 # define GET_COLUMN_1(NAME,TYPE) \
     3# define GET_COLUMN(NAME,TYPE) \
    174  TYPE *NAME; \
    185  NAME = (TYPE *) gfits_get_bintable_column_data (table.header, &table, #NAME, type, &Nrow, &Ncol); \
    … …  
    4633  int Ncol; // used in the GET_COLUMN_1,5 macros above
    4734 
     35  // XXX ??? what is this for?
    4836  if (in_theader == NULL) {
    4937    table.header = &theader;
    … …  
    6149  }
    6250
    63   //
    64   int hasDetID = gfits_scan (header, "DETECTID", "%s", 1, detectID);
    65   if (!hasDetID) {
     51  Coords coords;
     52  if (!GetCoords (&coords, &table.header)) {
     53    fprintf (stderr, "unable to read header WCS info\n");
     54    exit (3);
     55  }
     56
     57  // DETECTID tells us which chip supplied the data
     58  char detIDstr[80];
     59  if (!gfits_scan (header, "DETECTID", "%s", 1, detIDstr)) {
    6660    fprintf (stderr, "missing DETECTID\n");
    6761    exit (3);
    6862  }
    6963 
    70  
    71 
    72 w20050611_00700_st_cat.fits[APM-BINARYTABLE]  RSC:H2:60 
    73 w20050611_00700_st_cat.fits[APM-BINARYTABLE]  RSC:H2:63 
    74 w20050611_00700_st_cat.fits[APM-BINARYTABLE]  RSC:H2:76 
    75 w20050611_00700_st_cat.fits[APM-BINARYTABLE]  RSC:H2:41 
    76 
    77   strcpy (filtname[0], "u");
    78   strcpy (filtname[1], "g");
    79   strcpy (filtname[2], "r");
    80   strcpy (filtname[3], "i");
    81   strcpy (filtname[4], "z");
    82 
    83   NAMED_PHOTCODE_AND_ZP (photcode[0], zeropt[0], "U_SDSS");
    84   NAMED_PHOTCODE_AND_ZP (photcode[1], zeropt[1], "G_SDSS");
    85   NAMED_PHOTCODE_AND_ZP (photcode[2], zeropt[2], "R_SDSS");
    86   NAMED_PHOTCODE_AND_ZP (photcode[3], zeropt[3], "I_SDSS");
    87   NAMED_PHOTCODE_AND_ZP (photcode[4], zeropt[4], "Z_SDSS");
    88 
    89   // XXXYYYZZZ SDSS tables have special flags for undefined/unmeasured values and errors:
    90   // -9999 flags unmeasured values, and the corresponding error may or may not be meaningful
    91   // -1000 flags errors that are not determined, even though the corresponding quantity is.
    92   // These special values need to be trapped here to avoid averaging meaningful numbers with the flag values.
    93 
    94   // various header values needed to calculate per-star data below
    95   gfits_scan (header, "C_OBS", "%lf", 1, &clockRate); // value in header is usec / unbinned row
    96   clockRate *= 1e-6; // convert to seconds / unbinned row
    97 
    98   gfits_scan (table.header, "MJD_U", "%lf", 1, &mjd[0]);
    99   gfits_scan (table.header, "MJD_G", "%lf", 1, &mjd[1]);
    100   gfits_scan (table.header, "MJD_R", "%lf", 1, &mjd[2]);
    101   gfits_scan (table.header, "MJD_I", "%lf", 1, &mjd[3]);
    102   gfits_scan (table.header, "MJD_Z", "%lf", 1, &mjd[4]);
    103   tzero[0] = ohana_mjd_to_sec (mjd[0]);
    104   tzero[1] = ohana_mjd_to_sec (mjd[1]);
    105   tzero[2] = ohana_mjd_to_sec (mjd[2]);
    106   tzero[3] = ohana_mjd_to_sec (mjd[3]);
    107   tzero[4] = ohana_mjd_to_sec (mjd[4]);
    108 
    109   gfits_scan (table.header, "SEEING_U", "%f", 1, &seeing[0]);
    110   gfits_scan (table.header, "SEEING_G", "%f", 1, &seeing[1]);
    111   gfits_scan (table.header, "SEEING_R", "%f", 1, &seeing[2]);
    112   gfits_scan (table.header, "SEEING_I", "%f", 1, &seeing[3]);
    113   gfits_scan (table.header, "SEEING_Z", "%f", 1, &seeing[4]);
    114 
    115   gfits_scan (table.header, "PSFERR_U", "%f", 1, &photErr[0]);
    116   gfits_scan (table.header, "PSFERR_G", "%f", 1, &photErr[1]);
    117   gfits_scan (table.header, "PSFERR_R", "%f", 1, &photErr[2]);
    118   gfits_scan (table.header, "PSFERR_I", "%f", 1, &photErr[3]);
    119   gfits_scan (table.header, "PSFERR_Z", "%f", 1, &photErr[4]);
    120 
    121   gfits_scan (header, "CAMCOL", "%d", 1, &camcol); // value in header is usec / unbinned row
     64  int detID = 0;
     65  if (!strcmp (detIDstr, "RSC:H2:60")) { detID = 1; }
     66  if (!strcmp (detIDstr, "RSC:H2:63")) { detID = 2; }
     67  if (!strcmp (detIDstr, "RSC:H2:76")) { detID = 3; }
     68  if (!strcmp (detIDstr, "RSC:H2:41")) { detID = 4; }
     69  if (!detID) {
     70    fprintf (stderr, "unknown DETECTID %s\n", detIDstr);
     71    exit (3);
     72  }
     73
     74  char filter[80];
     75  if (!gfits_scan (header, "FILTER", "%s", 1, filter)) {
     76    fprintf (stderr, "missing FILTER\n");
     77    exit (3);
     78  }
     79
     80  int knownFilter = FALSE;
     81  if (!strcmp (filter, "Y")) { knownFilter = TRUE; }
     82  if (!strcmp (filter, "J")) { knownFilter = TRUE; }
     83  if (!strcmp (filter, "H")) { knownFilter = TRUE; }
     84  if (!strcmp (filter, "K")) { knownFilter = TRUE; }
     85  if (!knownFilter) {
     86    fprintf (stderr, "unknown FILTER %s\n", filter);
     87    exit (3);
     88  }
     89 
     90  sprintf (photcodeName, "UKIRT.%s.%d", filter, detID);
     91  PhotCode *code = GetPhotcodebyName (photcodeName);
     92  if (code == NULL) {
     93    fprintf (stderr, "ERROR:  photcode %s not found in photcode table\n", photcodeName);
     94    exit (3);
     95  }
     96  int photcode = code[0].code;
     97
     98  if (!gfits_scan (header, "MAGZPT", "%f", 1, &zeropt)) {
     99    fprintf (stderr, "ERROR:  missing MAGZPT\n");
     100    exit (3);
     101  }
     102
     103  gfits_scan (header, "MJD_OBS", "%lf", 1, &mjd);
     104  tzero = ohana_mjd_to_sec (mjd);
     105  double jd = ohana_sec_to_jd (tzero);
     106  double sidtime  = 15.0*ohana_lst (jd, Longitude); // sidtime in degrees
     107
     108  float seeing;
     109  gfits_scan (header, "SEEING", "%f", 1, &seeing);
     110 
     111
     112  float exptime;
     113  gfits_scan (header, "EXP_TIME", "%f", 1, &exptime);
     114  float magtime = 2.5*log10(exptime);
    122115
    123116  ZeroPt = GetZeroPoint();
    … …  
    127120  ALLOCATE (stars, Stars, NFILTER*Nstars);
    128121
    129   GET_COLUMN_5 (rowc, float);
    130   GET_COLUMN_5 (colc, float);
    131   GET_COLUMN_5 (sky, float);
    132   GET_COLUMN_5 (psfCounts, float);
    133   GET_COLUMN_5 (fiberCounts, float);
    134   GET_COLUMN_5 (offsetRa, float);
    135   GET_COLUMN_5 (offsetDec, float);
    136   GET_COLUMN_5 (flags, int);
    137   GET_COLUMN_5 (flags2, int);
    138 
    139 #ifdef notyet
    140   GET_COLUMN_5 (prob_psf, float);
    141 #else
    142   GET_COLUMN_5_NOASSERT (prob_psf, float);
    143 #endif
    144 
    145   GET_COLUMN_1 (ra, double);
    146   GET_COLUMN_1 (dec, double);
    147 
    148   GET_COLUMN_5 (rowcErr, float);
    149   GET_COLUMN_5 (colcErr, float);
    150   GET_COLUMN_5 (skyErr, float);
    151   GET_COLUMN_5 (psfCountsErr, float);
     122  GET_COLUMN (X_coordinate, float);
     123  GET_COLUMN (Y_coordinate, float);
     124
     125  GET_COLUMN (X_coordinate_err, float);
     126  GET_COLUMN (Y_coordinate_err, float);
     127
     128  GET_COLUMN (Aper_flux_8, float);
     129  GET_COLUMN (Aper_flux_8_err, float);
     130
     131  GET_COLUMN (Aper_flux_5, float);
     132  GET_COLUMN (Aper_flux_5_err, float);
     133
     134  GET_COLUMN (Kron_flux, float);
     135  GET_COLUMN (Kron_flux_err, float);
     136
     137  GET_COLUMN (Sky_level, float);
     138  GET_COLUMN (Sky_rms, float);
    152139
    153140  // the value of stars[].M is supposed to be the instrumental magnitude offset by the
    … …  
    157144
    158145  for (i = 0; i < Nstars; i++) {
    159     for (j = 0; j < NFILTER; j++) {
    160       N = NFILTER*i + j;
    161       InitStar (&stars[N]);
     146      InitStar (&stars[i]);
    162147     
    163148      // any values not explicitly set are left at 0.0
    164       stars[N].average.R         = ra[i] + dCOS(dec[i]) * offsetRa[N] / 3600.0;
    165       stars[N].average.D         = dec[i] + offsetDec[N] / 3600.0;
    166       stars[N].average.dR        = NAN;
    167       stars[N].average.dD        = NAN;
    168 
    169       stars[N].measure.Xccd      = colc[N];
    170       stars[N].measure.Yccd      = rowc[N];
    171       stars[N].measure.dXccd     = ToShortPixels(colcErr[N]);
    172       stars[N].measure.dYccd     = ToShortPixels(rowcErr[N]);
    173       stars[N].measure.M         = psfCounts[N] + ZeroPt - zeropt[j];
    174       stars[N].measure.dM        = psfCountsErr[N];
    175       stars[N].measure.Map       = fiberCounts[N] + ZeroPt - zeropt[j];
    176       stars[N].measure.Sky       = sky[N]; // adjust this to counts?
    177       stars[N].measure.dSky      = skyErr[N];
    178       stars[N].measure.FWx       = ToShortPixels(seeing[j]); // reported in arcsec?
    179       stars[N].measure.FWy       = ToShortPixels(seeing[j]);
    180       if (prob_psf) {
    181           stars[N].measure.psfChisq  = prob_psf[N]; // XXX not really the correct value...
    182       } else {
    183           stars[N].measure.psfChisq  = NAN;
    184       }
    185       stars[N].measure.detID     = N;
    186       stars[N].measure.t         = tzero[j] + clockRate*rowc[N]; // time since row 0
    187       stars[N].measure.dt        = 4.32912209; // 2.5 * log(53.907456) the sdss exposure time // old comment is 53907456 is this 2048*clockRate ?
    188 
    189       SetSDSSFlags (&stars[N], flags[N], flags2[N]);
    190 
    191       // longitude and latitude for SDSS:
    192       // Latitude 32° 46' 49.30" N, Longitude 105° 49' 13.50" W
    193       // longitude = 105.820419312 deg = 7.05469417
     149      XY_to_RD (&R, &D, X_coordinate[i], Y_coordinate[i], &coords);
     150
     151      stars[i].average.R         = R;
     152      stars[i].average.D         = D;
     153      stars[i].average.dR        = NAN;
     154      stars[i].average.dD        = NAN;
     155
     156      stars[i].measure.Xccd      = X_coordinate[i];
     157      stars[i].measure.Yccd      = Y_coordinate[i];
     158      stars[i].measure.dXccd     = ToShortPixels(X_coordinate_err[i]);
     159      stars[i].measure.dYccd     = ToShortPixels(Y_coordinate_err[i]);
     160      stars[i].measure.M         = -2.5*log10(Aper_flux_8[i]) + ZeroPt - zeropt;
     161      stars[i].measure.dM        = Aper_flux_8_err[i] / Aper_flux_8[i];
     162      stars[i].measure.Map       = -2.5*log10(Aper_flux_8[i]) + ZeroPt - zeropt;
     163      stars[i].measure.dMap      = Aper_flux_5_err[i] / Aper_flux_5[i];
     164      stars[i].measure.Mkron     = -2.5*log10(Aper_flux_8[i]) + ZeroPt - zeropt;
     165      stars[i].measure.dMkron    = Aper_flux_5_err[i] / Aper_flux_5[i];
     166      stars[i].measure.Sky       = Sky_level[i]; // adjust this to counts?
     167      stars[i].measure.dSky      = Sky_rms[i];
     168      stars[i].measure.FWx       = ToShortPixels(seeing); // reported in arcsec?
     169      stars[i].measure.FWy       = ToShortPixels(seeing);
     170
     171      stars[i].measure.detID     = i;
     172      stars[i].measure.t         = tzero; // time since row 0
     173      stars[i].measure.dt        = magtime; // 2.5 * log(exptime)
     174
     175      // longitude and latitude for UKIRT:
     176      // Latitude , Longitude ???
     177      // longitude = 105.820419312 deg = 7.05469417 h
    194178      // latitude = 32.7803611755 deg
    195179
    196       double Longitude = 7.05469417;   // hours (+ = W)
    197       double Latitude = 32.7803611755; // degrees
    198 
    199       jd = ohana_sec_to_jd (stars[N].measure.t);
    200       sidtime  = 15.0*ohana_lst (jd, Longitude); // sidtime in degrees
    201       altaz (&alt, &az, sidtime - stars[N].average.R, stars[N].average.D, Latitude);
    202 
    203       stars[N].measure.airmass   = 1.0 / dCOS(90.0 - alt);
    204       stars[N].measure.az        = az;
    205       stars[N].measure.photcode  = photcode[j];
    206       stars[N].measure.imageID   = j + *nimages; // set imageID to entry for this filter
     180      double Longitude = 7.05469417 * foo;   // hours (+ = W)
     181      double Latitude = 32.7803611755 * fooo; // degrees
     182
     183      altaz (&alt, &az, sidtime - stars[i].average.R, stars[i].average.D, Latitude);
     184
     185      stars[i].measure.airmass   = 1.0 / dCOS(90.0 - alt);
     186      stars[i].measure.az        = az;
     187      stars[i].measure.photcode  = photcode;
     188      stars[i].measure.imageID   = *nimages; // set imageID to entry for this filter
    207189    }
    208190  }   
    … …  
    278260  return (stars);
    279261}
    280 
    281 int SetSDSSFlags (Stars *star, unsigned int flags1, unsigned int flags2) {
    282 
    283   // XXX this is wrong, need to roll left to set the correct bit
    284   if (flags1 & 0x00000002) star[0].measure.photFlags |= 0x0001; // BRIGHT            - 1  1
    285   if (flags1 & 0x00000004) star[0].measure.photFlags |= 0x0002; // EDGE              - 1  2
    286   if (flags1 & 0x00000008) star[0].measure.photFlags |= 0x0004; // BLENDED           - 1  3
    287   if (flags1 & 0x00000010) star[0].measure.photFlags |= 0x0008; // CHILD             - 1  4
    288   if (flags1 & 0x00000020) star[0].measure.photFlags |= 0x0010; // PEAKCENTER        - 1  5
    289   if (flags1 & 0x00000040) star[0].measure.photFlags |= 0x0020; // NODEBLEND         - 1  6
    290   if (flags1 & 0x00040000) star[0].measure.photFlags |= 0x0040; // SATUR             - 1 18
    291   if (flags1 & 0x00080000) star[0].measure.photFlags |= 0x0080; // NOTCHECKED        - 1 19
    292   if (flags1 & 0x10000000) star[0].measure.photFlags |= 0x0100; // BINNED1           - 1 28
    293   if (flags1 & 0x20000000) star[0].measure.photFlags |= 0x0200; // BINNED2           - 1 29
    294   if (flags1 & 0x40000000) star[0].measure.photFlags |= 0x0400; // BINNED4           - 1 30
    295   if (flags2 & 0x00000040) star[0].measure.photFlags |= 0x0800; // LOCAL_EDGE        - 2  7
    296   if (flags2 & 0x00000800) star[0].measure.photFlags |= 0x1000; // INTERP_CENTER     - 2 12
    297   if (flags2 & 0x00002000) star[0].measure.photFlags |= 0x2000; // DEBLEND_NOPEAK    - 2 14
    298   if (flags2 & 0x02000000) star[0].measure.photFlags |= 0x4000; // NOTCHECKED_CENTER - 2 26
    299   return (TRUE);
    300 
    301 }
Note: See TracChangeset for help on using the changeset viewer.