Changeset 23941
- Timestamp:
- Apr 21, 2009, 1:34:06 AM (17 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/sj_branches/sj_SDSSaddstar_branch_r23928_20090420/Ohana/src/addstar/src/ReadStarsSDSS.c
r21508 r23941 29 29 // XXX NOTE : as of 2008.02.27, the zero point is still carried internally in millimags 30 30 31 // given a file with the pointer at the start of the table block and the 31 // given a file with the pointer at the start of the table block and the 32 32 // corresponding image header, load the stars from the table 33 33 Stars *ReadStarsSDSS (FILE *f, char *name, Header *header, Header *in_theader, Image *images, int *nimages, unsigned int *nstars) { … … 44 44 int photcode[5]; 45 45 int Nrow, Ncol; // used in the GET_COLUMN_1,5 macros above 46 46 47 47 if (in_theader == NULL) { 48 48 table.header = &theader; … … 51 51 table.header = in_theader; 52 52 Nskip = in_theader[0].size; 53 fseek (f, Nskip, SEEK_CUR); 53 fseek (f, Nskip, SEEK_CUR); 54 54 } 55 55 … … 66 66 strcpy (filtname[4], "z"); 67 67 68 // XXXYYY Why these? and not the 'ref' version? see 'Why these?' below, too. 68 69 NAMED_PHOTCODE_AND_ZP (photcode[0], zeropt[0], "U_SDSS"); 69 70 NAMED_PHOTCODE_AND_ZP (photcode[1], zeropt[1], "G_SDSS"); … … 92 93 tzero[4] = ohana_mjd_to_sec (mjd[4]); 93 94 95 // XXXYYY EDR paper says "psfWidth is is the effective width also determined at the center of each frame. It 96 // is a good generic number to quote for the seeing on each frame." It is stored in the tsField file as 97 // psf_width. 94 98 gfits_scan (table.header, "SEEING_U", "%f", 1, &seeing[0]); 95 99 gfits_scan (table.header, "SEEING_G", "%f", 1, &seeing[1]); … … 98 102 gfits_scan (table.header, "SEEING_Z", "%f", 1, &seeing[4]); 99 103 104 // XXXYYY What's the meaning of photErr? Where is it used again? 100 105 gfits_scan (table.header, "PSFERR_U", "%f", 1, &photErr[0]); 101 106 gfits_scan (table.header, "PSFERR_G", "%f", 1, &photErr[1]); … … 106 111 gfits_scan (header, "CAMCOL", "%d", 1, &camcol); // value in header is usec / unbinned row 107 112 113 // XXXYYY which number is read here? In other words, where is the zeropoint set that is returned by 114 // GetZeroPoint()? I can't find the corresponding SetZeroPoint() 108 115 ZeroPt = GetZeroPoint(); 109 116 … … 132 139 GET_COLUMN_5 (psfCountsErr, float); 133 140 141 GET_COLUMN_1 (status, int); 142 134 143 // the value of stars[].M is supposed to be the instrumental magnitude offset by the 135 144 // default zero point 25.0 (-2.5*log_10(counts/sec) + ZeroPt). The magnitude reported … … 141 150 N = NFILTER*i + j; 142 151 InitStar (&stars[N]); 143 152 144 153 // any values not explicitly set are left at 0.0 145 154 stars[N].average.R = ra[i] + dCOS(dec[i]) * offsetRa[N] / 3600.0; … … 152 161 stars[N].measure.dXccd = ShortPixels(colcErr[N]); 153 162 stars[N].measure.dYccd = ShortPixels(rowcErr[N]); 163 // XXXYYY again the question what ZeroPt-zeropt[j] evalutes to - according to what it says above, this 164 // should take .M back into instrumental magnitudes, but does it? In particular, the *actual* zeropoint 165 // that was used for the computation of this calibrated magnitude is in the tsField file; in general, it 166 // will differ by some small amount (less than 1% I'm guessing) from the 'default' zero points; for 167 // non-photometric stripe82 data, the actual zero point will differ A LOT from the default one. 154 168 stars[N].measure.M = psfCounts[N] + ZeroPt - zeropt[j]; 155 169 stars[N].measure.dM = psfCountsErr[N]; 156 170 stars[N].measure.Map = fiberCounts[N] + ZeroPt - zeropt[j]; 157 stars[N].measure.Sky = sky[N]; // adjust this to counts? 171 stars[N].measure.Sky = sky[N]; // adjust this to counts? XXXYYY It's given in maggies; to get back to counts, 158 172 stars[N].measure.dSky = skyErr[N]; 159 stars[N].measure.FWx = ShortPixels(seeing[j]); // reported in arcsec? 173 stars[N].measure.FWx = ShortPixels(seeing[j]); // reported in arcsec? XXXYYY: Yes! 160 174 stars[N].measure.FWy = ShortPixels(seeing[j]); 161 175 stars[N].measure.psfChisq = prob_psf[N]; // XXX not really the correct value... 162 176 stars[N].measure.detID = N; 163 177 stars[N].measure.t = tzero[j] + clockRate*rowc[N]; // time since row 0 164 stars[N].measure.dt = 53.907456; // is this 2048*clockRate ? 178 stars[N].measure.dt = 53.907456; // is this 2048*clockRate ? XXXYYY Nearly - 4 milliseconds more. shouldn't this be a #define? 165 179 166 180 SetSDSSFlags (&stars[N], flags[N], flags2[N]); … … 178 192 altaz (&alt, &az, sidtime - stars[N].average.R, stars[N].average.D, Latitude); 179 193 194 // XXXYYY airmass is in tsField file 180 195 stars[N].measure.airmass = 1.0 / dCOS(90.0 - alt); 181 196 stars[N].measure.az = az; 182 197 stars[N].measure.photcode = photcode[j]; 183 198 } 184 } 199 } 185 200 186 201 for (i = 0; i < NFILTER; i++) { 187 202 188 203 N = i + *nimages; 189 204 190 205 // XXX for now, we define a totally fake coordinate system centered on the first listed star 206 207 // XXXYYY what's happening here? *Exact* astrometry for SDSS can only be retrieved using the 208 // colour-dependent distortion coefficients from the tsField file. See 209 // www.sdss.org/dr7/products/general/astrometry.html 210 // 211 // WCS headers are of course present in the fpC image, 212 // but then, that would be downloading a lot more data. 191 213 strcpy (images[N].coords.ctype, "RA---TAN"); 192 214 193 215 images[N].coords.crval1 = stars[0].average.R; 194 216 images[N].coords.crval2 = stars[0].average.D; 195 217 196 218 images[N].coords.crpix1 = stars[0].measure.Xccd; 197 219 images[N].coords.crpix2 = stars[0].measure.Yccd; … … 207 229 208 230 images[N].NX = 2048; 231 // XXXYYY NZ should be 1361, since the last 128 rows are the duplicate of the neighbouring frame, and only 232 // objects in 64 <= row < 1361+64 (or maybe the <= is at the upper limit?) are counted as belonging to 233 // "this" field. I.e. we also need to read the status flag and use only primary and secondary objects 234 // (i.e. those that have OK_RUN set in status); otherwise, the SAME PHOTONS will appear more than once in 235 // the dvo database. 236 237 // From http://www.sdss.org/DR7/dm/flatFiles/tsObj.html#status 238 239 // AR_OBJECT_STATUS_OK_RUN = 0x10, /* Located within the primary range */ 240 /* of rows for this field. This is */ 241 /* usable object. This flag is set */ 242 /* by "setObjectStatus". */ 243 // OK_RUN is the union of the following two classes: 244 // AR_OBJECT_STATUS_SECONDARY = 0x1000, /* This is a secondary survey object.*/ 245 // 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 246 // objects rejects repeat observations by design (the SDSS philosophy: everything gets observed 247 // once). Adding 'secondary' objects includes repeat observations of the same object from a different 248 // run. This is desired in the case of dvo. Hence we need exactly those objects with OK_RUN set, not more, 249 // not less. 250 251 // And once more, NY = 1361, which led into the whole 'status' discussion. 209 252 images[N].NY = 1490; 210 253 211 254 images[N].tzero = tzero[i]; 212 255 images[N].cerror = 0.0; 213 256 214 257 // set photcodes for the 5 images (SDSS_U,G,R,I,Z) 258 // XXXYYY Why these? and not the 'u_sdss etc. 'ref' photcodes (see above where it says U_SDSS etc.) 215 259 images[N].photcode = photcode[i]; 216 260 217 261 // calculate this from : C_OBS, TRACKING, and NY 218 262 images[N].exptime = 2048*clockRate; 219 263 220 264 images[N].apmifit = 0.0; 221 265 images[N].dapmifit = 0.0; 222 images[N].detection_limit = 0.0; 266 images[N].detection_limit = 0.0; 223 267 images[N].saturation_limit = 0.0; 224 268 images[N].fwhm_x = seeing[i]; … … 249 293 250 294 images[N].nstar = Nstars; 251 295 252 296 images[N].imageID = N; 253 297 images[N].externID = 0; … … 285 329 int SetSDSSFlags (Stars *star, unsigned int flags1, unsigned int flags2) { 286 330 287 // XXX this is wrong, need to roll left to set the correct bit 331 // XXX this is wrong, need to roll left to set the correct bit 288 332 if (flags1 & 0x00000002) star[0].measure.photFlags |= 0x0001; // BRIGHT - 1 1 289 333 if (flags1 & 0x00000004) star[0].measure.photFlags |= 0x0002; // EDGE - 1 2
Note:
See TracChangeset
for help on using the changeset viewer.
