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 23940) +++ /branches/sj_branches/sj_SDSSaddstar_branch_r23928_20090420/Ohana/src/addstar/src/ReadStarsSDSS.c (revision 23941) @@ -29,5 +29,5 @@ // 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 +// 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, int *nimages, unsigned int *nstars) { @@ -44,5 +44,5 @@ int photcode[5]; int Nrow, Ncol; // used in the GET_COLUMN_1,5 macros above - + if (in_theader == NULL) { table.header = &theader; @@ -51,5 +51,5 @@ table.header = in_theader; Nskip = in_theader[0].size; - fseek (f, Nskip, SEEK_CUR); + fseek (f, Nskip, SEEK_CUR); } @@ -66,4 +66,5 @@ strcpy (filtname[4], "z"); + // XXXYYY Why these? and not the 'ref' version? see 'Why these?' below, too. NAMED_PHOTCODE_AND_ZP (photcode[0], zeropt[0], "U_SDSS"); NAMED_PHOTCODE_AND_ZP (photcode[1], zeropt[1], "G_SDSS"); @@ -92,4 +93,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 + // is a good generic number to quote for the seeing on each frame." It is stored in the tsField file as + // psf_width. gfits_scan (table.header, "SEEING_U", "%f", 1, &seeing[0]); gfits_scan (table.header, "SEEING_G", "%f", 1, &seeing[1]); @@ -98,4 +102,5 @@ gfits_scan (table.header, "SEEING_Z", "%f", 1, &seeing[4]); + // XXXYYY 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]); @@ -106,4 +111,6 @@ 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() ZeroPt = GetZeroPoint(); @@ -132,4 +139,6 @@ GET_COLUMN_5 (psfCountsErr, float); + GET_COLUMN_1 (status, int); + // 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 @@ -141,5 +150,5 @@ 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; @@ -152,15 +161,20 @@ 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 + // 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 + // will differ by some small amount (less than 1% I'm guessing) from the 'default' zero points; for + // non-photometric stripe82 data, the actual zero point will differ A LOT from the default one. 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.Sky = sky[N]; // adjust this to counts? XXXYYY It's given in maggies; to get back to counts, stars[N].measure.dSky = skyErr[N]; - stars[N].measure.FWx = ShortPixels(seeing[j]); // reported in arcsec? + stars[N].measure.FWx = ShortPixels(seeing[j]); // reported in arcsec? XXXYYY: 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 ? + stars[N].measure.dt = 53.907456; // is this 2048*clockRate ? XXXYYY Nearly - 4 milliseconds more. shouldn't this be a #define? SetSDSSFlags (&stars[N], flags[N], flags2[N]); @@ -178,20 +192,28 @@ altaz (&alt, &az, sidtime - stars[N].average.R, stars[N].average.D, Latitude); + // XXXYYY airmass is in tsField file stars[N].measure.airmass = 1.0 / dCOS(90.0 - alt); stars[N].measure.az = az; stars[N].measure.photcode = photcode[j]; } - } + } 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 + + // XXXYYY 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 + // + // WCS headers are of course present in the fpC image, + // but then, that would be downloading a lot more data. strcpy (images[N].coords.ctype, "RA---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; @@ -207,18 +229,40 @@ images[N].NX = 2048; + // XXXYYY NZ 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 + // (i.e. those that have OK_RUN set in status); otherwise, the SAME PHOTONS will appear more than once in + // the dvo database. + + // From http://www.sdss.org/DR7/dm/flatFiles/tsObj.html#status + + // AR_OBJECT_STATUS_OK_RUN = 0x10, /* Located within the primary range */ + /* of rows for this field. This is */ + /* usable object. This flag is set */ + /* by "setObjectStatus". */ + // OK_RUN is the union of the following two classes: + // AR_OBJECT_STATUS_SECONDARY = 0x1000, /* This is a secondary survey object.*/ + // AR_OBJECT_STATUS_PRIMARY = 0x2000, /* This is a primary survey object. */ // 'primary' means 'inside the area that defines a unique coverage of the sky', i.e. selecting primary + // objects rejects repeat observations by design (the SDSS philosophy: everything gets observed + // once). Adding 'secondary' objects includes repeat observations of the same object from a different + // 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. 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) + // XXXYYY Why these? and not the 'u_sdss etc. 'ref' photcodes (see above where it says U_SDSS etc.) 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].detection_limit = 0.0; images[N].saturation_limit = 0.0; images[N].fwhm_x = seeing[i]; @@ -249,5 +293,5 @@ images[N].nstar = Nstars; - + images[N].imageID = N; images[N].externID = 0; @@ -285,5 +329,5 @@ int SetSDSSFlags (Stars *star, unsigned int flags1, unsigned int flags2) { - // XXX this is wrong, need to roll left to set the correct bit + // 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