IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 38184


Ignore:
Timestamp:
Apr 24, 2015, 3:47:00 PM (11 years ago)
Author:
eugene
Message:

add code to read UKIRT data

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

Legend:

Unmodified
Added
Removed

  • branches/eam_branches/ipp-20150419/Ohana/src/addstar/Makefile

    r37807 r38184  
    7676$(SRC)/LoadData.$(ARCH).o \
    7777$(SRC)/LoadDataSDSS.$(ARCH).o \
     78$(SRC)/LoadDataUKIRT.$(ARCH).o \
    7879$(SRC)/LoadDataPMM.$(ARCH).o \
    7980$(SRC)/in_image.$(ARCH).o \
    … …  
    9091$(SRC)/ReadStarsTEXT.$(ARCH).o \
    9192$(SRC)/ReadStarsSDSS.$(ARCH).o \
     93$(SRC)/ReadStarsUKIRT.$(ARCH).o \
    9294$(SRC)/ReadXradFITS.$(ARCH).o \
    9395$(SRC)/FilterStars.$(ARCH).o \
    … …  
    204206$(SRC)/LoadData.$(ARCH).o \
    205207$(SRC)/LoadDataSDSS.$(ARCH).o \
     208$(SRC)/LoadDataUKIRT.$(ARCH).o \
    206209$(SRC)/LoadDataPMM.$(ARCH).o \
    207210$(SRC)/grefstars.$(ARCH).o \
    … …  
    213216$(SRC)/ReadStarsTEXT.$(ARCH).o \
    214217$(SRC)/ReadStarsSDSS.$(ARCH).o \
     218$(SRC)/ReadStarsUKIRT.$(ARCH).o \
    215219$(SRC)/ImageIndex.$(ARCH).o \
    216220$(SRC)/UpdateImageIDs.$(ARCH).o \

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

    r38167 r38184  
    332332
    333333
     334Stars *ReadStarsUKIRT (FILE *f, char *name, Header *header, Image *images, off_t *nimages, unsigned int *nstars);
     335int 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);
     336
    334337/**
    335338    there is an inconsistency to be resolved: fixed structures (like Image)

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

    r38167 r38184  
    3131  int have_TELESCOPE  = gfits_scan (header, "TELESCOP", "%s", 1, telescope);
    3232  int have_INSTRUMENT = gfits_scan (header, "INSTRUME", "%s", 1, instrument);
    33   if (haveTELESCOPE && have_INSTRUMENT && !strcmp(telescope, "UKIRT"), && !strcmp(instrument, "WFCAM")) {
     33  if (have_TELESCOPE && have_INSTRUMENT && !strcmp(telescope, "UKIRT") && !strcmp(instrument, "WFCAM")) {
    3434    if (VERBOSE) fprintf (stderr, "found UKIRT data\n");
    3535    return UKIRT_OBJ;

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

    r38179 r38184  
    22
    33// examine the header sets and set the Image entries for the the valid images
     4// UKIRT data has the WCS/image metadata header intermixed with the bintable header
    45int 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) {
    56
    67  off_t Nskip, Nvalid, NVALID;
    78  char *name;
    8   int j, Nhead, Ndata;
     9  int i, j, Nhead, Ndata;
    910  unsigned int Ninstars;
    1011  Stars *inStars;
    … …  
    2930
    3031    // advance the pointer to the start of the corresponding table block
    31     Nhead = headerSets[0].extnum_head;
     32    Nhead = headerSets[i].extnum_head;
    3233    Ndata = headerSets[i].extnum_data;
    3334    Nskip = 0;
    … …  
    3738    fseeko (f, Nskip, SEEK_SET);
    3839
    39     inStars = ReadStarsUKIRT (f, name, headers[Nhead], headers[Ndata], images[0], &Nvalid, &Ninstars);
     40    inStars = ReadStarsUKIRT (f, name, headers[Nhead], images[0], &Nvalid, &Ninstars);
    4041    *stars = MergeStars (*stars, Nstars, inStars, Ninstars);
    4142
    42     free (name);
    4343    *nimages = Nvalid;
    4444  }   

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

    r38169 r38184  
    5353      }
    5454    }
     55
     56    // the 0th entry is an empty PHU
     57    extsize[0][0] = headers[0][0].datasize + gfits_data_size (headers[0]);
    5558
    5659    for (i = 1; i < Nheaders; i++) {

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

    r38179 r38184  
    66  assert (NAME); assert (!strcmp (type, #TYPE)); assert (Nrow == Nstars); assert (Ncol == 1);
    77
    8 /* grab named photcode */
    9 # define NAMED_PHOTCODE_AND_ZP(CODE,ZP,NAME) { \
    10   PhotCode *code; \
    11   code = GetPhotcodebyName (NAME); \
    12   if (code == NULL) { \
    13     fprintf (stderr, "ERROR:  photcode %s not found in photcode table\n", NAME); \
    14     exit (0); } \
    15   CODE = code[0].code; \
    16   ZP = 0.001*code[0].C; }
    17 
    188// given a file with the pointer at the start of the table block and the
    199// corresponding image header, load the stars from the table
    20 Stars *ReadStarsSDSS (FILE *f, char *name, Header *header, Header *in_theader, Image *images, off_t *nimages, unsigned int *nstars) {
    21 
    22   off_t Nskip, Nrow;
    23   int i, j, N, Nstars, camcol;
     10Stars *ReadStarsUKIRT (FILE *f, char *name, Header *header, Image *images, off_t *nimages, unsigned int *nstars) {
     11
     12  off_t Nrow;
     13  int Ncol;
    2414  char type[80];
    25   Header theader;
    2615  FTable table;
    2716  Stars *stars;
    28   double clockRate, mjd[5], jd, sidtime, alt, az;
    29   float seeing[5], photErr[5], zeropt[5], ZeroPt;
    30   time_t tzero[5];
    31   char filtname[16][5];
    32   int photcode[5];
    33   int Ncol; // used in the GET_COLUMN_1,5 macros above
    34  
    35   // XXX ??? what is this for?
    36   if (in_theader == NULL) {
    37     table.header = &theader;
    38     if (!gfits_fread_header (f, table.header)) Shutdown ("ERROR: can't read table header");
    39   } else {
    40     table.header = in_theader;
    41     Nskip = in_theader[0].datasize;
    42     fseeko (f, Nskip, SEEK_CUR);
    43   }
     17 
     18  // ohana_memcheck_func (TRUE);
     19
     20  // the FITS binary table header is the same as the image header
     21  table.header = header;
     22
     23  // put the file pointer at the start of the fits table data
     24  off_t Nskip = header[0].datasize;
     25  fseeko (f, Nskip, SEEK_CUR);
    4426
    4527  /* load the table data */
    4628  if (!gfits_fread_ftable_data (f, &table, FALSE)) {
    47     fprintf (stderr, "ERROR: can't read table header\n");
     29    fprintf (stderr, "ERROR: can't read table data\n");
    4830    exit (1);
    4931  }
    5032
     33  // ohana_memcheck_func (TRUE);
     34
    5135  Coords coords;
    52   if (!GetCoords (&coords, &table.header)) {
     36  if (!GetCoords (&coords, header)) {
    5337    fprintf (stderr, "unable to read header WCS info\n");
    5438    exit (3);
    5539  }
     40
     41  // longitude and latitude for UKIRT:
     42  // Latitude = 19.822433, Longitude = -155.470289 W (from google maps)
     43  // longitude = 155.470289 deg = 10.3646859333 h
     44
     45  double Longitude = 10.3646859333;     // hours (+ = W)
     46  double Latitude  = 19.822433; // degrees
    5647
    5748  // DETECTID tells us which chip supplied the data
    … …  
    7970
    8071  int knownFilter = FALSE;
     72  if (!strcmp (filter, "Z")) { knownFilter = TRUE; }
    8173  if (!strcmp (filter, "Y")) { knownFilter = TRUE; }
    8274  if (!strcmp (filter, "J")) { knownFilter = TRUE; }
    … …  
    8880  }
    8981 
    90   sprintf (photcodeName, "UKIRT.%s.%d", filter, detID);
     82  char photcodeName[64];
     83  snprintf (photcodeName, 64, "UKIRT.%s.%02d", filter, detID);
    9184  PhotCode *code = GetPhotcodebyName (photcodeName);
    9285  if (code == NULL) {
    … …  
    9689  int photcode = code[0].code;
    9790
     91  float zeropt;
    9892  if (!gfits_scan (header, "MAGZPT", "%f", 1, &zeropt)) {
    9993    fprintf (stderr, "ERROR:  missing MAGZPT\n");
    10094    exit (3);
    10195  }
    102 
    103   gfits_scan (header, "MJD_OBS", "%lf", 1, &mjd);
    104   tzero = ohana_mjd_to_sec (mjd);
     96  // Mrel = Minst + K*(airmass - 1.0) + zp_nom - Mcal
     97  // Mrel = Minst + K*(airmass - 1.0) + zp
     98  // zp = zp_nom - Mcal
     99  // Mcal = zp_nom - zp
     100  float Mcal = code[0].C*0.001 - zeropt;
     101
     102  float ApCor4, ApCor7;
     103  if (!gfits_scan (header, "APCOR4", "%f", 1, &ApCor4)) {
     104    fprintf (stderr, "ERROR:  missing APCOR4\n");
     105    exit (3);
     106  }
     107  if (!gfits_scan (header, "APCOR7", "%f", 1, &ApCor7)) {
     108    fprintf (stderr, "ERROR:  missing APCOR7\n");
     109    exit (3);
     110  }
     111
     112  // ohana_memcheck_func (TRUE);
     113
     114  double mjd;
     115  gfits_scan (header, "MJD-OBS", "%lf", 1, &mjd);
     116  time_t tzero = ohana_mjd_to_sec (mjd);
    105117  double jd = ohana_sec_to_jd (tzero);
    106   double sidtime  = 15.0*ohana_lst (jd, Longitude); // sidtime in degrees
     118
     119  double sidtime  = ohana_lst (jd, Longitude); // sidtime in hours
     120  double sidtime_deg = 15.0*sidtime; // sidtime in degrees
    107121
    108122  float seeing;
    109123  gfits_scan (header, "SEEING", "%f", 1, &seeing);
    110  
    111124
    112125  float exptime;
    … …  
    114127  float magtime = 2.5*log10(exptime);
    115128
    116   ZeroPt = GetZeroPoint();
     129  float ZeroPt = GetZeroPoint();
     130
     131  // ohana_memcheck_func (TRUE);
    117132
    118133  // create a Star entry for each filter and detection
    119   Nstars = table.header[0].Naxis[1];
    120   ALLOCATE (stars, Stars, NFILTER*Nstars);
     134  int Nstars = table.header[0].Naxis[1];
     135  ALLOCATE (stars, Stars, Nstars);
    121136
    122137  GET_COLUMN (X_coordinate, float);
    … …  
    126141  GET_COLUMN (Y_coordinate_err, float);
    127142
    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);
     143  GET_COLUMN (Aper_flux_4, float); // 1 arcsec aperture
     144  GET_COLUMN (Aper_flux_4_err, float);
     145
     146  GET_COLUMN (Aper_flux_7, float); // 4 arcsec aperture
     147  GET_COLUMN (Aper_flux_7_err, float);
    133148
    134149  GET_COLUMN (Kron_flux, float);
    … …  
    138153  GET_COLUMN (Sky_rms, float);
    139154
    140   // the value of stars[].M is supposed to be the instrumental magnitude offset by the
    141   // default zero point 25.0 (-2.5*log_10(counts/sec) + ZeroPt).  The magnitude reported
    142   // by SDSS is the calibrated mag: -2.5*log_10(counts/sec) + C_0.  Adjust magnitudes to
    143   // compensate for the difference.
    144 
     155  double R, D, alt, az;
     156
     157  // ohana_memcheck_func (TRUE);
     158
     159  int i;
    145160  for (i = 0; i < Nstars; i++) {
    146       InitStar (&stars[i]);
     161    InitStar (&stars[i]);
    147162     
    148       // any values not explicitly set are left at 0.0
    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
    178       // latitude = 32.7803611755 deg
    179 
    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
    189     }
    190   }   
    191 
    192   for (i = 0; i < NFILTER; i++) {
    193 
    194     N = i + *nimages;
     163    // any values not explicitly set are left at 0.0
     164    XY_to_RD (&R, &D, X_coordinate[i], Y_coordinate[i], &coords);
     165
     166    stars[i].average.R         = R;
     167    stars[i].average.D         = D;
     168    stars[i].average.dR        = NAN;
     169    stars[i].average.dD        = NAN;
     170
     171    stars[i].measure.Xccd      = X_coordinate[i];
     172    stars[i].measure.Yccd      = Y_coordinate[i];
     173    stars[i].measure.dXccd     = ToShortPixels(X_coordinate_err[i]);
     174    stars[i].measure.dYccd     = ToShortPixels(Y_coordinate_err[i]);
     175
     176    stars[i].measure.M         = -2.5*log10(Aper_flux_4[i]) + ZeroPt + magtime + ApCor4;
     177    stars[i].measure.dM        = Aper_flux_4_err[i] / Aper_flux_4[i];
     178    stars[i].measure.Map       = -2.5*log10(Aper_flux_7[i]) + ZeroPt + magtime + ApCor7;
     179    stars[i].measure.dMap      = Aper_flux_7_err[i] / Aper_flux_7[i];
     180    stars[i].measure.Mkron     = -2.5*log10(Kron_flux[i]) + ZeroPt + magtime;
     181    stars[i].measure.dMkron    = Kron_flux_err[i] / Kron_flux[i];
     182    stars[i].measure.Sky       = Sky_level[i]; // adjust this to counts?
     183    stars[i].measure.dSky      = Sky_rms[i];
     184    stars[i].measure.FWx       = ToShortPixels(seeing); // reported in arcsec?
     185    stars[i].measure.FWy       = ToShortPixels(seeing);
     186
     187    stars[i].measure.Mcal      = Mcal;
     188
     189    stars[i].measure.detID     = i;
     190    stars[i].measure.t         = tzero; // time since row 0
     191    stars[i].measure.dt        = magtime; // 2.5 * log(exptime)
     192
     193    altaz (&alt, &az, sidtime_deg - stars[i].average.R, stars[i].average.D, Latitude);
     194
     195    stars[i].measure.airmass   = 1.0 / dCOS(90.0 - alt);
     196    stars[i].measure.az        = az;
     197    stars[i].measure.photcode  = photcode;
     198    stars[i].measure.imageID   = *nimages; // set imageID to entry for this filter
     199  }
     200
     201  if (1) {
     202    int N = *nimages;
    195203   
    196     // XXX for now, we define a totally fake coordinate system centered on the first listed star
    197     InitCoords (&images[N].coords, "DEC--TAN");
    198     images[N].coords.crval1 = stars[0].average.R;
    199     images[N].coords.crval2 = stars[0].average.D;
    200     images[N].coords.crpix1 = stars[0].measure.Xccd;
    201     images[N].coords.crpix2 = stars[0].measure.Yccd;
    202     images[N].coords.cdelt1 = images[N].coords.cdelt2 = 0.4 / 3600.0;
    203 
    204     images[N].NX = 2048;
    205     images[N].NY = 1490;
    206 
    207     images[N].tzero = tzero[i];
     204    // ohana_memcheck_func (TRUE);
     205
     206    images[N].coords = coords;
     207
     208    int Nx, Ny;
     209    gfits_scan (header, "NXOUT", "%d", 1, &Nx);
     210    gfits_scan (header, "NYOUT", "%d", 1, &Ny);
     211
     212    images[N].NX = Nx;
     213    images[N].NY = Ny;
     214 
     215    images[N].tzero = tzero;
    208216    images[N].cerror = 0.0;
    209217 
    210218    // set photcodes for the 5 images (SDSS_U,G,R,I,Z)
    211     images[N].photcode = photcode[i];
     219    images[N].photcode = photcode;
    212220
    213221    // calculate this from : C_OBS, TRACKING, and NY
    214     images[N].exptime = 2048*clockRate;
     222    images[N].exptime = exptime;
    215223 
    216224    images[N].apmifit = 0.0;
    … …  
    218226    images[N].detection_limit = 0.0;
    219227    images[N].saturation_limit = 0.0;
    220     images[N].fwhm_x = seeing[i];
    221     images[N].fwhm_y = seeing[i];
    222 
    223     // XXX longitude and latitude are known for SDSS
    224     // SDSS is at : Latitude 32° 46' 49.30" N, Longitude 105° 49' 13.50" W
    225     // longitude = 105.820419312 deg = 7.05469417
    226     // latitude = 32.7803611755 deg
    227 
    228     double Longitude = 7.05469417;   // hours (+ = W)
    229     double Latitude = 32.7803611755; // degrees
    230 
    231     jd = ohana_sec_to_jd (images[N].tzero);
    232     images[N].sidtime  = ohana_lst (jd, Longitude);
     228    images[N].fwhm_x = seeing * 25.0 * images[N].coords.cdelt1 * 3600.0;
     229    images[N].fwhm_y = seeing * 25.0 * images[N].coords.cdelt1 * 3600.0;
     230
     231    images[N].sidtime  = sidtime;
    233232    images[N].latitude = Latitude;
    234233
    235     altaz (&alt, &az, 15.0*images[N].sidtime - images[N].coords.crval1, images[N].coords.crval2, Latitude);
    236 
    237     images[N].trate = clockRate * 1e-4;
    238     images[N].secz = stars[0].measure.airmass;
    239     images[N].ccdnum = camcol;
     234    XY_to_RD (&R, &D, 0.5*images[N].NX, 0.5*images[N].NY, &coords);
     235    altaz (&alt, &az, sidtime_deg - R, D, Latitude);
     236
     237    images[N].trate = 0.0;
     238    images[N].secz = 1.0 / sin (RAD_DEG * alt);
     239    images[N].ccdnum = detID;
    240240
    241241    // secz is in units milli-airmass
    242     images[N].Mcal = 0.0;
     242    images[N].Mcal = Mcal;
    243243    images[N].Xm   = NAN_S_SHORT;
    244244    images[N].flags = 0;
    … …  
    252252
    253253    // save the filename
    254     snprintf (images[N].name, DVO_IMAGE_NAME_LEN, "%s[%s]", name, filtname[i]);
    255   }
    256 
    257   NSTAR_GROUP = NFILTER;
    258   *nimages += NFILTER;
    259   *nstars = Nstars*NFILTER;
     254    snprintf (images[N].name, DVO_IMAGE_NAME_LEN, "%s[%02d]", name, detID);
     255    *nimages = N + 1;
     256  }
     257 
     258  // ohana_memcheck_func (TRUE);
     259
     260  *nstars = Nstars;
    260261  return (stars);
    261262}

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

    r37357 r38184  
    1919    dvo_image_create (&db, GetZeroPoint());
    2020    ImageIndexFileInit ();
    21   isEmpty = TRUE;
     21    isEmpty = TRUE;
    2222  } else {
    2323    /* position to start of file */
Note: See TracChangeset for help on using the changeset viewer.