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