1 """module for performing calculations on Spectral Energy Distributions (SEDs)
2
3 (c) 2007-2013 Matt Hilton
4
5 U{http://astlib.sourceforge.net}
6
7 This module provides classes for manipulating SEDs, in particular the Bruzual &
8 Charlot 2003, Maraston 2005, and Percival et al 2009 stellar population
9 synthesis models are currently supported. Functions are provided for
10 calculating the evolution of colours and magnitudes in these models with
11 redshift etc., and for fitting broadband photometry using these models.
12
13 @var VEGA: The SED of Vega, used for calculation of magnitudes on the Vega system.
14 @type VEGA: L{SED} object
15 @var AB: Flat spectrum SED, used for calculation of magnitudes on the AB system.
16 @type AB: L{SED} object
17 @var SOL: The SED of the Sun.
18 @type SOL: L{SED} object
19
20 """
21
22
23 import numpy
24 import math
25 import operator
26 try:
27 from scipy import interpolate
28 from scipy import ndimage
29 except ImportError:
30 print("WARNING: astSED: failed to import scipy modules - some functions "
31 "will not work.")
32 import astLib
33 from astLib import astCalc
34 import os
35 try:
36 import matplotlib
37 from matplotlib import pylab
38 matplotlib.interactive(False)
39 except ImportError:
40 print("WARNING: astSED: failed to import matplotlib - some functions will "
41 "not work.")
42 import glob
43
44
46 """This class describes a filter transmission curve. Passband objects are
47 created by loading data from from text files containing wavelength in
48 angstroms in the first column, relative transmission efficiency in the
49 second column (whitespace delimited). For example, to create a Passband
50 object for the 2MASS J filter:
51
52 passband=astSED.Passband("J_2MASS.res")
53
54 where "J_2MASS.res" is a file in the current working directory that
55 describes the filter.
56
57 Wavelength units can be specified as 'angstroms', 'nanometres' or
58 'microns'; if either of the latter, they will be converted to angstroms.
59
60 """
61
62 - def __init__(self, fileName, normalise=True, inputUnits='angstroms'):
63
64 inFile = open(fileName, "r")
65 lines = inFile.readlines()
66
67 wavelength = []
68 transmission = []
69 for line in lines:
70
71 if line[0] != "#" and len(line) > 3:
72
73 bits = line.split()
74 transmission.append(float(bits[1]))
75 wavelength.append(float(bits[0]))
76
77 self.wavelength = numpy.array(wavelength)
78 self.transmission = numpy.array(transmission)
79
80 if inputUnits == 'angstroms':
81 pass
82 elif inputUnits == 'nanometres':
83 self.wavelength = self.wavelength * 10.0
84 elif inputUnits == 'microns':
85 self.wavelength = self.wavelength * 10000.0
86 elif inputUnits == 'mm':
87 self.wavelength = self.wavelength * 1e7
88 elif inputUnits == 'GHz':
89 self.wavelength = 3e8 / (self.wavelength * 1e9)
90 self.wavelength = self.wavelength * 1e10
91 else:
92 raise Exception("didn't understand passband input units")
93
94
95 merged = numpy.array([self.wavelength, self.transmission]).transpose()
96 sortedMerged = numpy.array(sorted(merged, key=operator.itemgetter(0)))
97 self.wavelength = sortedMerged[:, 0]
98 self.transmission = sortedMerged[:, 1]
99
100 if normalise:
101 self.transmission = self.transmission / numpy.trapz(
102 self.transmission, self.wavelength)
103
104
105 self.interpolator = interpolate.interp1d(self.wavelength,
106 self.transmission,
107 kind='linear')
108
110 """Returns a two dimensional list of [wavelength, transmission],
111 suitable for plotting by gnuplot.
112
113 @rtype: list
114 @return: list in format [wavelength, transmission]
115
116 """
117
118 listData = []
119 for l, f in zip(self.wavelength, self.transmission):
120 listData.append([l, f])
121
122 return listData
123
124 - def rescale(self, maxTransmission):
125 """Rescales the passband so that maximum value of the transmission is
126 equal to maxTransmission. Useful for plotting.
127
128 @type maxTransmission: float
129 @param maxTransmission: maximum value of rescaled transmission curve
130
131 """
132
133 self.transmission = self.transmission * (maxTransmission /
134 self.transmission.max())
135
136 - def plot(self, xmin='min', xmax='max', maxTransmission=None):
137 """Plots the passband, rescaling the maximum of the tranmission curve
138 to maxTransmission if required.
139
140 @type xmin: float or 'min'
141 @param xmin: minimum of the wavelength range of the plot
142 @type xmax: float or 'max'
143 @param xmax: maximum of the wavelength range of the plot
144 @type maxTransmission: float
145 @param maxTransmission: maximum value of rescaled transmission curve
146
147 """
148
149 if maxTransmission is not None:
150 self.rescale(maxTransmission)
151
152 pylab.matplotlib.interactive(True)
153 pylab.plot(self.wavelength, self.transmission)
154
155 if xmin == 'min':
156 xmin = self.wavelength.min()
157 if xmax == 'max':
158 xmax = self.wavelength.max()
159
160 pylab.xlim(xmin, xmax)
161 pylab.xlabel("Wavelength")
162 pylab.ylabel("Relative Flux")
163
165 """Calculates effective wavelength for the passband. This is the same
166 as equation (3) of Carter et al. 2009.
167
168 @rtype: float
169 @return: effective wavelength of the passband, in Angstroms
170
171 """
172
173 a = numpy.trapz(self.transmission * self.wavelength, self.wavelength)
174 b = numpy.trapz(self.transmission / self.wavelength, self.wavelength)
175 effWavelength = numpy.sqrt(a / b)
176
177 return effWavelength
178
179
180
182 """This class generates a passband with a top hat response between the
183 given wavelengths.
184
185 """
186
187 - def __init__(self, wavelengthMin, wavelengthMax, normalise=True):
188 """Generates a passband object with top hat response between
189 wavelengthMin, wavelengthMax. Units are assumed to be Angstroms.
190
191 @type wavelengthMin: float
192 @param wavelengthMin: minimum of the wavelength range of the passband
193 @type wavelengthMax: float
194 @param wavelengthMax: maximum of the wavelength range of the passband
195 @type normalise: bool
196 @param normalise: if True, scale such that total area under the
197 passband over the wavelength
198 range is 1.
199
200 """
201
202 self.wavelength = numpy.arange(
203 wavelengthMin, wavelengthMax + 10,
204 10, dtype=float)
205 self.transmission = numpy.ones(self.wavelength.shape, dtype=float)
206
207 if normalise:
208 self.transmission = self.transmission / numpy.trapz(
209 self.transmission, self.wavelength)
210
211
212 self.interpolator = interpolate.interp1d(self.wavelength,
213 self.transmission,
214 kind='linear')
215
216
217
219 """This class describes a Spectral Energy Distribution (SED).
220
221 To create a SED object, lists (or numpy arrays) of wavelength and relative
222 flux must be provided. The SED can optionally be redshifted. The wavelength
223 units of SEDs are assumed to be Angstroms - flux calculations using
224 Passband and SED objects specified with different wavelength units will be
225 incorrect.
226
227 The L{StellarPopulation} class (and derivatives) can be used to extract
228 SEDs for specified ages from e.g. the Bruzual & Charlot 2003 or Maraston
229 2005 models.
230
231 """
232
233 - def __init__(self,
234 wavelength=[],
235 flux=[],
236 z=0.0,
237 ageGyr=None,
238 normalise=False,
239 label=None):
240
241
242
243
244 self.z0wavelength = numpy.array(wavelength)
245 self.z0flux = numpy.array(flux)
246 self.wavelength = numpy.array(wavelength)
247 self.flux = numpy.array(flux)
248 self.z = z
249 self.label = label
250
251
252
253 self.EBMinusV = 0.0
254 self.intrinsic_z0flux = numpy.array(flux)
255
256 if normalise:
257 self.normalise()
258
259 if z != 0.0:
260 self.redshift(z)
261
262 self.ageGyr = ageGyr
263
265 """Copies the SED, returning a new SED object
266
267 @rtype: L{SED} object
268 @return: SED
269
270 """
271
272 newSED = SED(wavelength=self.z0wavelength,
273 flux=self.z0flux,
274 z=self.z,
275 ageGyr=self.ageGyr,
276 normalise=False,
277 label=self.label)
278
279 return newSED
280
282 """Loads SED from a white space delimited file in the format
283 wavelength, flux. Lines beginning with # are ignored.
284
285 @type fileName: string
286 @param fileName: path to file containing wavelength, flux data
287
288 """
289
290 inFile = open(fileName, "r")
291 lines = inFile.readlines()
292 inFile.close()
293 wavelength = []
294 flux = []
295 for line in lines:
296 if line[0] != "#" and len(line) > 3:
297 bits = line.split()
298 wavelength.append(float(bits[0]))
299 flux.append(float(bits[1]))
300
301
302 if wavelength[0] > wavelength[-1]:
303 wavelength.reverse()
304 flux.reverse()
305
306 self.z0wavelength = numpy.array(wavelength)
307 self.z0flux = numpy.array(flux)
308 self.wavelength = numpy.array(wavelength)
309 self.flux = numpy.array(flux)
310
312 """Writes SED to a white space delimited file in the format wavelength,
313 flux.
314
315 @type fileName: string
316 @param fileName: path to file
317
318 """
319
320 outFile = open(fileName, "w")
321 for l, f in zip(self.wavelength, self.flux):
322 outFile.write(str(l) + " " + str(f) + "\n")
323 outFile.close()
324
326 """Returns a two dimensional list of [wavelength, flux], suitable for
327 plotting by gnuplot.
328
329 @rtype: list
330 @return: list in format [wavelength, flux]
331
332 """
333
334 listData = []
335 for l, f in zip(self.wavelength, self.flux):
336 listData.append([l, f])
337
338 return listData
339
340 - def plot(self, xmin='min', xmax='max'):
341 """Produces a simple (wavelength, flux) plot of the SED.
342
343 @type xmin: float or 'min'
344 @param xmin: minimum of the wavelength range of the plot
345 @type xmax: float or 'max'
346 @param xmax: maximum of the wavelength range of the plot
347
348 """
349
350 pylab.matplotlib.interactive(True)
351 pylab.plot(self.wavelength, self.flux)
352
353 if xmin == 'min':
354 xmin = self.wavelength.min()
355 if xmax == 'max':
356 xmax = self.wavelength.max()
357
358
359 plotMask = numpy.logical_and(
360 numpy.greater(self.wavelength, xmin), numpy.less(self.wavelength,
361 xmax))
362 plotMax = self.flux[plotMask].max()
363 pylab.ylim(0, plotMax * 1.1)
364 pylab.xlim(xmin, xmax)
365 pylab.xlabel("Wavelength")
366 pylab.ylabel("Relative Flux")
367
368 - def integrate(self, wavelengthMin='min', wavelengthMax='max'):
369 """Calculates flux in SED within given wavelength range.
370
371 @type wavelengthMin: float or 'min'
372 @param wavelengthMin: minimum of the wavelength range
373 @type wavelengthMax: float or 'max'
374 @param wavelengthMax: maximum of the wavelength range
375 @rtype: float
376 @return: relative flux
377
378 """
379
380 if wavelengthMin == 'min':
381 wavelengthMin = self.wavelength.min()
382 if wavelengthMax == 'max':
383 wavelengthMax = self.wavelength.max()
384
385 mask = numpy.logical_and(numpy.greater(self.wavelength, wavelengthMin),
386 numpy.less(self.wavelength, wavelengthMax))
387 flux = numpy.trapz(self.flux[mask], self.wavelength[mask])
388
389 return flux
390
392 """Smooths SED.flux with a uniform (boxcar) filter of width smoothPix.
393 Cannot be undone.
394
395 @type smoothPix: int
396 @param smoothPix: size of uniform filter applied to SED, in pixels
397
398 """
399 smoothed = ndimage.uniform_filter1d(self.flux, smoothPix)
400 self.flux = smoothed
401
403 """Redshifts the SED to redshift z.
404
405 @type z: float
406 @param z: redshift
407
408 """
409
410
411
412
413
414 self.wavelength = numpy.zeros(self.z0wavelength.shape[0])
415 self.flux = numpy.zeros(self.z0flux.shape[0])
416 self.wavelength = self.wavelength + self.z0wavelength
417 self.flux = self.flux + self.z0flux
418
419 z0TotalFlux = numpy.trapz(self.z0wavelength, self.z0flux)
420 self.wavelength = self.wavelength * (1.0 + z)
421 zTotalFlux = numpy.trapz(self.wavelength, self.flux)
422 self.flux = self.flux * (z0TotalFlux / zTotalFlux)
423 self.z = z
424
425 - def normalise(self, minWavelength='min', maxWavelength='max'):
426 """Normalises the SED such that the area under the specified wavelength
427 range is equal to 1.
428
429 @type minWavelength: float or 'min'
430 @param minWavelength: minimum wavelength of range over which to
431 normalise SED
432 @type maxWavelength: float or 'max'
433 @param maxWavelength: maximum wavelength of range over which to
434 normalise SED
435
436 """
437 if minWavelength == 'min':
438 minWavelength = self.wavelength.min()
439 if maxWavelength == 'max':
440 maxWavelength = self.wavelength.max()
441
442 lowCut = numpy.greater(self.wavelength, minWavelength)
443 highCut = numpy.less(self.wavelength, maxWavelength)
444 totalCut = numpy.logical_and(lowCut, highCut)
445 sedFluxSlice = self.flux[totalCut]
446 sedWavelengthSlice = self.wavelength[totalCut]
447
448 self.flux = self.flux / numpy.trapz(
449 abs(sedFluxSlice), sedWavelengthSlice)
450
452 """Normalises the SED to match the flux equivalent to the given AB
453 magnitude in the given passband.
454
455 @type ABMag: float
456 @param ABMag: AB magnitude to which the SED is to be normalised at the
457 given passband
458 @type passband: an L{Passband} object
459 @param passband: passband at which normalisation to AB magnitude is
460 calculated
461
462 """
463
464 magFlux = mag2Flux(ABMag, 0.0, passband)
465 sedFlux = self.calcFlux(passband)
466 norm = magFlux[0] / sedFlux
467 self.flux = self.flux * norm
468 self.z0flux = self.z0flux * norm
469
470 - def matchFlux(self, matchSED, minWavelength, maxWavelength):
471 """Matches the flux in the wavelength range given by minWavelength,
472 maxWavelength to the flux in the same region in matchSED. Useful for
473 plotting purposes.
474
475 @type matchSED: L{SED} object
476 @param matchSED: SED to match flux to
477 @type minWavelength: float
478 @param minWavelength: minimum of range in which to match flux of
479 current SED to matchSED
480 @type maxWavelength: float
481 @param maxWavelength: maximum of range in which to match flux of
482 current SED to matchSED
483
484 """
485
486 interpMatch = interpolate.interp1d(matchSED.wavelength,
487 matchSED.flux,
488 kind='linear')
489 interpSelf = interpolate.interp1d(self.wavelength,
490 self.flux,
491 kind='linear')
492
493 wavelengthRange = numpy.arange(minWavelength, maxWavelength, 5.0)
494
495 matchFlux = numpy.trapz(interpMatch(wavelengthRange), wavelengthRange)
496 selfFlux = numpy.trapz(interpSelf(wavelengthRange), wavelengthRange)
497
498 self.flux = self.flux * (matchFlux / selfFlux)
499
501 """Calculates flux in the given passband.
502
503 @type passband: L{Passband} object
504 @param passband: filter passband through which to calculate the flux
505 from the SED
506 @rtype: float
507 @return: flux
508
509 """
510 lowCut = numpy.greater(self.wavelength, passband.wavelength.min())
511 highCut = numpy.less(self.wavelength, passband.wavelength.max())
512 totalCut = numpy.logical_and(lowCut, highCut)
513 sedFluxSlice = self.flux[totalCut]
514 sedWavelengthSlice = self.wavelength[totalCut]
515
516
517
518 sedInBand = passband.interpolator(sedWavelengthSlice) * sedFluxSlice
519 totalFlux = numpy.trapz(sedInBand * sedWavelengthSlice,
520 sedWavelengthSlice)
521 totalFlux = totalFlux /\
522 numpy.trapz(passband.interpolator(sedWavelengthSlice) *
523 sedWavelengthSlice, sedWavelengthSlice)
524
525 return totalFlux
526
527 - def calcMag(self, passband, addDistanceModulus=True, magType="Vega"):
528 """Calculates magnitude in the given passband. If addDistanceModulus ==
529 True, then the distance modulus (5.0*log10*(dl*1e5), where dl is the
530 luminosity distance in Mpc at the redshift of the L{SED}) is added.
531
532 @type passband: L{Passband} object
533 @param passband: filter passband through which to calculate the
534 magnitude from the SED
535 @type addDistanceModulus: bool
536 @param addDistanceModulus: if True, adds 5.0*log10*(dl*1e5) to the mag
537 returned, where dl is the luminosity distance (Mpc) corresponding to
538 the SED z
539 @type magType: string
540 @param magType: either "Vega" or "AB"
541 @rtype: float
542 @return: magnitude through the given passband on the specified
543 magnitude system
544
545 """
546 f1 = self.calcFlux(passband)
547 if magType == "Vega":
548 f2 = VEGA.calcFlux(passband)
549 elif magType == "AB":
550 f2 = AB.calcFlux(passband)
551
552 mag = -2.5 * math.log10(f1 / f2)
553 if magType == "Vega":
554
555 mag += 0.026
556
557 if self.z > 0.0 and addDistanceModulus:
558 appMag = 5.0 * math.log10(astCalc.dl(self.z) * 1e5) + mag
559 else:
560 appMag = mag
561
562 return appMag
563
564 - def calcColour(self, passband1, passband2, magType="Vega"):
565 """Calculates the colour passband1-passband2.
566
567 @type passband1: L{Passband} object
568 @param passband1: filter passband through which to calculate the first
569 magnitude
570 @type passband2: L{Passband} object
571 @param passband1: filter passband through which to calculate the second
572 magnitude
573 @type magType: string
574 @param magType: either "Vega" or "AB"
575 @rtype: float
576 @return: colour defined by passband1 - passband2 on the specified
577 magnitude system
578
579 """
580 mag1 = self.calcMag(passband1,
581 magType=magType,
582 addDistanceModulus=True)
583 mag2 = self.calcMag(passband2,
584 magType=magType,
585 addDistanceModulus=True)
586
587 colour = mag1 - mag2
588 return colour
589
591 """This is a convenience function for pulling out fluxes from a SED for
592 a given set of passbands
593 in the same format as made by L{mags2SEDDict} - designed to make
594 fitting code simpler.
595
596 @type passbands: list of L{Passband} objects
597 @param passbands: list of passbands through which fluxes will be
598 calculated
599
600 """
601
602 flux = []
603 wavelength = []
604 for p in passbands:
605 flux.append(self.calcFlux(p))
606 wavelength.append(p.effectiveWavelength())
607
608 SEDDict = {}
609 SEDDict['flux'] = numpy.array(flux)
610 SEDDict['wavelength'] = numpy.array(wavelength)
611
612 return SEDDict
613
615 """Applies the Calzetti et al. 2000 (ApJ, 533, 682) extinction law to
616 the SED with the given E(B-V) amount of extinction. R_v' = 4.05 is
617 assumed (see equation (5) of Calzetti et al.).
618
619 @type EBMinusV: float
620 @param EBMinusV: extinction E(B-V), in magnitudes
621
622 """
623
624 self.EBMinusV = EBMinusV
625
626
627 self.z0flux = self.intrinsic_z0flux
628
629
630 if EBMinusV > 0:
631
632
633
634 RvPrime = 4.05
635 shortWavelengthMask =\
636 numpy.logical_and(numpy.greater_equal(self.z0wavelength, 1200),
637 numpy.less(self.z0wavelength, 6300))
638 longWavelengthMask =\
639 numpy.logical_and(numpy.greater_equal(self.z0wavelength, 6300),
640 numpy.less_equal(self.z0wavelength, 22000))
641 wavelengthMicrons = numpy.array(self.z0wavelength / 10000.0,
642 dtype=numpy.float64)
643 kPrime = numpy.zeros(self.z0wavelength.shape[0],
644 dtype=numpy.float64)
645 kPrimeLong = (2.659 * (-1.857 + 1.040 / wavelengthMicrons)) +\
646 RvPrime
647 kPrimeShort = (2.659 * (-2.156 + 1.509 / wavelengthMicrons -
648 0.198 / wavelengthMicrons**2 + 0.011 /
649 wavelengthMicrons**3)) + RvPrime
650 kPrime[longWavelengthMask] = kPrimeLong[longWavelengthMask]
651 kPrime[shortWavelengthMask] = kPrimeShort[shortWavelengthMask]
652
653
654
655 try:
656 interpolator = interpolate.interp1d(self.z0wavelength,
657 kPrimeShort,
658 kind='linear')
659 slope = (interpolator(1100.0) - interpolator(1200.0)) / (
660 1100.0 - 1200.0)
661 intercept = interpolator(1200.0) - (slope * 1200.0)
662 mask = numpy.less(self.z0wavelength, 1200.0)
663 kPrime[mask] = slope * self.z0wavelength[mask] + intercept
664
665
666 interpolator = interpolate.interp1d(self.z0wavelength,
667 kPrimeLong,
668 kind='linear')
669 slope = (interpolator(21900.0) - interpolator(22000.0)) / (
670 21900.0 - 22000.0)
671 intercept = interpolator(21900.0) - (slope * 21900.0)
672 mask = numpy.greater(self.z0wavelength, 22000.0)
673 kPrime[mask] = slope * self.z0wavelength[mask] + intercept
674 except:
675 raise Exception("This SED has a wavelength range that doesn't "
676 "cover ~1200-22000 Angstroms")
677
678
679 kPrime[numpy.less_equal(kPrime, 0.0)] = 1e-6
680
681 reddening = numpy.power(10, 0.4 * EBMinusV * kPrime)
682 self.z0flux = self.z0flux / reddening
683
684 self.redshift(self.z)
685
686
687
689 """This class stores the SED of Vega, used for calculation of magnitudes on the Vega system.
690
691 The Vega SED used is taken from Bohlin 2007
692 (http://adsabs.harvard.edu/abs/2007ASPC..364..315B), and is available from
693 the STScI CALSPEC library
694 (http://www.stsci.edu/hst/observatory/cdbs/calspec.html).
695
696 """
697
699
700 VEGA_SED_PATH = astLib.__path__[
701 0] + os.path.sep + "data" + os.path.sep + "bohlin2006_Vega.sed"
702
703 inFile = open(VEGA_SED_PATH, "r")
704 lines = inFile.readlines()
705
706 wavelength = []
707 flux = []
708 for line in lines:
709
710 if line[0] != "#" and len(line) > 3:
711
712 bits = line.split()
713 flux.append(float(bits[1]))
714 wavelength.append(float(bits[0]))
715
716 self.wavelength = numpy.array(wavelength)
717 self.flux = numpy.array(flux, dtype=numpy.float64)
718
719
720
721 self.z0wavelength = numpy.array(wavelength)
722 self.z0flux = numpy.array(flux, dtype=numpy.float64)
723 self.z = 0.0
724
725
726
727
728
729
730
732 """This class describes a stellar population model, either a Simple Stellar
733 Population (SSP) or a Composite Stellar Population (CSP), such as the
734 models of Bruzual & Charlot 2003 or Maraston 2005.
735
736 The constructor for this class can be used for generic SSPs or CSPs stored
737 in white space delimited text files, containing columns for age,
738 wavelength, and flux. Columns are counted from 0 ... n. Lines starting
739 with # are ignored.
740
741 The classes L{M05Model} (for Maraston 2005 models), L{BC03Model} (for
742 Bruzual & Charlot 2003 models), and L{P09Model} (for Percival et al. 2009
743 models) are derived from this class. The only difference between them is
744 the code used to load in the model data.
745
746 """
747
748 - def __init__(self,
749 fileName,
750 ageColumn=0,
751 wavelengthColumn=1,
752 fluxColumn=2):
753
754 inFile = open(fileName, "r")
755 lines = inFile.readlines()
756 inFile.close()
757
758 self.fileName = fileName
759
760
761 self.ages = []
762 self.wavelengths = []
763 for line in lines:
764 if line[0] != "#" and len(line) > 3:
765 bits = line.split()
766 age = float(bits[ageColumn])
767 wavelength = float(bits[wavelengthColumn])
768 if age not in self.ages:
769 self.ages.append(age)
770 if wavelength not in self.wavelengths:
771 self.wavelengths.append(wavelength)
772
773
774 self.fluxGrid = numpy.zeros([len(self.ages), len(self.wavelengths)])
775 for line in lines:
776 if line[0] != "#" and len(line) > 3:
777 bits = line.split()
778 sedAge = float(bits[ageColumn])
779 sedWavelength = float(bits[wavelengthColumn])
780 sedFlux = float(bits[fluxColumn])
781
782 row = self.ages.index(sedAge)
783 column = self.wavelengths.index(sedWavelength)
784
785 self.fluxGrid[row][column] = sedFlux
786
787 - def getSED(self, ageGyr, z=0.0, normalise=False, label=None):
788 """Extract a SED for given age. Do linear interpolation between models
789 if necessary.
790
791 @type ageGyr: float
792 @param ageGyr: age of the SED in Gyr
793 @type z: float
794 @param z: redshift the SED from z = 0 to z = z
795 @type normalise: bool
796 @param normalise: normalise the SED to have area 1
797 @rtype: L{SED} object
798 @return: SED
799
800 """
801
802 if ageGyr in self.ages:
803
804 flux = self.fluxGrid[self.ages.index(ageGyr)]
805 sed = SED(self.wavelengths,
806 flux,
807 z=z,
808 normalise=normalise,
809 label=label)
810 return sed
811
812 else:
813
814
815 flux = []
816 for i in range(len(self.wavelengths)):
817 interpolator = interpolate.interp1d(self.ages,
818 self.fluxGrid[:, i],
819 kind='linear')
820 sedFlux = interpolator(ageGyr)
821 flux.append(sedFlux)
822 sed = SED(self.wavelengths,
823 flux,
824 z=z,
825 normalise=normalise,
826 label=label)
827 return sed
828
829 - def getColourEvolution(self,
830 passband1,
831 passband2,
832 zFormation,
833 zStepSize=0.05,
834 magType="Vega"):
835 """Calculates the evolution of the colour observed through passband1 -
836 passband2 for the StellarPopulation with redshift, from z = 0 to z =
837 zFormation.
838
839 @type passband1: L{Passband} object
840 @param passband1: filter passband through which to calculate the first
841 magnitude
842 @type passband2: L{Passband} object
843 @param passband2: filter passband through which to calculate the second
844 magnitude
845 @type zFormation: float
846 @param zFormation: formation redshift of the StellarPopulation
847 @type zStepSize: float
848 @param zStepSize: size of interval in z at which to calculate model
849 colours
850 @type magType: string
851 @param magType: either "Vega" or "AB"
852 @rtype: dictionary
853 @return: dictionary of numpy.arrays in format {'z', 'colour'}
854
855 """
856
857 zSteps = int(math.ceil(zFormation / zStepSize))
858 zData = []
859 colourData = []
860 for i in range(1, zSteps):
861 zc = i * zStepSize
862 age = astCalc.tl(zFormation) - astCalc.tl(zc)
863 sed = self.getSED(age, z=zc)
864 colour = sed.calcColour(passband1, passband2, magType=magType)
865 zData.append(zc)
866 colourData.append(colour)
867
868 zData = numpy.array(zData)
869 colourData = numpy.array(colourData)
870
871 return {'z': zData, 'colour': colourData}
872
873 - def getMagEvolution(self,
874 passband,
875 magNormalisation,
876 zNormalisation,
877 zFormation,
878 zStepSize=0.05,
879 onePlusZSteps=False,
880 magType="Vega"):
881 """Calculates the evolution with redshift (from z = 0 to z =
882 zFormation) of apparent magnitude in the observed frame through the
883 passband for the StellarPopulation, normalised to magNormalisation
884 (apparent) at z = zNormalisation.
885
886 @type passband: L{Passband} object
887 @param passband: filter passband through which to calculate the
888 magnitude
889 @type magNormalisation: float
890 @param magNormalisation: sets the apparent magnitude of the SED at
891 zNormalisation
892 @type zNormalisation: float
893 @param zNormalisation: the redshift at which the magnitude
894 normalisation is carried out
895 @type zFormation: float
896 @param zFormation: formation redshift of the StellarPopulation
897 @type zStepSize: float
898 @param zStepSize: size of interval in z at which to calculate model
899 magnitudes
900 @type onePlusZSteps: bool
901 @param onePlusZSteps: if True, zSteps are (1+z)*zStepSize, otherwise
902 zSteps are linear
903 @type magType: string
904 @param magType: either "Vega" or "AB"
905 @rtype: dictionary
906 @return: dictionary of numpy.arrays in format {'z', 'mag'}
907
908 """
909
910
911 zSteps = int(math.ceil(zFormation / zStepSize))
912 zData = []
913 magData = []
914 absMagData = []
915 zc0 = 0.0
916 for i in range(1, zSteps):
917 if not onePlusZSteps:
918 zc = i * zStepSize
919 else:
920 zc = zc0 + (1 + zc0) * zStepSize
921 zc0 = zc
922 if zc >= zFormation:
923 break
924 age = astCalc.tl(zFormation) - astCalc.tl(zc)
925 sed = self.getSED(age, z=zc)
926 mag = sed.calcMag(passband,
927 magType=magType,
928 addDistanceModulus=True)
929 zData.append(zc)
930 magData.append(mag)
931 absMagData.append(sed.calcMag(passband, addDistanceModulus=False))
932
933 zData = numpy.array(zData)
934 magData = numpy.array(magData)
935
936
937 interpolator = interpolate.interp1d(zData, magData, kind='linear')
938 modelNormMag = interpolator(zNormalisation)
939 normConstant = magNormalisation - modelNormMag
940 magData = magData + normConstant
941
942 return {'z': zData, 'mag': magData}
943
950 """Calculates the evolution correction in magnitudes in the rest frame
951 through the passband from redshift zFrom to redshift zTo, where the
952 stellarPopulation is assumed to be formed at redshift zFormation.
953
954 @type zFrom: float
955 @param zFormation: redshift to evolution correct from
956 @type zTo: float
957 @param zTo: redshift to evolution correct to
958 @type zFormation: float
959 @param zFormation: formation redshift of the StellarPopulation
960 @type passband: L{Passband} object
961 @param passband: filter passband through which to calculate magnitude
962 @type magType: string
963 @param magType: either "Vega" or "AB"
964 @rtype: float
965 @return: evolution correction in magnitudes in the rest frame
966
967 """
968
969 ageFrom = astCalc.tl(zFormation) - astCalc.tl(zFrom)
970 ageTo = astCalc.tl(zFormation) - astCalc.tl(zTo)
971
972 fromSED = self.getSED(ageFrom)
973 toSED = self.getSED(ageTo)
974
975 fromMag = fromSED.calcMag(passband,
976 magType=magType,
977 addDistanceModulus=False)
978 toMag = toSED.calcMag(passband,
979 magType=magType,
980 addDistanceModulus=False)
981
982 return fromMag - toMag
983
984
985
987 """This class describes a Maraston 2005 stellar population model. To load a
988 composite stellar population model (CSP) for a tau = 0.1 Gyr burst of star
989 formation, solar metallicity, Salpeter IMF:
990
991 m05csp = astSED.M05Model(M05_DIR+"/csp_e_0.10_z02_salp.sed_agb")
992
993 where M05_DIR is set to point to the directory where the Maraston 2005
994 models are stored on your system.
995
996 The file format of the Maraston 2005 simple stellar poulation (SSP) models
997 is different to the file format used for the CSPs, and this needs to be
998 specified using the fileType parameter. To load a SSP with solar
999 metallicity, red horizontal branch morphology:
1000
1001 m05ssp = astSED.M05Model(M05_DIR+"/sed.ssz002.rhb", fileType = "ssp")
1002
1003 The wavelength units of SEDs from M05 models are Angstroms, with flux in
1004 units of erg/s/Angstrom.
1005
1006 """
1007
1008 - def __init__(self, fileName, fileType="csp"):
1009
1010 self.modelFamily = "M05"
1011
1012 inFile = open(fileName, "r")
1013 lines = inFile.readlines()
1014 inFile.close()
1015
1016 self.fileName = fileName
1017
1018 if fileType == "csp":
1019 ageColumn = 0
1020 wavelengthColumn = 1
1021 fluxColumn = 2
1022 elif fileType == "ssp":
1023 ageColumn = 0
1024 wavelengthColumn = 2
1025 fluxColumn = 3
1026 else:
1027 raise Exception("fileType must be 'ssp' or 'csp'")
1028
1029
1030 self.ages = []
1031 self.wavelengths = []
1032 for line in lines:
1033 if line[0] != "#" and len(line) > 3:
1034 bits = line.split()
1035 age = float(bits[ageColumn])
1036 wavelength = float(bits[wavelengthColumn])
1037 if age not in self.ages:
1038 self.ages.append(age)
1039 if wavelength not in self.wavelengths:
1040 self.wavelengths.append(wavelength)
1041
1042
1043 self.fluxGrid = numpy.zeros([len(self.ages), len(self.wavelengths)])
1044 for line in lines:
1045 if line[0] != "#" and len(line) > 3:
1046 bits = line.split()
1047 sedAge = float(bits[ageColumn])
1048 sedWavelength = float(bits[wavelengthColumn])
1049 sedFlux = float(bits[fluxColumn])
1050
1051 row = self.ages.index(sedAge)
1052 column = self.wavelengths.index(sedWavelength)
1053
1054 self.fluxGrid[row][column] = sedFlux
1055
1056
1057
1059 """This class describes a Bruzual & Charlot 2003 stellar population model,
1060 extracted from a GALAXEV .ised file using the galaxevpl program that is
1061 included in GALAXEV. The file format is white space delimited, with
1062 wavelength in the first column. Subsequent columns contain the model fluxes
1063 for SEDs of different ages, as specified when running galaxevpl. The age
1064 corresponding to each flux column is taken from the comment line beginning
1065 "# Age (yr)", and is converted to Gyr.
1066
1067 For example, to load a tau = 0.1 Gyr burst of star formation, solar
1068 metallicity, Salpeter IMF model stored in a file (created by galaxevpl)
1069 called "csp_lr_solar_0p1Gyr.136":
1070
1071 bc03model = BC03Model("csp_lr_solar_0p1Gyr.136")
1072
1073 The wavelength units of SEDs from BC03 models are Angstroms. Flux is
1074 converted into units of erg/s/Angstrom (the units in the files output by
1075 galaxevpl are LSun/Angstrom).
1076
1077 """
1078
1080
1081 self.modelFamily = "BC03"
1082 self.fileName = fileName
1083
1084 inFile = open(fileName, "r")
1085 lines = inFile.readlines()
1086 inFile.close()
1087
1088
1089 self.ages = []
1090 for line in lines:
1091 if line.find("# Age (yr)") != -1:
1092 rawAges = line[line.find("# Age (yr)") + 10:].split()
1093 for age in rawAges:
1094 self.ages.append(float(age) / 1e9)
1095
1096
1097
1098 lambdaLinesCount = 0
1099 startFluxDataLine = None
1100 for i in range(len(lines)):
1101 line = lines[i]
1102 if "# Lambda(A)" in line:
1103 lambdaLinesCount = lambdaLinesCount + 1
1104 if line[0] != "#" and len(line) > 3 and startFluxDataLine is None:
1105 startFluxDataLine = i
1106 self.wavelengths = []
1107 for i in range(startFluxDataLine, len(lines), lambdaLinesCount):
1108 line = lines[i]
1109 bits = line.split()
1110 self.wavelengths.append(float(bits[0]))
1111
1112
1113 self.fluxGrid = numpy.zeros([len(self.ages), len(self.wavelengths)])
1114 for i in range(startFluxDataLine, len(lines), lambdaLinesCount):
1115 line = lines[i]
1116 bits = []
1117 for k in range(i, i + lambdaLinesCount):
1118 bits = bits + lines[k].split()
1119 ageFluxes = bits[1:]
1120 sedWavelength = float(bits[0])
1121 column = self.wavelengths.index(sedWavelength)
1122 for row in range(len(ageFluxes)):
1123 sedFlux = float(ageFluxes[row])
1124 self.fluxGrid[row][column] = sedFlux
1125
1126
1127 self.fluxGrid = self.fluxGrid * 3.826e33
1128
1129
1130
1132 """This class describes a Percival et al 2009 (BaSTI;
1133 http://albione.oa-teramo.inaf.it) stellar population model. We assume that
1134 the synthetic spectra for each model are unpacked under the directory
1135 pointed to by fileName.
1136
1137 The wavelength units of SEDs from P09 models are converted to Angstroms.
1138 Flux is converted into units of erg/s/Angstrom (the units in the BaSTI
1139 low-res spectra are 4.3607e-33 erg/s/m).
1140
1141 """
1142
1144
1145 self.modelFamily = "P09"
1146
1147 files = glob.glob(fileName + os.path.sep + "*.t??????")
1148
1149 self.fileName = fileName
1150
1151
1152 extensionAgeMap = {}
1153 self.ages = []
1154 for f in files:
1155 ext = f.split(".")[-1]
1156 ageGyr = float(f[-5:]) / 1e3
1157 self.ages.append(ageGyr)
1158 extensionAgeMap[ext] = ageGyr
1159 self.ages.sort()
1160
1161
1162 self.wavelengths = None
1163 self.fluxGrid = None
1164 for i in range(len(self.ages)):
1165 for e in extensionAgeMap.keys():
1166 if extensionAgeMap[e] == self.ages[i]:
1167 inFileName = glob.glob(fileName + os.path.sep + "*." + e)[
1168 0]
1169 inFile = open(inFileName, "r")
1170 lines = inFile.readlines()
1171 inFile.close()
1172 wavelength = []
1173 flux = []
1174 for line in lines:
1175 bits = line.split()
1176 wavelength.append(
1177 float(bits[0]) *
1178 10.0)
1179 flux.append(float(bits[1]))
1180 if self.wavelengths is None:
1181 self.wavelengths = wavelength
1182 if self.fluxGrid is None:
1183 self.fluxGrid = numpy.zeros(
1184 [len(self.ages), len(self.wavelengths)])
1185 self.fluxGrid[i] = flux
1186
1187
1188
1189 self.fluxGrid = self.fluxGrid / 4.3607e-33 / 1e10
1190
1191
1192
1193 -def makeModelSEDDictList(modelList,
1194 z,
1195 passbandsList,
1196 labelsList=[],
1197 EBMinusVList=[0.0],
1198 forceYoungerThanUniverse=True):
1199 """This routine makes a list of SEDDict dictionaries (see L{mags2SEDDict})
1200 for fitting using L{fitSEDDict}. This speeds up the fitting as this allows
1201 us to calculate model SED magnitudes only once, if all objects to be fitted
1202 are at the same redshift. We add some meta data to the modelSEDDicts (e.g.
1203 the model file names).
1204
1205 The effect of extinction by dust (assuming the Calzetti et al. 2000 law)
1206 can be included by giving a list of E(B-V) values.
1207
1208 If forceYoungerThanUniverse == True, ages which are older than the universe
1209 at the given z will not be included.
1210
1211 @type modelList: list of L{StellarPopulation} model objects
1212 @param modelList: list of StellarPopulation models to include
1213 @type z: float
1214 @param z: redshift to apply to all stellar population models in modelList
1215 @type EBMinusVList: list
1216 @param EBMinusVList: list of E(B-V) extinction values to apply to all
1217 models, in magnitudes
1218 @type labelsList: list
1219 @param labelsList: optional list used for labelling passbands in output
1220 SEDDicts
1221 @type forceYoungerThanUniverse: bool
1222 @param forceYoungerThanUniverse: if True, do not allow models that exceed
1223 the age of the universe at z
1224 @rtype: list
1225 @return: list of dictionaries containing model fluxes, to be used as input
1226 to L{fitSEDDict}.
1227
1228 """
1229
1230
1231 if EBMinusVList == []:
1232 EBMinusVList = [0.0]
1233
1234 modelSEDDictList = []
1235 for m in range(len(modelList)):
1236 testAges = numpy.array(modelList[m].ages)
1237 if forceYoungerThanUniverse:
1238 testAges = testAges[numpy.logical_and(
1239 numpy.less(testAges, astCalc.tz(z)), numpy.greater(testAges,
1240 0))]
1241 for t in testAges:
1242 s = modelList[m].getSED(
1243 t,
1244 z=z,
1245 label=modelList[m].fileName + " - age=" + str(t) + " Gyr")
1246 for EBMinusV in EBMinusVList:
1247 try:
1248 s.extinctionCalzetti(EBMinusV)
1249 except:
1250 raise Exception(
1251 "Model %s has a wavelength range that doesn't cover ~1200-22000 Angstroms"
1252 % (modelList[m].fileName))
1253 modelSEDDict = s.getSEDDict(passbandsList)
1254 modelSEDDict['labels'] = labelsList
1255 modelSEDDict['E(B-V)'] = EBMinusV
1256 modelSEDDict['ageGyr'] = t
1257 modelSEDDict['z'] = z
1258 modelSEDDict['fileName'] = modelList[m].fileName
1259 modelSEDDict['modelListIndex'] = m
1260 modelSEDDictList.append(modelSEDDict)
1261
1262 return modelSEDDictList
1263
1264
1265
1267 """Fits the given SED dictionary (made using L{mags2SEDDict}) with the
1268 given list of model SED dictionaries. The latter should be made using
1269 L{makeModelSEDDictList}, and entries for fluxes should correspond directly
1270 between the model and SEDDict.
1271
1272 Returns a dictionary with best fit values.
1273
1274 @type SEDDict: dictionary, in format of L{mags2SEDDict}
1275 @param SEDDict: dictionary of observed fluxes and uncertainties, in format
1276 of L{mags2SEDDict}
1277 @type modelSEDDictList: list of dictionaries, in format of
1278 L{makeModelSEDDictList}
1279 @param modelSEDDictList: list of dictionaries containing fluxes of models
1280 to be fitted to the observed fluxes listed in the SEDDict. This should be
1281 made using L{makeModelSEDDictList}.
1282 @rtype: dictionary
1283 @return: results of the fitting - keys:
1284 - 'minChiSq': minimum chi squared value of best fit
1285 - 'chiSqContrib': corresponding contribution at each passband to
1286 the minimum chi squared value
1287 - 'ageGyr': the age in Gyr of the best fitting model
1288 - 'modelFileName': the file name of the stellar population model
1289 corresponding to the best fit
1290 - 'modelListIndex': the index of the best fitting model in the
1291 input modelSEDDictList
1292 - 'norm': the normalisation that the best fit model should be
1293 multiplied by to match the SEDDict
1294 - 'z': the redshift of the best fit model
1295 - 'E(B-V)': the extinction, E(B-V), in magnitudes, of the best fit
1296 model
1297
1298 """
1299
1300 modelFlux = []
1301 for modelSEDDict in modelSEDDictList:
1302 modelFlux.append(modelSEDDict['flux'])
1303 modelFlux = numpy.array(modelFlux)
1304 sedFlux = numpy.array([SEDDict['flux']] * len(modelSEDDictList))
1305 sedFluxErr = numpy.array([SEDDict['fluxErr']] * len(modelSEDDictList))
1306
1307
1308 norm = numpy.sum(
1309 (modelFlux * sedFlux) /
1310 (sedFluxErr**2), axis=1) / numpy.sum(modelFlux**2 / sedFluxErr**2,
1311 axis=1)
1312 norms = numpy.array([norm] * modelFlux.shape[1]).transpose()
1313 chiSq = numpy.sum(((sedFlux - norms * modelFlux)**2) / sedFluxErr**2,
1314 axis=1)
1315 chiSq[numpy.isnan(
1316 chiSq)] = 1e6
1317 minChiSq = chiSq.min()
1318 bestMatchIndex = numpy.equal(chiSq, minChiSq).nonzero()[0][0]
1319 bestNorm = norm[bestMatchIndex]
1320 bestChiSq = minChiSq
1321 bestChiSqContrib = ((sedFlux[bestMatchIndex] - norms[bestMatchIndex] *
1322 modelFlux[bestMatchIndex])**2) /\
1323 sedFluxErr[bestMatchIndex]**2
1324
1325 resultsDict = {'minChiSq': bestChiSq,
1326 'chiSqContrib': bestChiSqContrib,
1327 'allChiSqValues': chiSq,
1328 'ageGyr': modelSEDDictList[bestMatchIndex]['ageGyr'],
1329 'modelFileName':
1330 modelSEDDictList[bestMatchIndex]['fileName'],
1331 'modelListIndex':
1332 modelSEDDictList[bestMatchIndex]['modelListIndex'],
1333 'norm': bestNorm,
1334 'z': modelSEDDictList[bestMatchIndex]['z'],
1335 'E(B-V)': modelSEDDictList[bestMatchIndex]['E(B-V)']}
1336
1337 return resultsDict
1338
1339
1340
1342 """Takes a set of corresponding AB magnitudes, uncertainties, and
1343 passbands, and returns a dictionary with keys 'flux', 'fluxErr'
1344 'wavelength'. Fluxes are in units of erg/s/cm^2/Angstrom, wavelength in
1345 Angstroms. These dictionaries are the staple diet of the L{fitSEDDict}
1346 routine.
1347
1348 @type ABMags: list or numpy array
1349 @param ABMags: AB magnitudes, specified in corresponding order to passbands
1350 and ABMagErrs
1351 @type ABMagErrs: list or numpy array
1352 @param ABMagErrs: AB magnitude errors, specified in corresponding order to
1353 passbands and ABMags
1354 @type passbands: list of L{Passband} objects
1355 @param passbands: passband objects, specified in corresponding order to
1356 ABMags and ABMagErrs
1357 @rtype: dictionary
1358 @return: dictionary with keys {'flux', 'fluxErr', 'wavelength'}, suitable
1359 for input to L{fitSEDDict}
1360
1361 """
1362
1363 flux = []
1364 fluxErr = []
1365 wavelength = []
1366 for m, e, p in zip(ABMags, ABMagErrs, passbands):
1367 f, err = mag2Flux(m, e, p)
1368 flux.append(f)
1369 fluxErr.append(err)
1370 wavelength.append(p.effectiveWavelength())
1371
1372 SEDDict = {}
1373 SEDDict['flux'] = numpy.array(flux)
1374 SEDDict['fluxErr'] = numpy.array(fluxErr)
1375 SEDDict['wavelength'] = numpy.array(wavelength)
1376
1377 return SEDDict
1378
1379
1380
1381 -def mag2Flux(ABMag, ABMagErr, passband):
1382 """Converts given AB magnitude and uncertainty into flux, in
1383 erg/s/cm^2/Angstrom.
1384
1385 @type ABMag: float
1386 @param ABMag: magnitude on AB system in passband
1387 @type ABMagErr: float
1388 @param ABMagErr: uncertainty in AB magnitude in passband
1389 @type passband: L{Passband} object
1390 @param passband: L{Passband} object at which ABMag was measured
1391 @rtype: list
1392 @return: [flux, fluxError], in units of erg/s/cm^2/Angstrom
1393
1394 """
1395
1396 fluxJy = (10**23.0) * 10**(-(ABMag + 48.6) / 2.5)
1397 aLambda = 3e-13
1398 effLMicron = passband.effectiveWavelength() * (1e-10 / 1e-6)
1399 fluxWLUnits = aLambda * fluxJy / effLMicron**2
1400
1401 fluxJyErr = (10**23.0) * 10**(-(ABMag - ABMagErr + 48.6) / 2.5)
1402 fluxWLUnitsErr = aLambda * fluxJyErr / effLMicron**2
1403 fluxWLUnitsErr = fluxWLUnitsErr - fluxWLUnits
1404
1405 return [fluxWLUnits, fluxWLUnitsErr]
1406
1407
1408
1410 """Converts given flux and uncertainty in erg/s/cm^2/Angstrom into AB
1411 magnitudes.
1412
1413 @type flux: float
1414 @param flux: flux in erg/s/cm^2/Angstrom in passband
1415 @type fluxErr: float
1416 @param fluxErr: uncertainty in flux in passband, in erg/s/cm^2/Angstrom
1417 @type passband: L{Passband} object
1418 @param passband: L{Passband} object at which ABMag was measured
1419 @rtype: list
1420 @return: [ABMag, ABMagError], in AB magnitudes
1421
1422 """
1423
1424
1425 aLambda = 3e-13
1426 effLMicron = passband.effectiveWavelength() * (1e-10 / 1e-6)
1427
1428 fluxJy = (flux * effLMicron**2) / aLambda
1429 mag = -2.5 * numpy.log10(fluxJy / 10**23) - 48.6
1430
1431 fluxErrJy = (fluxErr * effLMicron**2) / aLambda
1432 magErr = mag - (-2.5 * numpy.log10((fluxJy + fluxErrJy) / 10**23) - 48.6)
1433
1434 return [mag, magErr]
1435
1436
1437
1439 """Converts an AB magnitude into flux density in Jy
1440
1441 @type ABMag: float
1442 @param ABMag: AB magnitude
1443 @rtype: float
1444 @return: flux density in Jy
1445
1446 """
1447
1448 fluxJy = ((10**23) * 10**(-(float(ABMag) + 48.6) / 2.5))
1449
1450 return fluxJy
1451
1452
1453
1455 """Converts flux density in Jy into AB magnitude
1456
1457 @type fluxJy: float
1458 @param fluxJy: flux density in Jy
1459 @rtype: float
1460 @return: AB magnitude
1461
1462 """
1463
1464 ABMag = -2.5 * (numpy.log10(fluxJy) - 23.0) - 48.6
1465
1466 return ABMag
1467
1468
1469
1470
1471 VEGA = VegaSED()
1472
1473
1474 AB = SED(wavelength=numpy.logspace(1, 8, 1e5), flux=numpy.ones(1e6))
1475 AB.flux = (3e-5 * 3631) / (AB.wavelength**2)
1476 AB.z0flux = AB.flux[:]
1477
1478
1479 SOL = SED()
1480 SOL.loadFromFile(astLib.__path__[0] + os.path.sep + "data" + os.path.sep +
1481 "sun_reference_stis_001.ascii")
1482