cosmolopy.luminosityfunction

1 """Galaxy luminosity functions (Schechter functions). 2 3 The `LFHistory` class implements a luminosity function history, 4 encapsulating the changes in the galaxy luminosity distribution 5 parameters as a function of redshift. 6 7 """ 8 9 import os 10 import optparse 11 12 import numpy 13 import scipy.special 14 15 import cosmolopy.reionization as cr 16 import cosmolopy.distance as cd 17 import cosmolopy.parameters as cp 18 import cosmolopy.constants as cc 19 import cosmolopy.utils as utils 20 from cosmolopy.saveable import Saveable 21 import cosmolopy.magnitudes as magnitudes 22

23 -def mass_from_sfr(sfr):

24 """Use Labbe et al. (2009) relation between stellar mass and SFR. 25 26 See arXiv:0911.1356v4 27 28 Parameters 29 ---------- 30 31 sfr: 32 the star formation rate in M_Sun / year. 33 34 Returns 35 ------- 36 37 mass: 38 stellar mass in units of M_Sun 39 """ 40 return 10**(8.7 + 1.06 * numpy.log10(sfr))

41

42 -def sfr_from_mass(mass):

43 """Use Labbe et al. (2009) relation between stellar mass and SFR. 44 45 See arXiv:0911.1356v4 46 47 Parameters 48 ---------- 49 50 mass: 51 stellar mass in units of M_Sun 52 53 Returns 54 ------- 55 56 sfr: 57 the star formation rate in M_Sun / year. 58 59 """ 60 return 10**((numpy.log10(mass) - 8.7)/1.06)

61

62 -def sfr_from_L_nu(luminosity):

63 """Use Kennicutt (1998) conversion from UV luminosity to star formation rate. 64 65 Parameters 66 ---------- 67 68 luminosity: 69 the luminosity in units of ergs s^-1 Hz^-1 anywhere between 70 1500-2800 Angstroms. 71 72 Returns 73 ------- 74 75 The SFR in Msun/year. 76 77 Notes 78 ----- 79 80 Kennicutt (1998ARA&A..36..189K) says: 81 82 SFR/(MSun/year) = 1.4 * 10^-28 (L_nu/ergs s^-1 Hz^-1) 83 84 where L_nu is the UV luminosity anywhere between 1500-2800 Angstroms. 85 """ 86 return 1.4e-28 * luminosity

87

88 -def L_nu_from_sfr(sfr):

89 """Use Kennicutt (1998) conversion from UV luminosity to star formation rate. 90 91 Parameters 92 ---------- 93 94 sfr: 95 The SFR in Msun/year. 96 97 Returns 98 ------- 99 100 luminosity: 101 the luminosity in units of ergs s^-1 Hz^-1 anywhere between 102 1500-2800 Angstroms. 103 104 Notes 105 ----- 106 107 Kennicutt (1998ARA&A..36..189K) says: 108 109 SFR/(MSun/year) = 1.4 * 10^-28 (L_nu/ergs s^-1 Hz^-1) 110 111 where L_nu is the UV luminosity anywhere between 1500-2800 Angstroms. 112 """ 113 return sfr/(1.4e-28)

114

115 -def magnitudeAB_from_sfr(sfr):

116 """Use Kennicutt (1998) conversion from UV luminosity to AB magnitude. 117 118 Convenience function: uses L_nu_from_sfr and 119 magnitudes.magnitude_AB_from_L_nu. 120 """ 121 lnu = L_nu_from_sfr(sfr) 122 return magnitudes.magnitude_AB_from_L_nu(lnu)

123

124 -def schechterL(luminosity, phiStar, alpha, LStar):

125 """Schechter luminosity function.""" 126 LOverLStar = (luminosity/LStar) 127 return (phiStar/LStar) * LOverLStar**alpha * numpy.exp(- LOverLStar)

128

129 -def schechterM(magnitude, phiStar, alpha, MStar):

130 """Schechter luminosity function by magnitudes.""" 131 MStarMinM = 0.4 * (MStar - magnitude) 132 return (0.4 * numpy.log(10) * phiStar * 133 10.0**(MStarMinM * (alpha + 1.)) * numpy.exp(-10.**MStarMinM))

134

135 -def schechterCumuLL(luminosity, phiStar, alpha, LStar):

136 """Integrate luminosity in galaxies above luminosity=L. 137 138 Uses an analytical formula. 139 """ 140 # Note that the scipy.special definition of incomplete gamma is 141 # normalized to one and is the lower incomplete gamma function, so 142 # we have to use the complement and multiply by the unnormalized 143 # integral (which is computed in schechterTotLL). 144 return (schechterTotLL(phiStar, alpha, LStar) * 145 (1. - scipy.special.gammainc(alpha+2., luminosity/LStar)))

146

147 -def schechterCumuLM(magnitudeAB, phiStar, alpha, MStar):

148 """Integrate luminosity in galaxies brighter than magnitudeAB. 149 150 Uses an analytical formula. 151 """ 152 LStar = magnitudes.L_nu_from_magAB(MStar) 153 lum = magnitudes.L_nu_from_magAB(magnitudeAB) 154 return schechterCumuLL(lum, phiStar, alpha, LStar)

155

156 -def schechterTotLL(phiStar, alpha, LStar):

157 """Integrate total luminosity in galaxies. 158 159 Uses an analytical formula. 160 """ 161 return phiStar * LStar * scipy.special.gamma(alpha + 2.)

162

163 -def schechterTotLM(phiStar, alpha, MStar):

164 """Integrate total luminosity in galaxies. 165 166 Uses an analytical formula. 167 """ 168 LStar = magnitudes.L_nu_from_magAB(MStar) 169 return schechterTotLL(phiStar, alpha, LStar)

170

171 -def iPhotonRateDensity(schechterParams, 172 maglim=None, 173 sedParams={}, 174 wavelength=1500.):

175 """Ionizing photon rate density from a luminosity function. 176 177 in units of photons s^-1. 178 179 Given schecterParams, the parameters of a Schechter luminosity 180 function (in terms of AB Magnitudes), sedParams, the parameters of 181 the galactic Spectral Energy Distribution, and the wavelength of 182 the AB Magnitudes, calculate the emission rate density of ionizing 183 photons. 184 185 See Also 186 -------- 187 188 BrokenPowerlawSED 189 190 schechterTotLM 191 192 """ 193 if maglim is None: 194 lum = schechterTotLM(**schechterParams) 195 else: 196 lum = schechterCumuLM(maglim, **schechterParams) 197 rQL = BrokenPowerlawSED(**sedParams).iPhotonRateRatio(wavelength) 198 return lum * rQL

199 200 # From Bouwens et al. (2007) 201 B2007_z4 = {'MStar': -20.98, 202 'phiStar': 1.3e-3, 203 'alpha': -1.73, 204 'z':3.8} 205 206 B2007_z5 = {'MStar': -20.64, 207 'phiStar': 1.0e-3, 208 'alpha': -1.66, 209 'z':5.0} 210 211 B2007_z6 = {'MStar': -20.24, 212 'phiStar': 1.4e-3, 213 'alpha': -1.74, 214 'z':5.9} 215 216 B2007_z7 = {'MStar': -19.3, 217 'phiStar': 1.4e-3, 218 'alpha': -1.74, 219 'z':7.4 220 } 221 222 # From Oesch et al. (2010) 223 O2010_z7 = {'MStar': -19.9, 224 'phiStar': 1.4e-3, 225 'alpha': -1.77, 226 'z':6.8 # Not sure about this 227 } 228 229 # From Bouwens 2008ApJ...686..230B 230 B2008 = {'z' : [3.8, 5.0, 5.9, 7.3, 9.0], 231 'MStar': [-20.98, -20.64, -20.24, -19.8, -19.6], 232 'phiStar': [1.3e-3, 1.0e-3, 1.4e-3, 1.1e-3, 1.1e-3], 233 'alpha': [-1.73, -1.66, -1.74, -1.74, -1.74]} 234 235 # Like B2008, but with fixed phiStar, alpha 236 B2008_fixed = {'z' : [3.8, 5.0, 5.9, 7.3, 9.0], 237 'MStar': [-20.98, -20.64, -20.24, -19.8, -19.6], 238 'phiStar': [1.1e-3, 1.1e-3, 1.1e-3, 1.1e-3, 1.1e-3], 239 'alpha': [-1.74, -1.74, -1.74, -1.74, -1.74]} 240 241 # Like B2008, but with fixed phiStar, alpha 242 zlinear = numpy.array([3.8, 5.0, 5.9, 7.3, 9.0]) 243 B2008_linear_z5 = {'z' : zlinear, 244 'MStar': 0.36 * (zlinear - 5.0) - 20.64, 245 'phiStar': [1.1e-3, 1.1e-3, 1.1e-3, 1.1e-3, 1.1e-3], 246 'alpha': [-1.74, -1.74, -1.74, -1.74, -1.74]} 247 248 # From Trenti et al. 2010, EXCEPT the z=9 alpha value, which I lowered. 249 T2010ICLF = {'z' : [4., 5., 7., 8., 9.], 250 'MStar' : [-20.90, -20.55, -20.00, -19.70, -19.55], 251 'phiStar' : [1.3e-3, 1.5e-3, 1.0e-3, 0.60e-3, 0.22e-3], 252 'alpha' : [-1.57, -1.63, -1.84, -1.90, -1.999]} 253

254 -class LFHistory(Saveable):

255 """Interpolate / extrapolate the Schechter parameters. 256 257 By default, above the observed redshift range: 258 259 MStar is linearly extrapolated as a function of time (not z) to high z. 260 261 phiStar and alpha are constant at high z. 262 263 """ 264

265 - def __init__(self, params=B2008, 266 MStar_bounds = ['extrapolate', float('NaN')], 267 phiStar_bounds = ['constant', float('NaN')], 268 alpha_bounds = ['constant', float('NaN')], 269 extrap_args = {}, 270 extrap_var = 'z', 271 sedParams = {}, 272 wavelength = 1500., 273 **cosmo):

274 275 for (k, v) in params.iteritems(): 276 params[k] = numpy.asarray(v) 277 278 self.params = params 279 self.MStar_bounds = MStar_bounds 280 self.phiStar_bounds = phiStar_bounds 281 self.alpha_bounds = alpha_bounds 282 self.extrap_args = extrap_args 283 self.extrap_var = extrap_var 284 self.sedParams = sedParams 285 self.wavelength = wavelength 286 self.cosmo = cosmo 287 288 self.zobs = params['z'] 289 self.tobs = cd.age(self.zobs, **cosmo)[0] 290 self.MStar = params['MStar'] 291 self.phiStar = params['phiStar'] 292 self.alpha = params['alpha'] 293 294 if extrap_var == 't': 295 self.xobs = self.tobs 296 self._iPhotFunc = \ 297 lambda t1, mag: (self.iPhotonRateDensity_t(t1, 298 maglim=mag)) 299 300 elif extrap_var == 'z': 301 self.xobs = self.zobs 302 MStar_bounds = MStar_bounds[::-1] 303 phiStar_bounds = phiStar_bounds[::-1] 304 alpha_bounds = alpha_bounds[::-1] 305 self._iPhotFunc = \ 306 lambda z1, mag: (self.iPhotonRateDensity_z(z1, 307 maglim=mag)) 308 309 self._MStarfunc = utils.Extrapolate1d(self.xobs, self.MStar, 310 bounds_behavior=MStar_bounds, 311 **extrap_args 312 ) 313 print "M*:", 314 print self._MStarfunc.extrap_string() 315 316 self._phiStarfunc = utils.Extrapolate1d(self.xobs, self.phiStar, 317 bounds_behavior=phiStar_bounds, 318 **extrap_args 319 ) 320 print "phi*:", 321 print self._phiStarfunc.extrap_string() 322 self._alphafunc = utils.Extrapolate1d(self.xobs, self.alpha, 323 bounds_behavior=alpha_bounds, 324 **extrap_args 325 ) 326 print "alpha:", 327 print self._alphafunc.extrap_string() 328 329 self._SED = BrokenPowerlawSED(**sedParams) 330 self._rQL = self._SED.iPhotonRateRatio(wavelength) 331 332 for name, func in globals().items(): 333 if not name.startswith('schechter'): 334 continue 335 def newfunc(z, _name=name, _func=func, **kwargs): 336 params = self.params_z(z) 337 M = params['MStar'] 338 phi = params['phiStar'] 339 alpha = params['alpha'] 340 return _func(MStar=M, phiStar=phi, alpha=alpha, **kwargs)

341 newfunc.__name__ = func.__name__ 342 newfunc.__doc__ = func.__doc__ 343 self.__dict__[name] = newfunc

344

345 - def iPhotonRateDensity_z(self, z, maglim=None):

346 """Ionizing photon rate density from a luminosity function. 347 348 See the iPhotonRateRatio function. 349 """ 350 params = self.params_z(z) 351 if maglim is None: 352 lum = schechterTotLM(**params) 353 else: 354 lum = schechterCumuLM(maglim, **params) 355 return lum * self._rQL

356

357 - def iPhotonRateDensity_t(self, t, maglim=None):

358 """Ionizing photon rate density from a luminosity function. 359 360 See the iPhotonRateRatio function. 361 """ 362 params = self.params_t(t) 363 if maglim is None: 364 lum = schechterTotLM(**params) 365 else: 366 lum = schechterCumuLM(maglim, **params) 367 return lum * self._rQL

368

369 - def ionization(self, z, maglim=None):

370 xH = cr.ionization_from_luminosity(z, 371 self._iPhotFunc, 372 rate_is_tfunc = 373 self.extrap_var == 't', 374 ratedensityfunc_args=[maglim], 375 **self.cosmo) 376 return xH

377

378 - def params_t(self, t):

379 """Return interp/extrapolated Schechter function parameters.""" 380 if self.extrap_var == 't': 381 return {'MStar':self._MStarfunc(t), 382 'phiStar':self._phiStarfunc(t), 383 'alpha':self._alphafunc(t)} 384 elif self.extrap_var == 'z': 385 ### FIX THIS ### 386 raise NotImplementedError, "params_t not implemented for z interps!"

387

388 - def params_z(self, z):

389 """Return interp/extrapolated Schechter function parameters.""" 390 if self.extrap_var == 'z': 391 return {'MStar':self._MStarfunc(z), 392 'phiStar':self._phiStarfunc(z), 393 'alpha':self._alphafunc(z)} 394 elif self.extrap_var == 't': 395 z = numpy.atleast_1d(z) 396 t = cd.age(z, **self.cosmo)[0] 397 return self.params_t(t)

398

399 -def plotLFevo(hist=None, 400 params=B2008, 401 extrap_var='t', 402 #maglim = -21.07 - 2.5 * numpy.log10(0.2)): 403 maglim = -21. - 2.5 * numpy.log10(0.2), 404 z_max=20.0, 405 skipIon=True 406 ):

407 408 """Plot evolution of luminosity function params and total luminsity. 409 410 Schechter function at each redshift is integrated up to maglim to 411 find total luminsity. 412 """ 413 414 for (k, v) in params.iteritems(): 415 params[k] = numpy.asarray(v) 416 417 if hist is None: 418 cosmo = cp.WMAP7_BAO_H0_mean(flat=True) 419 hist = LFHistory(params, extrap_var=extrap_var, **cosmo) 420 else: 421 params = hist.params 422 extrap_var = hist.extrap_var 423 cosmo = hist.cosmo 424 425 z = params['z'] 426 t = cd.age(z, **cosmo)[0] / cc.yr_s 427 MStar = params['MStar'] 428 phiStar = params['phiStar'] 429 alpha = params['alpha'] 430 431 if maglim is None: 432 ltot = schechterTotLM(MStar=MStar, phiStar=phiStar, alpha=alpha) 433 else: 434 ltot = schechterCumuLM(magnitudeAB=maglim, 435 MStar=MStar, phiStar=phiStar, alpha=alpha) 436 437 print hist._MStarfunc.extrap_string() 438 439 zPlot = numpy.arange(z.min()-0.1, z_max, 0.1) 440 tPlot = cd.age(zPlot, **cosmo)[0] / cc.yr_s 441 newparams = hist.params_z(zPlot) 442 MPlot = newparams['MStar'] 443 phiPlot = newparams['phiStar'] 444 alphaPlot = newparams['alpha'] 445 446 if maglim is None: 447 ltotPlot = hist.schechterTotLM(zPlot) 448 else: 449 ltotPlot = hist.schechterCumuLM(zPlot, magnitudeAB=maglim) 450 # From Table 6 of 2008ApJ...686..230B 451 lB2008 =10.** numpy.array([26.18, 25.85, 25.72, 25.32, 25.14]) 452 453 iPhot = hist.iPhotonRateDensity_z(zPlot, maglim=maglim) 454 # iPhotFunc = lambda t1: cc.yr_s * hist.iPhotonRateDensity_t(t1, 455 # maglim=maglim) 456 457 if not skipIon: 458 xH = hist.ionization(zPlot, maglim) 459 import pylab 460 pylab.figure(1) 461 pylab.gcf().set_label('LFion_vs_z') 462 if skipIon: 463 pylab.subplot(111) 464 else: 465 pylab.subplot(211) 466 pylab.plot(zPlot, iPhot) 467 if not skipIon: 468 pylab.subplot(212) 469 pylab.plot(zPlot, xH) 470 pylab.ylim(0.0, 1.5) 471 472 pylab.figure(2) 473 pylab.gcf().set_label('LFparams_vs_z') 474 pylab.subplot(311) 475 pylab.plot(z, MStar, '-') 476 pylab.plot(z, MStar, 'o') 477 pylab.plot(zPlot, MPlot, ':') 478 479 pylab.subplot(312) 480 pylab.plot(z, phiStar, '-') 481 pylab.plot(z, phiStar, 'o') 482 pylab.plot(zPlot, phiPlot, ':') 483 484 pylab.subplot(313) 485 pylab.plot(z, alpha, '-') 486 pylab.plot(z, alpha, 'o') 487 pylab.plot(zPlot, alphaPlot, ':') 488 489 pylab.figure(3) 490 pylab.gcf().set_label('LFparams_vs_t') 491 pylab.subplot(311) 492 pylab.plot(t, MStar, '-') 493 pylab.plot(t, MStar, 'o') 494 pylab.plot(tPlot, MPlot, ':') 495 496 pylab.subplot(312) 497 pylab.plot(t, phiStar, '-') 498 pylab.plot(t, phiStar, '.') 499 pylab.plot(tPlot, phiPlot, ':') 500 501 pylab.subplot(313) 502 pylab.plot(t, alpha, '-') 503 pylab.plot(t, alpha, 'o') 504 pylab.plot(tPlot, alphaPlot, ':') 505 506 pylab.figure(4) 507 pylab.gcf().set_label('LFlum_vs_z') 508 pylab.subplot(121) 509 pylab.plot(z, ltot, 'o') 510 pylab.plot(z, lB2008, 'x') 511 pylab.plot(zPlot, ltotPlot) 512 513 pylab.subplot(122) 514 pylab.plot(t, ltot, 'o') 515 pylab.plot(t, lB2008, 'x') 516 pylab.plot(tPlot, ltotPlot)

517

518 -def test_plot_schechter():

519 520 phiStar = 1.8e-3 521 MStar = -20.04 522 alpha = -1.71 523 524 LStar = magnitudes.L_nu_from_magAB(MStar) 525 526 mags = numpy.arange(-22, -11.0, 0.5) 527 lums = magnitudes.L_nu_from_magAB(mags) 528 529 phi_L = schechterL(lums, phiStar, alpha, LStar) 530 phi_M = schechterM(mags, phiStar, alpha, MStar) 531 532 L_L = schechterCumuLL(lums, phiStar, alpha, LStar) 533 L_M = schechterCumuLM(mags, phiStar, alpha, MStar) 534 535 phi_L_func = lambda l: l * schechterL(l, phiStar, alpha, LStar) 536 L_L_num = utils.logquad(phi_L_func, lums, 1e35)[0] 537 L_L_num2 = utils.vecquad(phi_L_func, lums, 1e29)[0] 538 539 phi_M_func = lambda m: (magnitudes.L_nu_from_magAB(m) * 540 schechterM(m, phiStar, alpha, MStar)) 541 L_M_num2 = utils.vecquad(phi_M_func, -25, mags)[0] 542 543 Ltot_L = schechterTotLL(phiStar, alpha, LStar) 544 Ltot_M = schechterTotLM(phiStar, alpha, MStar) 545 546 pylab.figure() 547 pylab.subplot(221) 548 pylab.plot(lums, lums * lums * phi_L) 549 pylab.xscale('log') 550 pylab.yscale('log') 551 pylab.ylabel(r'$ L^2 \Phi_L$') 552 553 pylab.subplot(222) 554 pylab.plot(mags, -mags * lums * phi_M) 555 pylab.yscale('log') 556 pylab.ylabel(r'$ -M L \Phi_M$') 557 558 pylab.subplot(223) 559 pylab.plot(lums, Ltot_L - L_L) 560 pylab.plot(lums, L_M) 561 pylab.plot(lums, L_L_num, '--') 562 pylab.plot(lums, L_L_num2, ':') 563 pylab.plot(lums, L_M_num2, 'x') 564 pylab.axhline(y=Ltot_L) 565 pylab.axhline(y=Ltot_M) 566 pylab.xscale('log') 567 pylab.yscale('log') 568 569 pylab.subplot(224) 570 pylab.plot(mags, Ltot_M - L_M) 571 pylab.plot(mags, L_L) 572 pylab.plot(mags, L_L_num, '--') 573 pylab.plot(mags, L_L_num2, ':') 574 pylab.plot(mags, L_M_num2, 'x') 575 pylab.axhline(y=Ltot_L) 576 pylab.axhline(y=Ltot_M) 577 pylab.yscale('log')

578 579 #pylab.show() 580

581 -class BrokenPowerlawSED(object):

582 """Define an SED with a break at 912 Angstroms and different 583 slopes above and below the break. 584 """ 585 586 # From http://www.astro.wisc.edu/~dolan/constants.html 587 planck_erg_s = 6.6260755e-27 #erg s 588

589 - def lambdanu(self, wavelength):

590 """Convert between wavelength (Ang) and frequency (Hz) or vice-versa. 591 """ 592 if numpy.all(wavelength == 0): 593 wavelength = numpy.array(wavelength) 594 return cc.c_light_cm_s / (wavelength * cc.angstrom_cm)

595 596

597 - def luminosity_wavelength(wavelength0, wavelength1):

598 """Normalized luminosity in the defined band. 599 600 wavelength0 and wavelength1 in Angstroms. 601 602 This is the fraction of total luminosity in the band. 603 Multiply the result by the total luminosity (energy per 604 second) to get physical units. 605 """ 606 return self.sed.integrate(self.lambdanu(wavelength0), 607 self.lambdanu(wavelength1))

608

609 - def photonRate_wavelength(self, wavelength0, wavelength1):

610 """Normalized photon emission rate in the band between 611 wavelength0 and wavelength1. 612 613 Units are erg^-1 (which could also be expressed as s^-1 per 614 (erg/s)). 615 616 Multiply the result by the total luminosity (in ergs per unit 617 time), to get the actual photon emission rate. 618 619 Example 620 ------- 621 622 To get the ionizing photon emission rate (per unit luminosity): 623 >>> BrokenPowerlawSED().photonRate_wavelength(0., 912.) 624 3272819078.0292048 625 """ 626 return ((1./self.planck_erg_s) * 627 self.sed.integrate(self.lambdanu(wavelength0), 628 self.lambdanu(wavelength1), 629 weight_power=-1.))

630

631 - def __init__(self, s_ion=-3.0, s_red=0.0, break_factor=6.0):

632 """Return a model SED for a galaxy. 633 634 Parameters 635 ---------- 636 637 s_ion: 638 spectral index (f_nu ~ nu^s_ion) at lambda < 912 Ang 639 640 s_red: 641 spectral index (f_nu ~ nu^s_red) at lambda > 912 Ang 642 643 Notes 644 ----- 645 646 Bolton and Haehnelt (2007MNRAS.382..325B) use an SED with 647 648 eps_nu ~ v^0 for 912 < lambda < 3000 Ang. 649 ~ v^-3 for labmda < 912 Ang. 650 651 'with an additional break in the spectrum at the Lyman limit' 652 653 eps_L = eps(1500)/6 654 655 """ 656 657 limits_Ang=numpy.array((3000., 912., 0.)) 658 limits_nu = self.lambdanu(limits_Ang) 659 660 nu_912 = self.lambdanu(912.) 661 662 # Ratio of flux at 1500 Ang to flux at 912 Ang (redward of the break). 663 ratio_1500_912 = (912./1500.)**s_red 664 665 # Ratio of the ionizing to red coefficients. 666 coeff_ratio = ratio_1500_912 * nu_912**(s_red - s_ion) / break_factor 667 668 coefficients = numpy.array((1., coeff_ratio)) 669 670 # Define an SED: 671 powers = numpy.array((s_red, s_ion)) 672 sed = utils.PiecewisePowerlaw(limits_nu, powers, coefficients) 673 self.sed = sed 674 self.break_factor = break_factor 675 self.s_ion = s_ion 676 self.s_red = s_red

677

678 - def __call__(self, nu):

679 """The spectrum at the given frequency/frequencies. 680 681 Multiply by the total luminosity to get the luminosity per 682 unit frequency. 683 684 Units are (fraction of total luminosity) per Hz. 685 """ 686 return self.sed(nu)

687

688 - def iPhotonRateRatio(self, wavelength=1500.):

689 """The ratio of ionizing photon emission rate to luminosity. 690 691 Paramters 692 --------- 693 694 lambda: 695 wavelength in Angstroms 696 697 Returns 698 ------- 699 700 The ratio of ionizing photon emission rate to luminosity at 701 the given wavelength Q/L_nu(lambda) in units of photons s^-1 702 (erg s^-1 Hz^-1)^-1. 703 704 Notes 705 ----- 706 707 While this function takes an argument in angstroms, the ratio 708 is calculated using the luminosity per unit frequence, so as 709 to be easily comensurate with luminosities inferred from AB 710 magnitudes. 711 712 """ 713 714 return (self.photonRate_wavelength(0.,912.) / 715 self(self.lambdanu(wavelength)))

716 717 if __name__ == '__main__': 718 719 usage = """Run with a filename argument to produce image files, e.g.: 720 python luminosityfunction.py lum.png 721 """ 722 723 ### Argument parsing. ### 724 parser = optparse.OptionParser(usage) 725 parser.add_option("-f", "--file", action='store_true', dest="filename", 726 default=None) 727 (options, args) = parser.parse_args() 728 if (len(args) > 0) and (options.filename is None): 729 options.filename = args[0] 730 731 if options.filename is None: 732 print "No filename given." 733 print usage 734 else: 735 prefix, extension = os.path.splitext(options.filename) 736 737 import pylab 738 739 plotLFevo() # Canonical observed LF 740 plotLFevo(extrap_var='z') 741 plotLFevo(params=B2008_fixed) 742 plotLFevo(params=B2008_fixed, extrap_var='z') 743 744 test_plot_schechter() 745 746 # plotLFevo(maglim=None) # Integrate LF all the way to zero luminosity 747 # plotLFevo(B2008_fixed) # Use fixed alpha, phiStar 748 # plotLFevo(B2008_fixed, maglim=None) # ... integrated to L=0 749 750 # ... with alpha raised to -1.4 (from -1.74) 751 # newparams = B2008_fixed 752 # newparams['alpha'] = numpy.array(newparams['alpha']) 753 # newparams['alpha'][:] = -1.4 754 755 # plotLFevo(newparams, maglim=None) 756 757 # test_plot_schechter() 758 759 ### Plot output code. ### 760 if options.filename is None: 761 pylab.show() 762 else: 763 from matplotlib import _pylab_helpers 764 for manager in _pylab_helpers.Gcf.get_all_fig_managers(): 765 fig = manager.canvas.figure 766 if len(fig.get_label()) > 0: 767 label = fig.get_label() 768 else: 769 label = '_Fig' + str(manager.num) 770 newfilename = prefix + '_' + label + extension 771 fig.savefig(newfilename, dpi=75)#, bbox_inches="tight") 772

Source Code for Module cosmolopy.luminosityfunction