Index: /branches/eam_branches/ipp-20150419/Ohana/src/addstar/include/addstar.h =================================================================== --- /branches/eam_branches/ipp-20150419/Ohana/src/addstar/include/addstar.h (revision 38166) +++ /branches/eam_branches/ipp-20150419/Ohana/src/addstar/include/addstar.h (revision 38167) @@ -77,5 +77,5 @@ enum {ADDSTAR_MODE_NONE, ADDSTAR_MODE_IMAGE, ADDSTAR_MODE_REFLIST, ADDSTAR_MODE_REFCAT, ADDSTAR_MODE_FAKEIMAGE, ADDSTAR_MODE_RESORT, ADDSTAR_MODE_CREATE_ID}; -enum {NONE, SIMPLE_CMP, SIMPLE_CMF, SIMPLE_MEF, MOSAIC_CMP, MOSAIC_CMF, MOSAIC_MEF, MOSAIC_PHU, SDSS_OBJ}; +enum {NONE, SIMPLE_CMP, SIMPLE_CMF, SIMPLE_MEF, MOSAIC_CMP, MOSAIC_CMF, MOSAIC_MEF, MOSAIC_PHU, SDSS_OBJ, UKIRT_OBJ}; /* note: MEF implies CMF */ Index: /branches/eam_branches/ipp-20150419/Ohana/src/addstar/src/GetFileMode.c =================================================================== --- /branches/eam_branches/ipp-20150419/Ohana/src/addstar/src/GetFileMode.c (revision 38166) +++ /branches/eam_branches/ipp-20150419/Ohana/src/addstar/src/GetFileMode.c (revision 38167) @@ -24,4 +24,14 @@ if (VERBOSE) fprintf (stderr, "found SDSS objects\n"); return SDSS_OBJ; + } + + + // UKIRT files have TELESCOP = UKIRT and INSTRUME = WFCAM: + char telescope[81], instrument[81]; + int have_TELESCOPE = gfits_scan (header, "TELESCOP", "%s", 1, telescope); + int have_INSTRUMENT = gfits_scan (header, "INSTRUME", "%s", 1, instrument); + if (haveTELESCOPE && have_INSTRUMENT && !strcmp(telescope, "UKIRT"), && !strcmp(instrument, "WFCAM")) { + if (VERBOSE) fprintf (stderr, "found UKIRT data\n"); + return UKIRT_OBJ; } Index: /branches/eam_branches/ipp-20150419/Ohana/src/addstar/src/LoadDataUKIRT.c =================================================================== --- /branches/eam_branches/ipp-20150419/Ohana/src/addstar/src/LoadDataUKIRT.c (revision 38167) +++ /branches/eam_branches/ipp-20150419/Ohana/src/addstar/src/LoadDataUKIRT.c (revision 38167) @@ -0,0 +1,48 @@ +# include "addstar.h" + +// 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 LoadDataSDSS (FILE *f, char *file, Image **images, off_t *nvalid, Stars **stars, unsigned int *Nstars, Header **headers, off_t *extsize, HeaderSet *headerSets, off_t Nimages) { + + off_t Nskip, Nvalid, NVALID; + char *name; + int j, Nhead, Ndata; + unsigned int Ninstars; + Stars *inStars; + + if (images[0] == NULL) { + Nvalid = 0; + NVALID = 4; + ALLOCATE (images[0], Image, NVALID); + } else { + Nvalid = *nvalid; + NVALID = Nvalid + 4; + REALLOCATE (images[0], Image, NVALID); + } + + // find image rootname + name = filebasename (file); + + // validate the number of headers sets == 4 + for (i = 0; i < 4; i++) { + if (VERBOSE) fprintf (stderr, "reading header for %s (%s)\n", headerSets[i].exthead, headerSets[i].extdata); + + // advance the pointer to the start of the corresponding table block + Nhead = headerSets[0].extnum_head; + Ndata = headerSets[i].extnum_data; + Nskip = 0; + for (j = 0; j < Ndata; j++) { + Nskip += extsize[j]; + } + fseeko (f, Nskip, SEEK_SET); + + inStars = ReadStarsUKIRT (f, name, headers[Nhead], headers[Ndata], images[0], &Nvalid, &Ninstars); + *stars = MergeStars (*stars, Nstars, inStars, Ninstars); + + free (name); + *nvalid = Nvalid; + } + return (TRUE); +} + Index: /branches/eam_branches/ipp-20150419/Ohana/src/addstar/src/LoadHeaders.c =================================================================== --- /branches/eam_branches/ipp-20150419/Ohana/src/addstar/src/LoadHeaders.c (revision 38166) +++ /branches/eam_branches/ipp-20150419/Ohana/src/addstar/src/LoadHeaders.c (revision 38167) @@ -24,4 +24,5 @@ if (i == 0) { *mode = GetFileMode (headers[0]); + // CMP mode has no binary table, just image header and text blocks if ((*mode == SIMPLE_CMP) || (*mode == MOSAIC_CMP)) { *Nheaders = i; Index: /branches/eam_branches/ipp-20150419/Ohana/src/addstar/src/LoadStars.c =================================================================== --- /branches/eam_branches/ipp-20150419/Ohana/src/addstar/src/LoadStars.c (revision 38166) +++ /branches/eam_branches/ipp-20150419/Ohana/src/addstar/src/LoadStars.c (revision 38167) @@ -100,4 +100,10 @@ } + // if these are SDSS data, load with SDSS-specific wrapper + if (headerSets[0].exttype && !strcmp (headerSets[0].exttype, "UKIRT_OBJ")) { + LoadDataUKIRT (f, file[i], images, Nimages, &stars, Nstars, headers, extsize, headerSets, NheaderSets); + continue; + } + LoadData (f, file[i], images, Nimages, &stars, Nstars, headers, extsize, headerSets, NheaderSets, options); Index: /branches/eam_branches/ipp-20150419/Ohana/src/addstar/src/MatchHeaders.c =================================================================== --- /branches/eam_branches/ipp-20150419/Ohana/src/addstar/src/MatchHeaders.c (revision 38166) +++ /branches/eam_branches/ipp-20150419/Ohana/src/addstar/src/MatchHeaders.c (revision 38167) @@ -39,4 +39,35 @@ } + // check for ukirt mode + if (Nheaders == 5) { + int myMode = GetFileMode (headers[1]); + if (myMode != UKIRT_OBJ) goto other_modes; + for (i = 2; i < 5; i++) { + myMode = GetFileMode (headers[i]); + if (myMode != UKIRT_OBJ) { + fprintf (stderr, "inconsistent headers for UKIRT\n"); + exit (2); + } + } + + for (i = 1; i < Nheaders; i++) { + extsize[0][i] = headers[i][0].datasize + gfits_data_size (headers[i]); + + gfits_scan (headers[i], ExtnameKeyword, "%s", 1, extname); + headerSets[Nimage].exttype = strcreate ("UKIRT_OBJ"); + headerSets[Nimage].extdata = strcreate (extname); + headerSets[Nimage].exthead = strcreate (extname); + headerSets[Nimage].extxrad = NULL; + headerSets[Nimage].extnum_data = i; + headerSets[Nimage].extnum_head = i; + headerSets[Nimage].extnum_xrad = -1; + Nimage ++; + } + *nimage = Nimage; + return headerSets; + } + +other_modes: + // now examine the headers, count the table entries, find corresponding headers for (i = 0; i < Nheaders; i++) { @@ -71,4 +102,5 @@ if (!strcmp (exttype, "PS1_DV4")) goto keep; if (!strcmp (exttype, "PS1_DV5")) goto keep; + continue; Index: /branches/eam_branches/ipp-20150419/Ohana/src/addstar/src/ReadStarsUKIRT.c =================================================================== --- /branches/eam_branches/ipp-20150419/Ohana/src/addstar/src/ReadStarsUKIRT.c (revision 38167) +++ /branches/eam_branches/ipp-20150419/Ohana/src/addstar/src/ReadStarsUKIRT.c (revision 38167) @@ -0,0 +1,303 @@ +# 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) \ + 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 == 1); + +/* grab named photcode */ +# define NAMED_PHOTCODE_AND_ZP(CODE,ZP,NAME) { \ + PhotCode *code; \ + code = GetPhotcodebyName (NAME); \ + if (code == NULL) { \ + fprintf (stderr, "ERROR: photcode %s not found in photcode table\n", NAME); \ + exit (0); } \ + CODE = code[0].code; \ + ZP = 0.001*code[0].C; } + +// XXX NOTE : as of 2008.02.27, the zero point is still carried internally in millimags + +// given a file with the pointer at the start of the table block and the +// corresponding image header, load the stars from the table +Stars *ReadStarsSDSS (FILE *f, char *name, Header *header, Header *in_theader, Image *images, off_t *nimages, unsigned int *nstars) { + + off_t Nskip, Nrow; + int i, j, N, Nstars, camcol; + char type[80]; + Header theader; + FTable table; + Stars *stars; + double clockRate, mjd[5], jd, sidtime, alt, az; + float seeing[5], photErr[5], zeropt[5], ZeroPt; + time_t tzero[5]; + char filtname[16][5]; + int photcode[5]; + int Ncol; // used in the GET_COLUMN_1,5 macros above + + if (in_theader == NULL) { + table.header = &theader; + if (!gfits_fread_header (f, table.header)) Shutdown ("ERROR: can't read table header"); + } else { + table.header = in_theader; + Nskip = in_theader[0].datasize; + fseeko (f, Nskip, SEEK_CUR); + } + + /* load the table data */ + if (!gfits_fread_ftable_data (f, &table, FALSE)) { + fprintf (stderr, "ERROR: can't read table header\n"); + exit (1); + } + + // + int hasDetID = gfits_scan (header, "DETECTID", "%s", 1, detectID); + if (!hasDetID) { + 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 + + ZeroPt = GetZeroPoint(); + + // create a Star entry for each filter and detection + Nstars = table.header[0].Naxis[1]; + 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); + + // the value of stars[].M is supposed to be the instrumental magnitude offset by the + // default zero point 25.0 (-2.5*log_10(counts/sec) + ZeroPt). The magnitude reported + // by SDSS is the calibrated mag: -2.5*log_10(counts/sec) + C_0. Adjust magnitudes to + // compensate for the difference. + + for (i = 0; i < Nstars; i++) { + for (j = 0; j < NFILTER; j++) { + N = NFILTER*i + j; + InitStar (&stars[N]); + + // 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 + // 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 + } + } + + for (i = 0; i < NFILTER; i++) { + + N = i + *nimages; + + // XXX for now, we define a totally fake coordinate system centered on the first listed star + InitCoords (&images[N].coords, "DEC--TAN"); + images[N].coords.crval1 = stars[0].average.R; + images[N].coords.crval2 = stars[0].average.D; + images[N].coords.crpix1 = stars[0].measure.Xccd; + images[N].coords.crpix2 = stars[0].measure.Yccd; + images[N].coords.cdelt1 = images[N].coords.cdelt2 = 0.4 / 3600.0; + + images[N].NX = 2048; + images[N].NY = 1490; + + images[N].tzero = tzero[i]; + images[N].cerror = 0.0; + + // set photcodes for the 5 images (SDSS_U,G,R,I,Z) + images[N].photcode = photcode[i]; + + // calculate this from : C_OBS, TRACKING, and NY + images[N].exptime = 2048*clockRate; + + images[N].apmifit = 0.0; + images[N].dapmifit = 0.0; + images[N].detection_limit = 0.0; + images[N].saturation_limit = 0.0; + images[N].fwhm_x = seeing[i]; + images[N].fwhm_y = seeing[i]; + + // XXX longitude and latitude are known for SDSS + // SDSS is at : Latitude 32° 46' 49.30" N, Longitude 105° 49' 13.50" W + // longitude = 105.820419312 deg = 7.05469417 + // latitude = 32.7803611755 deg + + double Longitude = 7.05469417; // hours (+ = W) + double Latitude = 32.7803611755; // degrees + + jd = ohana_sec_to_jd (images[N].tzero); + images[N].sidtime = ohana_lst (jd, Longitude); + images[N].latitude = Latitude; + + altaz (&alt, &az, 15.0*images[N].sidtime - images[N].coords.crval1, images[N].coords.crval2, Latitude); + + images[N].trate = clockRate * 1e-4; + images[N].secz = stars[0].measure.airmass; + images[N].ccdnum = camcol; + + // secz is in units milli-airmass + images[N].Mcal = 0.0; + images[N].Xm = NAN_S_SHORT; + images[N].flags = 0; + + images[N].nstar = Nstars; + + images[N].imageID = N; + images[N].externID = 0; + images[N].sourceID = 0; + images[N].parentID = UINT32_MAX; // UpdateImageIDs sets parentID = 0 + + // save the filename + snprintf (images[N].name, DVO_IMAGE_NAME_LEN, "%s[%s]", name, filtname[i]); + } + + NSTAR_GROUP = NFILTER; + *nimages += NFILTER; + *nstars = Nstars*NFILTER; + 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); + +}