Index: /branches/eam_branches/ipp-20150419/Ohana/src/addstar/src/LoadDataUKIRT.c =================================================================== --- /branches/eam_branches/ipp-20150419/Ohana/src/addstar/src/LoadDataUKIRT.c (revision 38178) +++ /branches/eam_branches/ipp-20150419/Ohana/src/addstar/src/LoadDataUKIRT.c (revision 38179) @@ -2,6 +2,4 @@ // examine the header sets and set the Image entries for the the valid images -// there should only be a single data set (phu + table) in this file -// each SDSS data set corresponds to 5 images (ugriz) int 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) { @@ -15,9 +13,9 @@ if (images[0] == NULL) { Nvalid = 0; - NVALID = 4; + NVALID = NheaderSets; ALLOCATE (images[0], Image, NVALID); } else { Nvalid = *nimages; - NVALID = Nvalid + 4; + NVALID = Nvalid + NheaderSets; REALLOCATE (images[0], Image, NVALID); } Index: /branches/eam_branches/ipp-20150419/Ohana/src/addstar/src/ReadStarsUKIRT.c =================================================================== --- /branches/eam_branches/ipp-20150419/Ohana/src/addstar/src/ReadStarsUKIRT.c (revision 38178) +++ /branches/eam_branches/ipp-20150419/Ohana/src/addstar/src/ReadStarsUKIRT.c (revision 38179) @@ -1,18 +1,5 @@ # include "addstar.h" -int SetSDSSFlags (Stars *star, unsigned int flags1, unsigned int flags2); - -# define NFILTER 5 - -# define GET_COLUMN_5(NAME,TYPE) \ - TYPE *NAME; \ - NAME = (TYPE *) gfits_get_bintable_column_data (table.header, &table, #NAME, type, &Nrow, &Ncol); \ - assert (NAME); assert (!strcmp (type, #TYPE)); assert (Nrow == Nstars); assert (Ncol == NFILTER); - -# define GET_COLUMN_5_NOASSERT(NAME,TYPE) \ - TYPE *NAME; \ - NAME = (TYPE *) gfits_get_bintable_column_data (table.header, &table, #NAME, type, &Nrow, &Ncol); - -# define GET_COLUMN_1(NAME,TYPE) \ +# define GET_COLUMN(NAME,TYPE) \ TYPE *NAME; \ NAME = (TYPE *) gfits_get_bintable_column_data (table.header, &table, #NAME, type, &Nrow, &Ncol); \ @@ -46,4 +33,5 @@ int Ncol; // used in the GET_COLUMN_1,5 macros above + // XXX ??? what is this for? if (in_theader == NULL) { table.header = &theader; @@ -61,63 +49,68 @@ } - // - int hasDetID = gfits_scan (header, "DETECTID", "%s", 1, detectID); - if (!hasDetID) { + Coords coords; + if (!GetCoords (&coords, &table.header)) { + fprintf (stderr, "unable to read header WCS info\n"); + exit (3); + } + + // DETECTID tells us which chip supplied the data + char detIDstr[80]; + if (!gfits_scan (header, "DETECTID", "%s", 1, detIDstr)) { fprintf (stderr, "missing DETECTID\n"); exit (3); } - - -w20050611_00700_st_cat.fits[APM-BINARYTABLE] RSC:H2:60 -w20050611_00700_st_cat.fits[APM-BINARYTABLE] RSC:H2:63 -w20050611_00700_st_cat.fits[APM-BINARYTABLE] RSC:H2:76 -w20050611_00700_st_cat.fits[APM-BINARYTABLE] RSC:H2:41 - - strcpy (filtname[0], "u"); - strcpy (filtname[1], "g"); - strcpy (filtname[2], "r"); - strcpy (filtname[3], "i"); - strcpy (filtname[4], "z"); - - NAMED_PHOTCODE_AND_ZP (photcode[0], zeropt[0], "U_SDSS"); - NAMED_PHOTCODE_AND_ZP (photcode[1], zeropt[1], "G_SDSS"); - NAMED_PHOTCODE_AND_ZP (photcode[2], zeropt[2], "R_SDSS"); - NAMED_PHOTCODE_AND_ZP (photcode[3], zeropt[3], "I_SDSS"); - NAMED_PHOTCODE_AND_ZP (photcode[4], zeropt[4], "Z_SDSS"); - - // XXXYYYZZZ SDSS tables have special flags for undefined/unmeasured values and errors: - // -9999 flags unmeasured values, and the corresponding error may or may not be meaningful - // -1000 flags errors that are not determined, even though the corresponding quantity is. - // These special values need to be trapped here to avoid averaging meaningful numbers with the flag values. - - // various header values needed to calculate per-star data below - gfits_scan (header, "C_OBS", "%lf", 1, &clockRate); // value in header is usec / unbinned row - clockRate *= 1e-6; // convert to seconds / unbinned row - - gfits_scan (table.header, "MJD_U", "%lf", 1, &mjd[0]); - gfits_scan (table.header, "MJD_G", "%lf", 1, &mjd[1]); - gfits_scan (table.header, "MJD_R", "%lf", 1, &mjd[2]); - gfits_scan (table.header, "MJD_I", "%lf", 1, &mjd[3]); - gfits_scan (table.header, "MJD_Z", "%lf", 1, &mjd[4]); - tzero[0] = ohana_mjd_to_sec (mjd[0]); - tzero[1] = ohana_mjd_to_sec (mjd[1]); - tzero[2] = ohana_mjd_to_sec (mjd[2]); - tzero[3] = ohana_mjd_to_sec (mjd[3]); - tzero[4] = ohana_mjd_to_sec (mjd[4]); - - gfits_scan (table.header, "SEEING_U", "%f", 1, &seeing[0]); - gfits_scan (table.header, "SEEING_G", "%f", 1, &seeing[1]); - gfits_scan (table.header, "SEEING_R", "%f", 1, &seeing[2]); - gfits_scan (table.header, "SEEING_I", "%f", 1, &seeing[3]); - gfits_scan (table.header, "SEEING_Z", "%f", 1, &seeing[4]); - - gfits_scan (table.header, "PSFERR_U", "%f", 1, &photErr[0]); - gfits_scan (table.header, "PSFERR_G", "%f", 1, &photErr[1]); - gfits_scan (table.header, "PSFERR_R", "%f", 1, &photErr[2]); - gfits_scan (table.header, "PSFERR_I", "%f", 1, &photErr[3]); - gfits_scan (table.header, "PSFERR_Z", "%f", 1, &photErr[4]); - - gfits_scan (header, "CAMCOL", "%d", 1, &camcol); // value in header is usec / unbinned row + int detID = 0; + if (!strcmp (detIDstr, "RSC:H2:60")) { detID = 1; } + if (!strcmp (detIDstr, "RSC:H2:63")) { detID = 2; } + if (!strcmp (detIDstr, "RSC:H2:76")) { detID = 3; } + if (!strcmp (detIDstr, "RSC:H2:41")) { detID = 4; } + if (!detID) { + fprintf (stderr, "unknown DETECTID %s\n", detIDstr); + exit (3); + } + + char filter[80]; + if (!gfits_scan (header, "FILTER", "%s", 1, filter)) { + fprintf (stderr, "missing FILTER\n"); + exit (3); + } + + int knownFilter = FALSE; + if (!strcmp (filter, "Y")) { knownFilter = TRUE; } + if (!strcmp (filter, "J")) { knownFilter = TRUE; } + if (!strcmp (filter, "H")) { knownFilter = TRUE; } + if (!strcmp (filter, "K")) { knownFilter = TRUE; } + if (!knownFilter) { + fprintf (stderr, "unknown FILTER %s\n", filter); + exit (3); + } + + sprintf (photcodeName, "UKIRT.%s.%d", filter, detID); + PhotCode *code = GetPhotcodebyName (photcodeName); + if (code == NULL) { + fprintf (stderr, "ERROR: photcode %s not found in photcode table\n", photcodeName); + exit (3); + } + int photcode = code[0].code; + + if (!gfits_scan (header, "MAGZPT", "%f", 1, &zeropt)) { + fprintf (stderr, "ERROR: missing MAGZPT\n"); + exit (3); + } + + gfits_scan (header, "MJD_OBS", "%lf", 1, &mjd); + tzero = ohana_mjd_to_sec (mjd); + double jd = ohana_sec_to_jd (tzero); + double sidtime = 15.0*ohana_lst (jd, Longitude); // sidtime in degrees + + float seeing; + gfits_scan (header, "SEEING", "%f", 1, &seeing); + + + float exptime; + gfits_scan (header, "EXP_TIME", "%f", 1, &exptime); + float magtime = 2.5*log10(exptime); ZeroPt = GetZeroPoint(); @@ -127,27 +120,21 @@ ALLOCATE (stars, Stars, NFILTER*Nstars); - GET_COLUMN_5 (rowc, float); - GET_COLUMN_5 (colc, float); - GET_COLUMN_5 (sky, float); - GET_COLUMN_5 (psfCounts, float); - GET_COLUMN_5 (fiberCounts, float); - GET_COLUMN_5 (offsetRa, float); - GET_COLUMN_5 (offsetDec, float); - GET_COLUMN_5 (flags, int); - GET_COLUMN_5 (flags2, int); - -#ifdef notyet - GET_COLUMN_5 (prob_psf, float); -#else - GET_COLUMN_5_NOASSERT (prob_psf, float); -#endif - - GET_COLUMN_1 (ra, double); - GET_COLUMN_1 (dec, double); - - GET_COLUMN_5 (rowcErr, float); - GET_COLUMN_5 (colcErr, float); - GET_COLUMN_5 (skyErr, float); - GET_COLUMN_5 (psfCountsErr, float); + GET_COLUMN (X_coordinate, float); + GET_COLUMN (Y_coordinate, float); + + GET_COLUMN (X_coordinate_err, float); + GET_COLUMN (Y_coordinate_err, float); + + GET_COLUMN (Aper_flux_8, float); + GET_COLUMN (Aper_flux_8_err, float); + + GET_COLUMN (Aper_flux_5, float); + GET_COLUMN (Aper_flux_5_err, float); + + GET_COLUMN (Kron_flux, float); + GET_COLUMN (Kron_flux_err, float); + + GET_COLUMN (Sky_level, float); + GET_COLUMN (Sky_rms, float); // the value of stars[].M is supposed to be the instrumental magnitude offset by the @@ -157,52 +144,47 @@ for (i = 0; i < Nstars; i++) { - for (j = 0; j < NFILTER; j++) { - N = NFILTER*i + j; - InitStar (&stars[N]); + InitStar (&stars[i]); // any values not explicitly set are left at 0.0 - stars[N].average.R = ra[i] + dCOS(dec[i]) * offsetRa[N] / 3600.0; - stars[N].average.D = dec[i] + offsetDec[N] / 3600.0; - stars[N].average.dR = NAN; - stars[N].average.dD = NAN; - - stars[N].measure.Xccd = colc[N]; - stars[N].measure.Yccd = rowc[N]; - stars[N].measure.dXccd = ToShortPixels(colcErr[N]); - stars[N].measure.dYccd = ToShortPixels(rowcErr[N]); - stars[N].measure.M = psfCounts[N] + ZeroPt - zeropt[j]; - stars[N].measure.dM = psfCountsErr[N]; - stars[N].measure.Map = fiberCounts[N] + ZeroPt - zeropt[j]; - stars[N].measure.Sky = sky[N]; // adjust this to counts? - stars[N].measure.dSky = skyErr[N]; - stars[N].measure.FWx = ToShortPixels(seeing[j]); // reported in arcsec? - stars[N].measure.FWy = ToShortPixels(seeing[j]); - if (prob_psf) { - stars[N].measure.psfChisq = prob_psf[N]; // XXX not really the correct value... - } else { - stars[N].measure.psfChisq = NAN; - } - stars[N].measure.detID = N; - stars[N].measure.t = tzero[j] + clockRate*rowc[N]; // time since row 0 - stars[N].measure.dt = 4.32912209; // 2.5 * log(53.907456) the sdss exposure time // old comment is 53907456 is this 2048*clockRate ? - - SetSDSSFlags (&stars[N], flags[N], flags2[N]); - - // longitude and latitude for SDSS: - // Latitude 32° 46' 49.30" N, Longitude 105° 49' 13.50" W - // longitude = 105.820419312 deg = 7.05469417 + XY_to_RD (&R, &D, X_coordinate[i], Y_coordinate[i], &coords); + + stars[i].average.R = R; + stars[i].average.D = D; + stars[i].average.dR = NAN; + stars[i].average.dD = NAN; + + stars[i].measure.Xccd = X_coordinate[i]; + stars[i].measure.Yccd = Y_coordinate[i]; + stars[i].measure.dXccd = ToShortPixels(X_coordinate_err[i]); + stars[i].measure.dYccd = ToShortPixels(Y_coordinate_err[i]); + stars[i].measure.M = -2.5*log10(Aper_flux_8[i]) + ZeroPt - zeropt; + stars[i].measure.dM = Aper_flux_8_err[i] / Aper_flux_8[i]; + stars[i].measure.Map = -2.5*log10(Aper_flux_8[i]) + ZeroPt - zeropt; + stars[i].measure.dMap = Aper_flux_5_err[i] / Aper_flux_5[i]; + stars[i].measure.Mkron = -2.5*log10(Aper_flux_8[i]) + ZeroPt - zeropt; + stars[i].measure.dMkron = Aper_flux_5_err[i] / Aper_flux_5[i]; + stars[i].measure.Sky = Sky_level[i]; // adjust this to counts? + stars[i].measure.dSky = Sky_rms[i]; + stars[i].measure.FWx = ToShortPixels(seeing); // reported in arcsec? + stars[i].measure.FWy = ToShortPixels(seeing); + + stars[i].measure.detID = i; + stars[i].measure.t = tzero; // time since row 0 + stars[i].measure.dt = magtime; // 2.5 * log(exptime) + + // longitude and latitude for UKIRT: + // Latitude , Longitude ??? + // longitude = 105.820419312 deg = 7.05469417 h // latitude = 32.7803611755 deg - double Longitude = 7.05469417; // hours (+ = W) - double Latitude = 32.7803611755; // degrees - - jd = ohana_sec_to_jd (stars[N].measure.t); - sidtime = 15.0*ohana_lst (jd, Longitude); // sidtime in degrees - altaz (&alt, &az, sidtime - stars[N].average.R, stars[N].average.D, Latitude); - - stars[N].measure.airmass = 1.0 / dCOS(90.0 - alt); - stars[N].measure.az = az; - stars[N].measure.photcode = photcode[j]; - stars[N].measure.imageID = j + *nimages; // set imageID to entry for this filter + double Longitude = 7.05469417 * foo; // hours (+ = W) + double Latitude = 32.7803611755 * fooo; // degrees + + altaz (&alt, &az, sidtime - stars[i].average.R, stars[i].average.D, Latitude); + + stars[i].measure.airmass = 1.0 / dCOS(90.0 - alt); + stars[i].measure.az = az; + stars[i].measure.photcode = photcode; + stars[i].measure.imageID = *nimages; // set imageID to entry for this filter } } @@ -278,24 +260,2 @@ return (stars); } - -int SetSDSSFlags (Stars *star, unsigned int flags1, unsigned int flags2) { - - // XXX this is wrong, need to roll left to set the correct bit - if (flags1 & 0x00000002) star[0].measure.photFlags |= 0x0001; // BRIGHT - 1 1 - if (flags1 & 0x00000004) star[0].measure.photFlags |= 0x0002; // EDGE - 1 2 - if (flags1 & 0x00000008) star[0].measure.photFlags |= 0x0004; // BLENDED - 1 3 - if (flags1 & 0x00000010) star[0].measure.photFlags |= 0x0008; // CHILD - 1 4 - if (flags1 & 0x00000020) star[0].measure.photFlags |= 0x0010; // PEAKCENTER - 1 5 - if (flags1 & 0x00000040) star[0].measure.photFlags |= 0x0020; // NODEBLEND - 1 6 - if (flags1 & 0x00040000) star[0].measure.photFlags |= 0x0040; // SATUR - 1 18 - if (flags1 & 0x00080000) star[0].measure.photFlags |= 0x0080; // NOTCHECKED - 1 19 - if (flags1 & 0x10000000) star[0].measure.photFlags |= 0x0100; // BINNED1 - 1 28 - if (flags1 & 0x20000000) star[0].measure.photFlags |= 0x0200; // BINNED2 - 1 29 - if (flags1 & 0x40000000) star[0].measure.photFlags |= 0x0400; // BINNED4 - 1 30 - if (flags2 & 0x00000040) star[0].measure.photFlags |= 0x0800; // LOCAL_EDGE - 2 7 - if (flags2 & 0x00000800) star[0].measure.photFlags |= 0x1000; // INTERP_CENTER - 2 12 - if (flags2 & 0x00002000) star[0].measure.photFlags |= 0x2000; // DEBLEND_NOPEAK - 2 14 - if (flags2 & 0x02000000) star[0].measure.photFlags |= 0x4000; // NOTCHECKED_CENTER - 2 26 - return (TRUE); - -}