Package cosmolopy ::
Module 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
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
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
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
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
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
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
155
157 """Integrate total luminosity in galaxies.
158
159 Uses an analytical formula.
160 """
161 return phiStar * LStar * scipy.special.gamma(alpha + 2.)
162
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)