Index: /branches/sj_branches/sj_SDSSaddstar_branch_r23928_20090420/Ohana/src/addstar/src/ReadStarsSDSS.c =================================================================== --- /branches/sj_branches/sj_SDSSaddstar_branch_r23928_20090420/Ohana/src/addstar/src/ReadStarsSDSS.c (revision 24981) +++ /branches/sj_branches/sj_SDSSaddstar_branch_r23928_20090420/Ohana/src/addstar/src/ReadStarsSDSS.c (revision 24982) @@ -67,5 +67,6 @@ strcpy (filtname[4], "z"); - // XXXYYY Why these? and not the 'ref' version? see 'Why these?' below, too. + // XXXJester: SDSS photometry is read with 'dep' photcodes - recalibratable in principle, not standard-star + // photometry (which would get 'ref'). NAMED_PHOTCODE_AND_ZP (photcode[0], zeropt[0], "U_SDSS"); NAMED_PHOTCODE_AND_ZP (photcode[1], zeropt[1], "G_SDSS"); @@ -74,8 +75,13 @@ 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. + // XXXJester:ZZZ 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. However, we can't just skip an entire star + // just because one measurement is missing in one quantity. Perhaps the right thing to do is to replace + // -9999 and -1000 by NaN; however, does the averaging routine (i.e. liststats() in + // Ohana/src/relphot/src/liststats.c, or maybe rather the sort functions in there) correctly ignore NaNs? I don't + // think so. In fact, I think not ignoring NaNs is a generic problem, because NaNs tend to be contagious, + // and just a single bad/missing measurement should not invalidate the average quantity of anything. // various header values needed to calculate per-star data below @@ -94,7 +100,7 @@ tzero[4] = ohana_mjd_to_sec (mjd[4]); - // XXXYYY EDR paper says "psfWidth is is the effective width also determined at the center of each frame. It + // XXXJester: EDR paper says "psfWidth is is the effective width also determined at the center of each frame. It // is a good generic number to quote for the seeing on each frame." It is stored in the tsField file as - // psf_width. + // psf_width. If we read the actual per-field zero points and extinction coefficients, we could also read psf_width. gfits_scan (table.header, "SEEING_U", "%f", 1, &seeing[0]); gfits_scan (table.header, "SEEING_G", "%f", 1, &seeing[1]); @@ -103,5 +109,5 @@ gfits_scan (table.header, "SEEING_Z", "%f", 1, &seeing[4]); - // XXXYYY What's the meaning of photErr? Where is it used again? + // XXXJester: What's the meaning of photErr? Where is it used again? gfits_scan (table.header, "PSFERR_U", "%f", 1, &photErr[0]); gfits_scan (table.header, "PSFERR_G", "%f", 1, &photErr[1]); @@ -112,6 +118,5 @@ gfits_scan (header, "CAMCOL", "%d", 1, &camcol); // value in header is usec / unbinned row - // XXXYYY which number is read here? In other words, where is the zeropoint set that is returned by - // GetZeroPoint()? I can't find the corresponding SetZeroPoint() + // XXXJester: The corresponding SetZeroPoint() is in ConfigInit.c, where ZERO_PT from ~/.ptolemyrc is read ZeroPt = GetZeroPoint(); @@ -147,17 +152,26 @@ // compensate for the difference. - Nstars_used = Nstars; + // XXXYYYJester: the actual C_0 and the airmass coefficient are reported in the tsField file that go with + // every tsObj file (i.e. C_0 = - aa - kk * airmass). Worse, SDSS reports *luptitudes* (asinh magnitudes) + // not conventional ("Pogson") magnitudes, so in fact SDSS reports -(2.5/ln(10))*[asinh((f/f0)/2b)+ln(b)] + // where f/f0 = counts/sec * 10^[0.4*(aa + kk * airmass)]. See + // http://www.sdss.org/dr7/algorithms/fluxcal.html#counts2mag + + Nstars_used = 0; for (i = 0; i < Nstars; i++) { - // Skip objects that are not OK_RUN - if ( ((status[i]) & AR_OBJECT_STATUS_OK_RUN) != 0) { - Nstars_used--; + // Skip objects that are not OK_RUN (including other objects leads to counting the same photons more than + // once - see comment near images[N].NY = ... below) + if ( ((status[i]) & AR_OBJECT_STATUS_OK_RUN) == 0) { continue; } + Nstars_used += 1; for (j = 0; j < NFILTER; j++) { - N = NFILTER*i + j; + N = NFILTER * (Nstars_used-1) + 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; + // XXXYYYJester: the offsetRa here used to be multiplied by dCOS(dec[i]), but that's already been done + // by SDSS + stars[N].average.R = ra[i] + offsetRa[N] / 3600.0; stars[N].average.D = dec[i] + offsetDec[N] / 3600.0; stars[N].average.dR = NAN; @@ -168,5 +182,5 @@ stars[N].measure.dXccd = ShortPixels(colcErr[N]); stars[N].measure.dYccd = ShortPixels(rowcErr[N]); - // XXXYYY again the question what ZeroPt-zeropt[j] evalutes to - according to what it says above, this + // XXXJester: again the question what ZeroPt-zeropt[j] evalutes to - according to what it says above, this // should take .M back into instrumental magnitudes, but does it? In particular, the *actual* zeropoint // that was used for the computation of this calibrated magnitude is in the tsField file; in general, it @@ -176,12 +190,12 @@ 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? XXXYYY It's given in maggies; to get back to counts, + stars[N].measure.Sky = sky[N]; // adjust this to counts? XXXJester: The SDSS Sky is given in maggies/arcsec^2; a maggie is a linear flux measure, 1 maggie corresponds to 0 magnitude, a surface brightness of 20 mag/square arcsec corresponds to 10-8 maggies per square arcsec. stars[N].measure.dSky = skyErr[N]; - stars[N].measure.FWx = ShortPixels(seeing[j]); // reported in arcsec? XXXYYY: Yes! + stars[N].measure.FWx = ShortPixels(seeing[j]); // reported in arcsec? XXXJester:: Yes! stars[N].measure.FWy = ShortPixels(seeing[j]); stars[N].measure.psfChisq = prob_psf[N]; // XXX not really the correct value... stars[N].measure.detID = N; stars[N].measure.t = tzero[j] + clockRate*rowc[N]; // time since row 0 - stars[N].measure.dt = 53.907456; // is this 2048*clockRate ? XXXYYY Nearly - 4 milliseconds more. shouldn't this be a #define? + stars[N].measure.dt = 53.907456; // is this 2048*clockRate ? XXXJester: Nearly - 4 milliseconds more. shouldn't this be a #define? SetSDSSFlags (&stars[N], flags[N], flags2[N]); @@ -199,5 +213,5 @@ altaz (&alt, &az, sidtime - stars[N].average.R, stars[N].average.D, Latitude); - // XXXYYY airmass is in tsField file + // XXXJester: airmass is in tsField file stars[N].measure.airmass = 1.0 / dCOS(90.0 - alt); stars[N].measure.az = az; @@ -206,6 +220,4 @@ } - // Throw away unneeded trailing (I hope) part of Stars - REALLOCATE (stars, Stars, NFILTER*Nstars_used); for (i = 0; i < NFILTER; i++) { @@ -215,5 +227,5 @@ // XXX for now, we define a totally fake coordinate system centered on the first listed star - // XXXYYY what's happening here? *Exact* astrometry for SDSS can only be retrieved using the + // XXXJester: what's happening here? *Exact* astrometry for SDSS can only be retrieved using the // colour-dependent distortion coefficients from the tsField file. See // www.sdss.org/dr7/products/general/astrometry.html @@ -242,5 +254,6 @@ images[N].NX = 2048; - // XXXYYY NZ should be 1361, since the last 128 rows are the duplicate of the neighbouring frame, and only + images[N].NY = 1490; + // XXXJester: NY should be 1361, since the last 128 rows are the duplicate of the neighbouring frame, and only // objects in 64 <= row < 1361+64 (or maybe the <= is at the upper limit?) are counted as belonging to // "this" field. I.e. we also need to read the status flag and use only primary and secondary objects @@ -261,8 +274,7 @@ // run. This is desired in the case of dvo. Hence we need exactly those objects with OK_RUN set, not more, // not less. - - // And once more, NY = 1361, which led into the whole 'status' discussion. In any case, why is this 1490 - // not 1489? - images[N].NY = 1490; + // + // And once more, I think we need NY = 1361, which led into the whole 'status' discussion. In any case, + // why is this 1490 not 1489? images[N].tzero = tzero[i]; @@ -270,5 +282,7 @@ // set photcodes for the 5 images (SDSS_U,G,R,I,Z) - // XXXYYY Why these? and not the 'u_sdss etc. 'ref' photcodes (see above where it says U_SDSS etc.) + // + // We use the 'dep' photcodes and not 'ref' because we can in principle recalibrate the SDSS photometry in + // dvo; looks like 'ref' is reserved for real standard-star observations images[N].photcode = photcode[i]; @@ -343,5 +357,4 @@ 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