Package astLib :: Module astSED
[hide private]
[frames] | no frames]

Source Code for Module astLib.astSED

   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  #------------------------------------------------------------------------------------------------------------ 
44 -class Passband:
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', wavelengthColumn = 0, transmissionColumn = 1):
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[transmissionColumn])) 71 wavelength.append(float(bits[wavelengthColumn])) 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 # Sort into ascending order of wavelength otherwise normalisation will be wrong 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 # Store a ready-to-go interpolation object to speed calculation of fluxes up 100 self.interpolator=interpolate.interp1d(self.wavelength, self.transmission, kind='linear')
101
102 - def asList(self):
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
155 - def effectiveWavelength(self):
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 #------------------------------------------------------------------------------------------------------------
171 -class TopHatPassband(Passband):
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 # Store a ready-to-go interpolation object to speed calculation of fluxes up 197 self.interpolator=interpolate.interp1d(self.wavelength, self.transmission, kind='linear')
198 199 200 #------------------------------------------------------------------------------------------------------------
201 -class SED:
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 # We keep a copy of the wavelength, flux at z = 0, as it's more robust to copy that 216 # to self.wavelength, flux and redshift it, rather than repeatedly redshifting the same 217 # arrays back and forth 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 # plain text label, handy for using in photo-z codes 224 225 # Store the intrinsic (i.e. unextincted) flux in case we change extinction 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
238 - def copy(self):
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
252 - def loadFromFile(self, fileName):
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 # Sort SED so wavelength is in ascending order 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
283 - def writeToFile(self, fileName):
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
296 - def asList(self):
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 # Sensible y scale 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
359 - def smooth(self, smoothPix):
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
369 - def redshift(self, z):
370 """Redshifts the SED to redshift z. 371 372 @type z: float 373 @param z: redshift 374 375 """ 376 377 # We have to conserve energy so the area under the redshifted SED has to be equal to 378 # the area under the unredshifted SED, otherwise magnitude calculations will be wrong 379 # when comparing SEDs at different zs 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)#self.wavelength)
412
413 - def normaliseToMag(self, ABMag, passband):
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
453 - def calcFlux(self, passband):
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 # Use linear interpolation to rebin the passband to the same dimensions as the 469 # part of the SED we're interested in 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 # Add 0.026 because Vega has V=0.026 (e.g. Bohlin & Gilliland 2004) 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
529 - def getSEDDict(self, passbands):
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
550 - def extinctionCalzetti(self, EBMinusV):
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 # All done in rest frame 562 self.z0flux=self.intrinsic_z0flux 563 564 # Allow us to set EBMinusV == 0 to turn extinction off 565 if EBMinusV > 0: 566 # Note that EBMinusV is assumed to be Es as in equations (2) - (5) 567 # Note here wavelength units have to be microns for constants to make sense 568 RvPrime=4.05 # equation (5) of Calzetti et al. 2000 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 # Here we extrapolate kPrime in similar way to what HYPERZ does 587 # Short wavelengths 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 # Long wavelengths 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 # Never let go negative 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 #------------------------------------------------------------------------------------------------------------
613 -class VegaSED(SED):
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
621 - def __init__(self, normalise = False):
622 623 VEGA_SED_PATH=astLib.__path__[0]+os.path.sep+"data"+os.path.sep+"bohlin2006_Vega.sed" # from HST CALSPEC 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 # We may want to redshift reference SEDs to calculate rest-frame colors from SEDs at different zs 642 self.z0wavelength=numpy.array(wavelength) 643 self.z0flux=numpy.array(flux, dtype=numpy.float64) 644 self.z=0.0
645 646 #if normalise == True: 647 #self.flux=self.flux/numpy.trapz(self.flux, self.wavelength) 648 #self.z0flux=self.z0flux/numpy.trapz(self.z0flux, self.z0wavelength) 649 650 #------------------------------------------------------------------------------------------------------------
651 -class StellarPopulation:
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 # Extract a list of model ages and valid wavelengths from the file 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 # Construct a grid of flux - rows correspond to each wavelength, columns to age 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 # Use interpolation, iterating over each wavelength column 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 # Count upwards in z steps as interpolation doesn't work if array ordered z decreasing 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 # Do the normalisation 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
822 - def calcEvolutionCorrection(self, zFrom, zTo, zFormation, passband, magType = "Vega"):
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 #------------------------------------------------------------------------------------------------------------
854 -class M05Model(StellarPopulation):
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 # Extract a list of model ages and valid wavelengths from the file 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 # Construct a grid of flux - rows correspond to each wavelength, columns to age 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 #------------------------------------------------------------------------------------------------------------
920 -class BC03Model(StellarPopulation):
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
937 - def __init__(self, fileName):
938 939 self.modelFamily="BC03" 940 self.fileName=fileName 941 942 inFile=open(fileName, "r") 943 lines=inFile.readlines() 944 inFile.close() 945 946 # Extract a list of model ages - BC03 ages are in years, so convert to Gyr 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 # Extract a list of valid wavelengths from the file 955 # If we have many ages in the file, this is more complicated... 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 # Construct a grid of flux - rows correspond to each wavelength, columns to age 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 # Convert flux into erg/s/Angstrom - native units of galaxevpl files are LSun/Angstrom 985 self.fluxGrid=self.fluxGrid*3.826e33
986 987 #------------------------------------------------------------------------------------------------------------
988 -class P09Model(StellarPopulation):
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
998 - def __init__(self, fileName):
999 1000 self.modelFamily="P09" 1001 1002 files=glob.glob(fileName+os.path.sep+"*.t??????") 1003 1004 self.fileName=fileName 1005 1006 # Map end of filenames to ages in Gyr 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 # Construct a grid of flux - rows correspond to each wavelength, columns to age 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) # units in file are nm, not angstroms 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 # Convert flux into erg/s/Angstrom - native units in BaSTI files are 4.3607e-33 erg/s/m 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 # Otherwise if this is the case we won't actually make any model SEDDicts ... 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 #------------------------------------------------------------------------------------------------------------
1097 -def fitSEDDict(SEDDict, modelSEDDictList):
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 # Analytic expression below is for normalisation at minimum chi squared (see note book) 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 # throw these out, should check this out and handle more gracefully 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 #------------------------------------------------------------------------------------------------------------
1154 -def mags2SEDDict(ABMags, ABMagErrs, passbands):
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) # AB mag 1203 aLambda=3e-13 # for conversion to erg s-1 cm-2 angstrom-1 with lambda in microns 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) # AB mag 1208 fluxWLUnitsErr=aLambda*fluxJyErr/effLMicron**2 1209 fluxWLUnitsErr=fluxWLUnitsErr-fluxWLUnits 1210 1211 return [fluxWLUnits, fluxWLUnitsErr]
1212 1213 #------------------------------------------------------------------------------------------------------------
1214 -def flux2Mag(flux, fluxErr, passband):
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 # aLambda = 3x10-5 for effective wavelength in angstroms 1229 aLambda=3e-13 # for conversion to erg s-1 cm-2 angstrom-1 with lambda in microns 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 #------------------------------------------------------------------------------------------------------------
1241 -def mag2Jy(ABMag):
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 #------------------------------------------------------------------------------------------------------------
1257 -def Jy2Mag(fluxJy):
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 # Data 1273 VEGA=VegaSED() 1274 1275 # AB SED has constant flux density 3631 Jy 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 # Solar SED from HST CALSPEC (http://www.stsci.edu/hst/observatory/cdbs/calspec.html) 1281 SOL=SED() 1282 SOL.loadFromFile(astLib.__path__[0]+os.path.sep+"data"+os.path.sep+"sun_reference_stis_001.ascii") 1283