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