001 package util.sdss;
002
003 import static java.lang.Math.*;
004 import java.io.*;
005 import java.text.DecimalFormat;
006 import java.util.*;
007 import org.eso.fits.*;
008
009 /**
010 * A data container to hold the quasar instances from
011 * <a href="http://www.sdss.org">Sloan Digital Sky Survey</a>
012 * DR3 QSO Catalog
013 * that appear in the file dr1qso.txt along with emission lines locations.
014 *
015 * <h3>Acknowledgment</h3>
016 * Funding for the creation and distribution of the SDSS Archive has been provided
017 * by the Alfred P. Sloan Foundation, the Participating Institutions, the National
018 * Aeronautics and Space Administration, the National Science Foundation, the U.S.
019 * Department of Energy, the Japanese Monbukagakusho, and the Max Planck Society.
020 * The SDSS Web site is http://www.sdss.org/.
021 *
022 * The SDSS is managed by the Astrophysical Research Consortium (ARC) for the
023 * Participating Institutions. The Participating Institutions are The University of
024 * Chicago, Fermilab, the Institute for Advanced Study, the Japan Participation
025 * Group, The Johns Hopkins University, the Korean Scientist Group, Los Alamos
026 * National Laboratory, the Max-Planck-Institute for Astronomy (MPIA), the Max-
027 * Planck-Institute for Astrophysics (MPA), New Mexico State University, University
028 * of Pittsburgh, Princeton University, the United States Naval Observatory, and
029 * the University of Washington.
030 */
031 public class Quasar implements Star, Comparable<Quasar> {
032
033 public static DecimalFormat zFormatter = new DecimalFormat("0.00");
034
035 public String name; // columns 1-19
036 public String ra; // columns 20-30
037 public String dec; // columns 31-41
038 public String zem; // columns 42-48
039 public String umag; // columns 51-57
040 public String umagErr; // columns 58-63
041 public String gmag; // columns 64-70
042 public String gmagErr; // columns 71-76
043 public String rmag; // columns 77-83
044 public String rmagErr; // columns 84-89
045 public String imag; // columns 90-96
046 public String imagErr; // columns 97-102
047 public String zmag; // columns 103-109
048 public String zmagErr; // columns 110-115
049 public String au; // columns 116-122
050 public String first; // columns 123-129
051 public String firstSN; // columns 130-137
052 public String firstSep;// columns 138-144
053 public String rosat; // columns 145-152
054 public String rosatSN; // columns 153-159
055 public String rosatSep;// columns 160-166
056 public String j2mass; // columns 167-173
057 public String j2massErr;//columns 174-179
058 public String h2mass; // columns 180-186
059 public String h2massErr;//columns 187-192
060 public String k2mass; // columns 193-186
061 public String k2massErr;//columns 200-205
062 public int s_mjd; // columns 281-286
063 public int plate; // columns 287-291
064 public int fiber; // columns 292-296
065 public String altName; // columns 297-322
066
067 protected List<util.sdss.Line> emissionLines;
068 protected List<util.sdss.Line> absorptionLines;
069
070 public Quasar(double z) {
071 zem = "" + z;
072 }
073
074 public static Integer getPrimaryKey(int aPlate, int aFiber) {
075 return new Integer(aPlate * 1000 + aFiber);
076 }
077
078 public Integer getPrimaryKey() {
079 return getPrimaryKey(plate, fiber);
080 }
081
082 public double getRa() {
083 return Double.parseDouble(ra);
084 }
085
086 public double getDec() {
087 return Double.parseDouble(dec);
088 }
089
090 public double getZem() {
091 return Double.parseDouble(zem);
092 }
093
094 public double[] getMagnitude() {
095 return new double[] {
096 Double.parseDouble(rosat),
097 Double.parseDouble(umag),
098 Double.parseDouble(gmag),
099 Double.parseDouble(rmag),
100 Double.parseDouble(imag),
101 Double.parseDouble(zmag),
102 Double.parseDouble(j2mass),
103 Double.parseDouble(h2mass),
104 Double.parseDouble(k2mass),
105 Double.parseDouble(first)};
106 }
107
108 public double[] getMagnitudeError(int aBand) {
109 return new double[] {
110 Double.parseDouble(rosatSN),
111 Double.parseDouble(umagErr),
112 Double.parseDouble(gmagErr),
113 Double.parseDouble(rmagErr),
114 Double.parseDouble(imagErr),
115 Double.parseDouble(zmagErr),
116 Double.parseDouble(j2massErr),
117 Double.parseDouble(h2massErr),
118 Double.parseDouble(k2massErr),
119 Double.parseDouble(firstSN)};
120 }
121
122 public Quasar(String aRow) {
123 // columns are the 1 offset position up to 1-plus end position
124 name= aRow.substring(0, 18); // columns 1-19
125
126 ra = aRow.substring(19, 29); // columns 20-30
127 dec = aRow.substring(30, 40); // columns 31-41
128 zem = aRow.substring(41, 47); // columns 42-48
129
130 umag = aRow.substring(50, 56); // columns 51-57
131 umagErr= aRow.substring(57, 62); // columns 58-63
132
133 gmag = aRow.substring(63, 69); // columns 64-70
134 gmagErr= aRow.substring(70, 75); // columns 71-76
135
136 rmag = aRow.substring(76, 82); // columns 77-83
137 rmagErr= aRow.substring(83, 88); // columns 84-89
138
139 imag = aRow.substring(89, 95); // columns 90-96
140 imagErr= aRow.substring(96, 101); // columns 97-102
141
142 zmag = aRow.substring(102, 108); // columns 103-109
143 zmagErr= aRow.substring(109, 114); // columns 110-115
144
145 au = aRow.substring(115, 121); // columns 116-122
146
147 first = aRow.substring(122, 128); // columns 123-129
148 firstSN = aRow.substring(129, 136); // columns 130-137
149 firstSep= aRow.substring(137, 143); // columns 138-144
150
151 rosat = aRow.substring(144, 151); // columns 145-152
152 rosatSN = aRow.substring(152, 158); // columns 153-159
153 rosatSep= aRow.substring(159, 165); // columns 160-166
154
155 j2mass = aRow.substring(166, 172); // columns 167-173
156 j2massErr= aRow.substring(173, 178); // columns 174-179
157
158 h2mass = aRow.substring(179, 185); // columns 180-186
159 h2massErr= aRow.substring(186, 191); // columns 187-192
160
161 k2mass = aRow.substring(192, 198); // columns 193-186
162 k2massErr= aRow.substring(199, 204); // columns 200-205
163
164 s_mjd = Integer.parseInt(aRow.substring(280, 285).trim()); // columns 281-286
165 plate = Integer.parseInt(aRow.substring(286, 290).trim()); // columns 287-291
166 fiber = Integer.parseInt(aRow.substring(291, 295).trim()); // columns 292-296
167 altName = aRow.substring(296, 321); // columns 297-322
168
169 }
170
171
172 public int compareTo(Quasar aQuasar) {
173 if (aQuasar == null) {
174 System.out.println("null pointer for quasar");
175 throw(new NullPointerException("Compare to null quasar instance"));
176 } else {
177 double w = aQuasar.getZem();
178 if (getZem() > w) return 1;
179 else if (getZem() < w) return -1;
180 }
181 return 0;
182 }
183
184 /** Get the spectra file for this quasar */
185 public File getSpectraFile() {
186 return getSpectraFile(plate, fiber, s_mjd);
187 }
188
189 /** Get the spectra file for the given plate fiber and julian date */
190 public static File getSpectraFile(int aPlate, int aFiber, int aJulianDate) {
191 StringBuffer buffer = new StringBuffer(50);
192 buffer.append("H:/sdss/0"); //M:/sdss2/0
193 buffer.append(aPlate);
194 buffer.append("/spSpec/");
195 buffer.append("spSpec-");
196 buffer.append(aJulianDate);
197 buffer.append("-0");
198 buffer.append(aPlate);
199 buffer.append("-");
200 if (aFiber < 10) buffer.append("00");
201 else if (aFiber < 100) buffer.append("0");
202 buffer.append(aFiber);
203 buffer.append(".fit");
204 return new File(buffer.toString());
205 }
206
207 // also do for ps.gz
208 public static String getSpectraImageFileName(int aPlate, int aFiber, int aJulianDate) {
209 StringBuffer buffer = new StringBuffer(50);
210 buffer.append("0");
211 buffer.append(aPlate);
212 buffer.append("/gif/");
213 buffer.append("spPlot-");
214 buffer.append(aJulianDate);
215 buffer.append("-0");
216 buffer.append(aPlate);
217 buffer.append("-");
218 if (aFiber < 10)
219 buffer.append("00");
220 else if (aFiber < 100)
221 buffer.append("0");
222 buffer.append(aFiber);
223 buffer.append(".ps.gz"); // or .gif
224 return buffer.toString();
225 }
226
227 // http://das.sdss.org/DR1/data/spectro/1d_20/0266/gif/
228 public static java.net.URL getSpectraImageURL(int aPlate, int aFiber, int aJulianDate) {
229 StringBuffer buffer = new StringBuffer(50);
230 buffer.append("http://das.sdss.org/DR1/data/spectro/1d_20/");
231 buffer.append(getSpectraImageFileName(aPlate, aFiber, aJulianDate));
232 java.net.URL gif = null;
233 try {
234 gif = new java.net.URL(buffer.toString());
235 } catch (java.net.MalformedURLException e) {
236 System.out.println("Malformed URL " + buffer.toString());
237 }
238 return gif;
239 }
240
241 /** Get the image file for this quasar */
242 public File getImageFile() {
243 return getImageFile(plate, fiber, s_mjd);
244 }
245
246 /** Get the image file for the given plate fiber and julian date */
247 public static File getImageFile(int aPlate, int aFiber, int aJulianDate) {
248 StringBuffer buffer = new StringBuffer(50);
249 buffer.append("R:/SDSS/0");
250 buffer.append(aPlate);
251 buffer.append("/spAtlas/");
252 buffer.append("spAtlas-0");
253 buffer.append(aPlate);
254 buffer.append("-");
255 buffer.append(aJulianDate);
256 buffer.append("-");
257 buffer.append(aFiber);
258 buffer.append(".fit");
259 return new File(buffer.toString());
260 }
261
262 public static int entries = 16713; //5349;
263
264 /** Get a list of quasars from SDSS
265 * @return a list of quasars
266 */
267 public static List<util.sdss.Quasar> readQuasars() throws java.io.IOException {
268 FileReader reader = new FileReader("data/sdss/dr1qso-by-plate.txt");
269 BufferedReader in = new BufferedReader(reader);
270 int count = 1;
271 List<util.sdss.Quasar> quasars = new ArrayList<util.sdss.Quasar>(entries);
272 while (true) {
273 String row = in.readLine();
274 if (count > entries) break; // last line
275 Quasar qso = new Quasar(row);
276 quasars.add(qso);
277 //System.out.println(qso + " file=" + qso.getSpectraFile());
278 count++;
279 }
280 reader.close();
281 return quasars;
282 }
283
284 public String toString() {
285 String lines = "";
286 String blank = " ";
287 // 13 entries: 3057.13 3057.13 3057.13 3057.13 3057.13 3057.13 3057.13 3057.13 3057.13 3057.13 3057.13 3057.13 3057.13
288 if (emissionLines != null) { // create a concatenated list of emission lines
289 for (Iterator x = getEmissionLines().iterator(); x.hasNext(); ) {
290 Line line = (Line) x.next();
291 String wavelengthString = zFormatter.format(line.getWavelength());
292 if (wavelengthString.length() == 7) wavelengthString = " " + wavelengthString;
293 lines += wavelengthString;
294 }
295 if (lines.length() > blank.length()) {
296 lines = lines.substring(0, blank.length());
297 } else {
298 lines = lines + blank.substring(0, blank.length()-lines.length());
299 }
300 } else {
301 lines = blank;
302 }
303 StringBuffer fiberString = new StringBuffer();
304 fiberString.append(" ");
305 if (fiber < 10)
306 fiberString.append("00");
307 else if (fiber < 100)
308 fiberString.append("0");
309 fiberString.append(fiber);
310
311 return name + " " + zem + lines + " " + plate + " " + fiberString.toString() + " " + altName + " " + ra + " " + dec + " " + umag;
312 }
313
314 public List<util.sdss.Line> getEmissionLines() {
315 // TODO: LAZY evaluation of emission lines from each quasar FITS spectra file
316 if (emissionLines == null) {
317 File spectraFile = getSpectraFile();
318 FitsFile file = null;
319 try {
320 file = new FitsFile(spectraFile);
321 int noHDU = file.getNoHDUnits();
322 boolean gotLines = false;
323 for (int i = 0; i < noHDU && !gotLines; i++) {
324 FitsHDUnit hdu = file.getHDUnit(i);
325 FitsHeader hdr = hdu.getHeader();
326 int type = hdr.getType();
327 if (type==Fits.BTABLE || type==Fits.ATABLE) {
328 gotLines = true; // ensures that only the first binary table is read
329 FitsTable table = (FitsTable) hdu.getData();
330 emissionLines = Line.getLines(this, table);
331 }
332 }
333 if (file != null) file.closeFile();
334 } catch (FitsException e) {
335 System.out.println("Fits exception:" + e);
336 } catch (IOException e) {
337 System.out.println("Error: cannot open file >" + file + "<");
338 }
339 }
340 return emissionLines;
341 }
342
343 public List<util.sdss.Line> getAbsorptionLines() {
344 return null;
345 }
346
347 // lazy evaluation paramaters
348 double coeff0 = 0.0d, coeff1 = 0.0d; // Wavelength parameters
349 int n = 0; // Number of sample points
350
351 /** Get the wavelength parameters n, coeff0 and coeff1 from the FITS spectra file.
352 */
353 protected void getWavelengthParameters() {
354 File spectraFile = getSpectraFile();
355 FitsFile spSpecFile = null;
356 try {
357 spSpecFile = new FitsFile(spectraFile);
358 int noHDU = spSpecFile.getNoHDUnits();
359 for (int i = 0; i < noHDU; i++) {
360 FitsHDUnit hdu = spSpecFile.getHDUnit(i);
361 FitsHeader hdr = hdu.getHeader();
362 int type = hdr.getType();
363 if (type == Fits.IMAGE) {
364 FitsMatrix dm = (FitsMatrix) hdu.getData();
365 FitsKeyword kw = hdr.getKeyword("COEFF0");
366 coeff0 = kw.getReal();
367 kw = hdr.getKeyword("COEFF1");
368 coeff1 = kw.getReal();
369 int naxis[] = dm.getNaxis();
370 if (dm.getNoValues() > 0) n = naxis[0];
371 }
372 }
373 if (spSpecFile != null) spSpecFile.closeFile();
374 } catch (FitsException e) {
375 System.out.println("Fits exception:" + e);
376 } catch (IOException e) {
377 System.out.println("Error: cannot open file >" + spSpecFile + "<");
378 }
379 }
380
381 /** Get the vacuum wavelength corresponding to the fractional sample point,
382 * possibly in between sample points (in Angstroms).
383 * @param samplePoint a floating point number representing the fractional sample point
384 * @return vacuum wavelength of the sample point or zero if out of range
385 */
386 public float getWavelength(float samplePoint) {
387 float wavelength = 0.0f;
388 if (coeff1 == 0.0d) getWavelengthParameters(); // get coeff0 and coeff1 on first invocation
389 if (0 <= samplePoint && samplePoint < n)
390 wavelength = (float) Math.pow(10.0d, coeff0 + ((double) samplePoint)*coeff1);
391
392 return wavelength;
393 }
394
395 /** Get the vacuum wavelength corresponding to the sample point (in Angstroms).
396 * @param samplePoint the sample point
397 * @return vacuum wavelength of the sample point or zero if out of range
398 */
399 public float getWavelength(int samplePoint) {
400 float wavelength = 0.0f;
401 if (coeff1 == 0.0d) getWavelengthParameters(); // get coeff0 and coeff1 on first invocation
402 if (0 <= samplePoint && samplePoint < n)
403 wavelength = (float) Math.pow(10.0d, coeff0 + ((double) samplePoint)*coeff1);
404
405 return wavelength;
406 }
407
408 /** Get the vacuum wavelengths corresponding all the spectra points (in Angstroms).
409 * @return the vacuum wavelengths of all the sample points
410 */
411 public float[] getWavelengths() {
412 float[] wavelength = null;
413 if (coeff1 == 0.0d) getWavelengthParameters(); // get coeff0 and coeff1 on first invocation
414 wavelength = new float[n];
415 for (int j = 0; j < n; j++)
416 wavelength[j] = getWavelength(j);
417
418 return wavelength;
419 }
420
421 /** Convert from vacuum wavelength to air wavelength (in angstroms).
422 * The conversion equation is from the
423 * <a href="http://www.sdss.org/dr3/products/spectra/vacwavelength.html">SDSS DR3</a> web page.
424 * @param vacuumWavelength the vacuum wavelength in angstroms
425 * @return air wavelength of the given vacuum wavelength
426 */
427 public static float convertToAirWavelength(float vacuumWavelength) {
428 float vacuumWavelength2 = vacuumWavelength * vacuumWavelength;
429 float vacuumWavelength4 = vacuumWavelength2* vacuumWavelength2;
430 return vacuumWavelength / (1f + 2.735182E-4f + 131.4182f / vacuumWavelength2 + 2.76249E8f / vacuumWavelength4);
431 }
432
433 /** Convert from air wavelength to vacuum wavelength (in angstroms).
434 * Equation is from the IAU standard for conversion from air to vacuum wavelengths.
435 * (<a href="http://adsabs.harvard.edu/cgi-bin/nph-bib_query?bibcode=1991ApJS..77...119M">Morton: 1991, ApJS, 77, 119</a>).
436 * Pre-condition: air wavelength > 2000 Å
437 * @param airWavelength the air wavelength in angstroms
438 * @return vacuum wavelength of the given air wavelength
439 */
440 public static float convertToVacuumWavelength(float airWavelength) {
441 float sigma = 1e-4f / airWavelength;
442 float sigma2 = sigma * sigma;
443 return airWavelength * (1f + 6.4328E-5f + 2.94981E-2f / (146f - sigma2) + 2.5540E-4f / (41f - sigma2));
444 }
445
446 float[] spectra; // lazy evaluation of spectra
447 public float[] getSpectra() {
448 if (spectra == null)
449 spectra = getArray(0);
450 return spectra;
451 }
452
453 /** Get the spectra intensity at sample point.
454 * @param samplePoint the sample point
455 * @return the spectra intensity at the sample point
456 */
457 public float getIntensity(int samplePoint) {
458 getSpectra();
459 return spectra[samplePoint];
460 }
461
462 /** Get the spectra intensity at the given positive wavelength, linear interpolation is
463 * performed between sample points.
464 * @param wavelength the wavelength at which the intensity is requested
465 * @return the spectra intensity at the wavelength
466 */
467 public float getIntensity(float wavelength) {
468 if (wavelength <= 0.0f) return 0.0f;
469 getSpectra();
470 double jfractional = (log((double)wavelength) / log(10.0d) - coeff0) / coeff1;
471 int j = (int) floor(jfractional);
472 float delta = (float) (jfractional - (double)j);
473 if (0 <= j && j < n - 1)
474 return spectra[j] * (1.0f - delta) + spectra[j + 1] * delta;
475 else
476 return 0.0f;
477 }
478
479 float[] noise; // lazy evaluation of noise
480 public float[] getNoise() {
481 if (noise == null)
482 noise= getArray(2);
483 return noise;
484 }
485
486 public float[] getArray(int row) {
487 if (coeff1 == 0.0d) getWavelengthParameters(); // get coeff0 and coeff1 on first invocation
488 float[] spectra = null;
489 File spectraFile = getSpectraFile();
490 FitsFile spSpecFile = null;
491 try {
492 spSpecFile = new FitsFile(spectraFile);
493 int noHDU = spSpecFile.getNoHDUnits();
494 for (int i = 0; i < noHDU; i++) {
495 FitsHDUnit hdu = spSpecFile.getHDUnit(i);
496 FitsHeader hdr = hdu.getHeader();
497 int type = hdr.getType();
498 if (type == Fits.IMAGE) {
499 FitsMatrix dm = (FitsMatrix) hdu.getData();
500 int nval = dm.getNoValues();
501 if (nval > 0) {
502 int nrow = nval / n;
503 spectra = new float[n]; // mandatory !!
504 try {
505 dm.getFloatValues(n * row, n, spectra);
506 } catch (FitsException e) {
507 spectra = null;
508 }
509 }
510 }
511 }
512 if (spSpecFile != null) spSpecFile.closeFile();
513 } catch (FitsException e) {
514 System.out.println("Fits exception:" + e);
515 } catch (IOException e) {
516 System.out.println("Error: cannot open file >" + spSpecFile + "<");
517 }
518 return spectra;
519 }
520
521 public void showSpectra() {
522 File spectraFile = getSpectraFile();
523 FitsFile spSpecFile = null;
524 try {
525 System.out.println(spectraFile);
526 spSpecFile = new FitsFile(spectraFile);
527 String name = spSpecFile.getName();
528 int noHDU = spSpecFile.getNoHDUnits();
529 for (int i = 0; i < noHDU; i++) {
530 FitsHDUnit hdu = spSpecFile.getHDUnit(i);
531 FitsHeader hdr = hdu.getHeader();
532 int noKw = hdr.getNoKeywords();
533 int type = hdr.getType();
534 int size = (int) hdr.getDataSize();
535 if (type == Fits.IMAGE) {
536 FitsMatrix dm = (FitsMatrix) hdu.getData();
537 DecimalFormat xFormatter = new DecimalFormat("0.000");
538 DecimalFormat yFormatter = new DecimalFormat("0.000");
539 int naxis[] = dm.getNaxis();
540 double crval[] = dm.getCrval();
541 double crpix[] = dm.getCrpix();
542 double cdelt[] = dm.getCdelt();
543 FitsKeyword kw = hdr.getKeyword("COEFF0");
544 double coeff0 = kw.getReal();
545 kw = hdr.getKeyword("COEFF1");
546 double coeff1 = kw.getReal();
547 int nv, off, npix;
548 int nval = dm.getNoValues();
549 if (0<nval) {
550 int ncol = naxis[0];
551 int nrow = nval/ncol;
552 float data[] = new float[ncol];
553 double val;
554 off = nv = npix = 0;
555 //nrow = 1; // Only get first row
556 for (int nr=0; nr<nrow; nr++) {
557 try {
558 dm.getFloatValues(off, ncol, data);
559 StringBuffer s;
560 for (int n = 0; n<ncol; n++) {
561 double wavelength = Math.pow(10.0, coeff0 + ((double) n)*coeff1);
562 val = data[n];
563 s = new StringBuffer(xFormatter.format(wavelength));
564 s.append(" ");
565 String y = yFormatter.format(val);
566 s.append(y);
567 s.append("\n");
568 System.out.print(s.toString());
569 npix++;
570 }
571 } catch (FitsException e) {
572 }
573 off += ncol;
574 }
575 }
576 }
577 }
578 if (spSpecFile != null) spSpecFile.closeFile();
579 } catch (FitsException e) {
580 System.out.println("Fits exception:" + e);
581 } catch (IOException e) {
582 System.out.println("Error: cannot open file >" + spSpecFile + "<");
583 }
584 }
585
586 }
587
588