cosmolopy.perturbation

1 """Perturbation theory and the power spectrum routines. 2 3 This module relies largely on power.c from Eisenstein & Hu (1999 ApJ 511 5) 4 5 http://background.uchicago.edu/~whu/transfer/transferpage.html 6 7 See notes in `cosmolopy.EH`. 8 9 """ 10 11 import math 12 import warnings 13 14 import numpy 15 import scipy 16 import scipy.special 17 import scipy.integrate as si 18 19 import constants as cc 20 import density as cden 21 22 powererror = None 23 try: 24 import EH.power as power 25 havepower = True 26 except ImportError as ie: 27 havepower = False 28 powererror = ie 29 pass 30 31 tffiterror = None 32 try: 33 import EH.tf_fit as tf_fit 34 havetffit = True 35 except ImportError as ie: 36 havetffit = False 37 tffiterror = ie 38 pass

39 40 # Turn on printing of special function error messages. 41 #scipy.special.errprint(1) 42 43 @numpy.vectorize 44 -def _vec_transfer_func(k,baryonic_effects=False):

45 if baryonic_effects: 46 return (tf_fit.TFfit_onek(k), power.TFmdm_onek_mpc(k)) 47 else: 48 return (power.TFmdm_onek_mpc(k), power.cvar.tf_cbnu)

49 -def transfer_function_EH(k, **cosmology):

50 """The transfer function as a function of wavenumber k. 51 52 Parameters 53 ---------- 54 55 cosmology : dict 56 Specify the cosmological parameters with the keys 'omega_M_0', 57 'omega_b_0', 'omega_n_0', 'N_nu', 'omega_lambda_0', 'h' and 58 'baryonic_effects'. 59 60 k : array 61 Wavenumber in Mpc^-1. 62 63 Returns 64 ------- 65 66 If baryonic_effects is true, returns a tuple of arrays matching 67 the shape of k: 68 69 (the transfer function for CDM + Baryons with baryonic effects, 70 the transfer function for CDM + Baryons without baryonic effects) 71 72 Otherwise, returns a tuple of arrays matching the shape of k: 73 74 (the transfer function for CDM + Baryons, 75 the transfer function for CDM + Baryons + Neutrinos). 76 77 Notes 78 ----- 79 80 Uses transfer function code power.c from Eisenstein & Hu (1999 ApJ 511 5). 81 For baryonic effects, uses tf_fit.c from Eisenstein & Hu (1997 ApJ 496 605). 82 83 http://background.uchicago.edu/~whu/transfer/transferpage.html 84 85 """ 86 baryonic_effects = cosmology['baryonic_effects'] 87 if baryonic_effects: 88 if not havetffit: 89 raise ImportError, "Could not import EH.tf_fit module. Transfer function cannot be calculated." 90 if not havepower: 91 raise ImportError, "Could not import EH.power module. Transfer function cannot be calculated." 92 else: 93 if not havepower: 94 raise ImportError, "Could not import EH.power module. Transfer function cannot be calculated." 95 96 z_val=0 97 98 if not baryonic_effects: 99 # Baryonic effect are not used. Default to more general CDM variants transfer function. 100 # 101 # /* TFmdm_set_cosm() -- User passes all the cosmological parameters as 102 # arguments; the routine sets up all of the scalar quantites needed 103 # computation of the fitting formula. The input parameters are: 104 # 1) omega_matter -- Density of CDM, baryons, and massive neutrinos, 105 # in units of the critical density. 106 # 2) omega_baryon -- Density of baryons, in units of critical. 107 # 3) omega_hdm -- Density of massive neutrinos, in units of critical 108 # 4) degen_hdm -- (Int) Number of degenerate massive neutrino species 109 # 5) omega_lambda -- Cosmological constant 110 # 6) hubble -- Hubble constant, in units of 100 km/s/Mpc 111 # 7) redshift -- The redshift at which to evaluate */ 112 if int(cosmology['N_nu']) != cosmology['N_nu']: 113 raise TypeError('N_nu must be an integer.') 114 power.TFmdm_set_cosm(cosmology['omega_M_0'], cosmology['omega_b_0'], 115 cosmology['omega_n_0'], int(cosmology['N_nu']), 116 cosmology['omega_lambda_0'], cosmology['h'], 117 z_val) 118 119 # Given a wavenumber in Mpc^-1, return the transfer function for 120 # the cosmology held in the global variables. 121 if numpy.isscalar(k): 122 return (power.TFmdm_onek_mpc(k), power.cvar.tf_cbnu) 123 else: 124 return _vec_transfer_func(k) 125 else: 126 # Baryonic effects are in use. This reduces the range of validity of the 127 # transfer function, for instance by not including effects from neutrinos. 128 # 129 # /* TFset_parameters() -- User passes certain cosmological parameters are 130 # arguments; the routine sets up all of the scalar quantities needed 131 # in the computation of the fitting formula. The input parameters are: 132 # 1) omega_mhh -- Density of (CDM+Baryons), in units of the critical 133 # density, times the hubble parameter squared. 134 # 2) f_baryon -- The fraction of baryons in all matter 135 # 3) Tcmb -- the CMB temperature in Kelvin. Set to 2.728. 136 omhh = cosmology['omega_M_0'] * cosmology['h'] * cosmology['h'] 137 fbaryon = cosmology['omega_b_0'] / cosmology['omega_M_0'] 138 Tcmb = 2.728 139 tf_fit.TFset_parameters(omhh, fbaryon, Tcmb) 140 141 # Given a wavenumber in Mpc^-1, return the transfer function for 142 # the cosmology held in the global variables. 143 if numpy.isscalar(k): 144 return (tf_fit.TFfit_onek(k), power.TFmdm_onek_mpc(k)) 145 else: 146 return _vec_transfer_func(k, baryonic_effects)

147

148 149 -def fgrowth(z, omega_M_0, unnormed=False):

150 r"""Cosmological perturbation growth factor, normalized to 1 at z = 0. 151 152 Approximate forumla from Carol, Press, & Turner (1992, ARA&A, 30, 153 499), "good to a few percent in regions of plausible Omega_M, 154 Omega_Lambda". 155 156 This is proportional to D_1(z) from Eisenstein & Hu (1999 ApJ 511 157 5) equation 10, but the normalization is different: fgrowth = 1 at 158 z = 0 and ``D_1(z) = \frac{1+z_\mathrm{eq}}{1+z}`` as z goes 159 to infinity. 160 161 To get D_1 one would just use 162 163 :: 164 165 D_1(z) = (1+z_\mathrm{eq}) \mathtt{fgrowth}(z,\Omega_{M0}, 1) 166 167 (see \EH\ equation 1 for z_eq). 168 169 :: 170 171 \mathtt{fgrowth} = \frac{D_1(z)}{D_1(0)} 172 173 Setting unnormed to true turns off normalization. 174 175 Note: assumes Omega_lambda_0 = 1 - Omega_M_0! 176 177 """ 178 #if cden.get_omega_k_0(**) != 0: 179 # raise ValueError, "Not valid for non-flat (omega_k_0 !=0) cosmology." 180 181 omega = cden.omega_M_z(z, omega_M_0=omega_M_0, omega_lambda_0=1.-omega_M_0) 182 lamb = 1 - omega 183 a = 1/(1 + z) 184 185 if unnormed: 186 norm = 1.0 187 else: 188 norm = 1.0 / fgrowth(0.0, omega_M_0, unnormed=True) 189 return (norm * (5./2.) * a * omega / 190 (omega**(4./7.) - lamb + (1. + omega/2.) * (1. + lamb/70.)) 191 )

192

193 -def w_tophat(k, r):

194 r"""The k-space Fourier transform of a spherical tophat. 195 196 Parameters 197 ---------- 198 199 k: array 200 wavenumber 201 202 r: array 203 radius of the 3-D spherical tophat 204 205 Note: k and r need to be in the same units. 206 207 Returns 208 ------- 209 210 ``\tilde{w}``: array 211 the value of the transformed function at wavenumber k. 212 213 """ 214 return (3. * ( numpy.sin(k * r) - k * r * numpy.cos(k * r) ) / 215 ((k * r)**3.))

216

217 -def w_gauss(k, r):

218 r"""The k-space Fourier transform of an isotropic three-dimensional gaussian 219 220 Parameters 221 ---------- 222 223 k: array 224 wavenumber 225 226 r: array 227 width of the 3-D gaussian 228 229 Note: k and r need to be in the same units. 230 231 Returns 232 ------- 233 234 ``\tilde{w}``: array 235 the value of the transformed function at wavenumber k. 236 237 """ 238 return numpy.exp( -(k * r)**2./2. )

239

240 -def _sigmajsq_integrand_log(logk, r, j, cosmology):

241 """Integrand used internally by the sigma_j function. 242 """ 243 k = numpy.exp(logk) 244 # The 1e-10 factor in the integrand is added to avoid roundoff 245 # error warnings. It is divided out later. 246 return (k * 247 (1.e-10 / (2. * math.pi**2.)) * k**(2.*(j+1.)) * 248 w_gauss(k, r)**2. * 249 power_spectrum(k, 0.0, **cosmology))

250

251 -def _sigmasq_integrand_log(logk, r, cosmology):

252 """Integrand used internally by the sigma_r function. 253 """ 254 k = numpy.exp(logk) 255 # The 1e-10 factor in the integrand is added to avoid roundoff 256 # error warnings. It is divided out later. 257 return (k * 258 (1.e-10 / (2. * math.pi**2.)) * k**2. * 259 w_tophat(k, r)**2. * 260 power_spectrum(k, 0.0, **cosmology))

261

262 -def _klimsj(r, j, cosmology):

263 """Integration limits used internally by the sigma_j function.""" 264 logk = numpy.arange(-20., 20., 0.1) 265 integrand = _sigmajsq_integrand_log(logk, r, j, cosmology) 266 267 maxintegrand = numpy.max(integrand) 268 factor = 1.e-4 269 highmask = integrand > maxintegrand * factor 270 while highmask.ndim > logk.ndim: 271 highmask = numpy.logical_or.reduce(highmask) 272 273 mink = numpy.min(logk[highmask]) 274 maxk = numpy.max(logk[highmask]) 275 276 return mink, maxk

277

278 -def _klims(r, cosmology):

279 """Integration limits used internally by the sigma_r function.""" 280 logk = numpy.arange(-20., 20., 0.1) 281 integrand = _sigmasq_integrand_log(logk, r, cosmology) 282 283 maxintegrand = numpy.max(integrand) 284 factor = 1.e-4 285 highmask = integrand > maxintegrand * factor 286 while highmask.ndim > logk.ndim: 287 highmask = numpy.logical_or.reduce(highmask) 288 289 mink = numpy.min(logk[highmask]) 290 maxk = numpy.max(logk[highmask]) 291 292 return mink, maxk

293

294 -def _sigmasq_r_scalar(r, 295 n, deltaSqr, omega_M_0, omega_b_0, omega_n_0, N_nu, 296 omega_lambda_0, h, baryonic_effects):

297 """sigma_r^2 at z=0. Works only for scalar r. 298 299 Used internally by the sigma_r function. 300 301 Parameters 302 ---------- 303 304 r : array 305 radius in Mpc. 306 307 n, omega_M_0, omega_b_0, omega_n_0, N_nu, omega_lambda_0, h, baryonic_effecs: 308 cosmological parameters, specified like this to allow this 309 function to be vectorized (see source code of sigma_r). 310 311 Returns 312 ------- 313 314 sigma^2, error(sigma^2) 315 316 """ 317 # r is in Mpc, so k will also by in Mpc for the integration. 318 319 cosmology = {'n':n, 320 'deltaSqr':deltaSqr, 321 'omega_M_0':omega_M_0, 322 'omega_b_0':omega_b_0, 323 'omega_n_0':omega_n_0, 324 'N_nu':N_nu, 325 'omega_lambda_0':omega_lambda_0, 326 'h':h, 327 'baryonic_effects':baryonic_effects} 328 329 logk_lim = _klims(r, cosmology) 330 #print "Integrating from logk = %.1f to %.1f." % logk_lim 331 332 # Integrate over logk from -infinity to infinity. 333 integral, error = si.quad(_sigmasq_integrand_log, 334 logk_lim[0], 335 logk_lim[1], 336 args=(r, cosmology), 337 limit=10000)#, epsabs=1e-9, epsrel=1e-9) 338 return 1.e10 * integral, 1.e10 * error

339 340 _sigmasq_r_vec = numpy.vectorize(_sigmasq_r_scalar)

341 342 -def _sigmasq_j_scalar(r, j, 343 n, deltaSqr, omega_M_0, omega_b_0, omega_n_0, N_nu, 344 omega_lambda_0, h, baryonic_effects):

345 """sigma_j^2(r) at z=0. Works only for scalar r. 346 347 Used internally by the sigma_j function. 348 349 Parameters 350 ---------- 351 352 r : array 353 radius in Mpc. 354 355 j : array 356 order of sigma statistic. 357 358 n, omega_M_0, omega_b_0, omega_n_0, N_nu, omega_lambda_0, h: 359 cosmological parameters, specified like this to allow this 360 function to be vectorized (see source code of sigma_r). 361 362 Returns 363 ------- 364 365 sigma^2, error(sigma^2) 366 367 """ 368 # r is in Mpc, so k will also by in Mpc for the integration. 369 370 cosmology = {'n':n, 371 'deltaSqr':deltaSqr, 372 'omega_M_0':omega_M_0, 373 'omega_b_0':omega_b_0, 374 'omega_n_0':omega_n_0, 375 'N_nu':N_nu, 376 'omega_lambda_0':omega_lambda_0, 377 'h':h, 378 'baryonic_effects':baryonic_effects} 379 logk_lim = _klimsj(r, j, cosmology) 380 #print "Integrating from logk = %.1f to %.1f." % logk_lim 381 382 # Integrate over logk from -infinity to infinity. 383 integral, error = si.quad(_sigmajsq_integrand_log, 384 logk_lim[0], 385 logk_lim[1], 386 args=(r, j, cosmology), 387 limit=10000)#, epsabs=1e-9, epsrel=1e-9) 388 return 1.e10 * integral, 1.e10 * error

389 390 _sigmasq_j_vec = numpy.vectorize(_sigmasq_j_scalar)

391 392 -def sigma_j(r, j, z, **cosmology):

393 r"""Sigma statistic of order j for gaussian field of variancea r at redshift z. 394 395 Returns sigma and the error on sigma. 396 397 Parameters 398 ---------- 399 400 r : array 401 radius of sphere in Mpc 402 403 j : array 404 order of the sigma statistic (0, 1, 2, 3, ...) 405 406 z : array 407 redshift 408 409 Returns 410 ------- 411 412 sigma: 413 j-th order variance of the field smoothed by gaussian with with r 414 415 error: 416 An estimate of the numerical error on the calculated value of sigma. 417 418 Notes 419 ----- 420 :: Eq. (152) of Matsubara (2003) 421 422 \sigma_j(R,z) = \sqrt{\int_0^\infty \frac{k^2}{2 \pi^2}~P(k, z)~k^{2j} 423 \tilde{w}_k^2(k, R)~dk} = \sigma_j(R,0) \left(\frac{D_1(z)}{D_1(0)}\right) 424 425 """ 426 omega_M_0 = cosmology['omega_M_0'] 427 428 fg = fgrowth(z, omega_M_0) 429 430 if 'deltaSqr' not in cosmology: 431 cosmology['deltaSqr'] = norm_power(**cosmology) 432 433 #Uses 'n', as well as (for transfer_function_EH), 'omega_M_0', 434 #'omega_b_0', 'omega_n_0', 'N_nu', 'omega_lambda_0', and 'h'. 435 436 if numpy.isscalar(r): 437 sigmasq_0, errorsq_0 = _sigmasq_j_scalar(r, j, 438 cosmology['n'], 439 cosmology['deltaSqr'], 440 cosmology['omega_M_0'], 441 cosmology['omega_b_0'], 442 cosmology['omega_n_0'], 443 cosmology['N_nu'], 444 cosmology['omega_lambda_0'], 445 cosmology['h'], 446 cosmology['baryonic_effects'],) 447 else: 448 sigmasq_0, errorsq_0 = _sigmasq_j_vec(r, j, 449 cosmology['n'], 450 cosmology['deltaSqr'], 451 cosmology['omega_M_0'], 452 cosmology['omega_b_0'], 453 cosmology['omega_n_0'], 454 cosmology['N_nu'], 455 cosmology['omega_lambda_0'], 456 cosmology['h'], 457 cosmology['baryonic_effects'],) 458 sigma = numpy.sqrt(sigmasq_0) * fg 459 460 # Propagate the error on sigmasq_0 to sigma. 461 error = fg * errorsq_0 / (2. * sigmasq_0) 462 463 return sigma, error

464

465 -def sigma_r(r, z, **cosmology):

466 r"""RMS mass fluctuations of a sphere of radius r at redshift z. 467 468 Returns sigma and the error on sigma. 469 470 Parameters 471 ---------- 472 473 r : array 474 radius of sphere in Mpc 475 476 z : array 477 redshift 478 479 Returns 480 ------- 481 482 sigma: 483 RMS mass fluctuations of a sphere of radius r at redshift z. 484 485 error: 486 An estimate of the numerical error on the calculated value of sigma. 487 488 Notes 489 ----- 490 :: 491 492 \sigma(R,z) = \sqrt{\int_0^\infty \frac{k^2}{2 \pi^2}~P(k, z)~ 493 \tilde{w}_k^2(k, R)~dk} = \sigma(R,0) \left(\frac{D_1(z)}{D_1(0)}\right) 494 495 """ 496 omega_M_0 = cosmology['omega_M_0'] 497 498 fg = fgrowth(z, omega_M_0) 499 500 if 'deltaSqr' not in cosmology: 501 cosmology['deltaSqr'] = norm_power(**cosmology) 502 503 #Uses 'n', as well as (for transfer_function_EH), 'omega_M_0', 504 #'omega_b_0', 'omega_n_0', 'N_nu', 'omega_lambda_0', and 'h'. 505 506 if numpy.isscalar(r): 507 sigmasq_0, errorsq_0 = _sigmasq_r_scalar(r, 508 cosmology['n'], 509 cosmology['deltaSqr'], 510 cosmology['omega_M_0'], 511 cosmology['omega_b_0'], 512 cosmology['omega_n_0'], 513 cosmology['N_nu'], 514 cosmology['omega_lambda_0'], 515 cosmology['h'], 516 cosmology['baryonic_effects'],) 517 else: 518 sigmasq_0, errorsq_0 = _sigmasq_r_vec(r, 519 cosmology['n'], 520 cosmology['deltaSqr'], 521 cosmology['omega_M_0'], 522 cosmology['omega_b_0'], 523 cosmology['omega_n_0'], 524 cosmology['N_nu'], 525 cosmology['omega_lambda_0'], 526 cosmology['h'], 527 cosmology['baryonic_effects'],) 528 sigma = numpy.sqrt(sigmasq_0) * fg 529 530 # Propagate the error on sigmasq_0 to sigma. 531 error = fg * errorsq_0 / (2. * sigmasq_0) 532 533 return sigma, error

534

535 -def norm_power(**cosmology):

536 """Normalize the power spectrum to the specified sigma_8. 537 538 Returns the factor deltaSqr. 539 540 """ 541 cosmology['deltaSqr'] = 1.0 542 deltaSqr = (cosmology['sigma_8'] / 543 sigma_r(8.0 / cosmology['h'], 0.0, **cosmology)[0] 544 )**2.0 545 #print " deltaSqr = %.3g" % deltaSqr 546 547 del cosmology['deltaSqr'] 548 sig8 = sigma_r(8.0 / cosmology['h'], 0.0, deltaSqr=deltaSqr, 549 **cosmology)[0] 550 #print " Input sigma_8 = %.3g" % cosmology['sigma_8'] 551 #print " Numerical sigma_8 = %.3g" % sig8 552 sigma_8_error = (sig8 - cosmology['sigma_8'])/cosmology['sigma_8'] 553 if sigma_8_error > 1e-4: 554 warnings.warn("High sigma_8 fractional error = %.3g" % sigma_8_error) 555 return deltaSqr

556

557 #def rho_crit(h): 558 # return 3.0 * (h * H100)**2. / ( 8.0 * math.pi * G) 559 560 -def power_spectrum(k, z, **cosmology):

561 r"""The matter power spectrum P(k,z). 562 563 Uses equation 25 of Eisenstein & Hu (1999 ApJ 511 5). 564 565 Parameters 566 ---------- 567 568 k should be in Mpc^-1 569 570 Cosmological Parameters 571 ----------------------- 572 573 Uses 'n', and either 'sigma_8' or 'deltaSqr', as well as, for 574 transfer_function_EH, 'omega_M_0', 'omega_b_0', 'omega_n_0', 575 'N_nu', 'omega_lambda_0', and 'h'. 576 577 578 Notes 579 ----- 580 581 :: 582 583 P(k,z) = \delta^2 \frac{2 \pi^2}{k^3} \left(\frac{c k}{h 584 H_{100}}\right)^{3+n} \left(T(k,z) \frac{D_1(z)}{D_1(0)}\right)^2 585 586 Using the non-dependence of the transfer function on redshift, we can 587 rewrite this as 588 589 :: 590 591 P(k,z) = P(k,0) \left( \frac{D_1(z)}{D_1(0)} \right)^2 592 593 which is used by sigma_r to the z-dependence out of the integral. 594 595 """ 596 597 omega_M_0 = cosmology['omega_M_0'] 598 n = cosmology['n'] 599 h = cosmology['h'] 600 if 'deltaSqr' in cosmology: 601 deltaSqr = cosmology['deltaSqr'] 602 else: 603 deltaSqr = norm_power(**cosmology) 604 605 transFunc = transfer_function_EH(k, **cosmology)[0] 606 607 # This equals D1(z)/D1(0) 608 growthFact = fgrowth(z, omega_M_0) 609 610 # Expression: 611 # (sigma8/sig8)^2 * (2 pi^2) * k^-3. * (c * k / H_0)^(3 + n) * (tF * D)^2) 612 613 # Simplifies to: 614 # k^n (sigma8/sig8)^2 * (2 pi^2) * (c / H_0)^(3 + n) * (tF * D)^2 615 616 # compare w/ icosmo: 617 # k^n * tk^2 * (2 pi^2) * (sigma8/sig8)^2 * D^2 618 # Just a different normalization. 619 620 ps = (deltaSqr * (2. * math.pi**2.) * k**n * 621 (cc.c_light_Mpc_s / (h * cc.H100_s))**(3. + n) * 622 (transFunc * growthFact)**2.) 623 624 return ps

625

626 -def volume_radius_dmdr(mass, **cosmology):

627 """The volume, radius, and dm/dr for a sphere of the given mass. 628 629 Uses the mean density of the universe. 630 631 Parameters 632 ---------- 633 634 mass: array 635 mass of the sphere in Solar Masses, M_sun. 636 637 Returns 638 ------- 639 640 volume in Mpc^3 641 radius in Mpc 642 dmdr in Msun / Mpc 643 644 """ 645 rho_crit, rho_0 = cden.cosmo_densities(**cosmology) 646 647 volume = mass / rho_0 648 r = (volume / ((4. / 3.) * math.pi))**(1./3.) 649 650 dmdr = 4. * math.pi * r**2. * rho_0 651 652 return volume, r, dmdr

653

654 -def mass_to_radius(mass, **cosmology):

655 """The radius in Mpc of a sphere of the given mass. 656 657 Parameters 658 ----------- 659 660 mass in Msun 661 662 Returns 663 ------- 664 665 radius in Mpc 666 667 Notes 668 ----- 669 670 This is a convenience function that calls volume_radius_dmdr and 671 returns only the radius. 672 673 """ 674 volume, r, dmdr = volume_radius_dmdr(mass, **cosmology) 675 return r

676

677 -def radius_to_mass(r, **cosmology):

678 """The mass of a sphere of radius r in Mpc. 679 680 Uses the mean density of the universe. 681 682 """ 683 volume = (4./3.) * math.pi * r**3. 684 685 if 'rho_0' in cosmology: 686 rho_0 = cosmology['rho_0'] 687 else: 688 rho_crit, rho_0 = cden.cosmo_densities(**cosmology) 689 690 mass = volume * rho_0 691 return mass

692

693 694 -def virial_temp(mass, z, mu=None, **cosmology):

695 r"""The Virial temperature for a halo of a given mass. 696 697 Calculates the Virial temperature in Kelvin for a halo of a given 698 mass using equation 26 of Barkana & Loeb. 699 700 The transition from neutral to ionized is assumed to occur at temp 701 = 1e4K. At temp >= 10^4 k, the mean partical mass drops from 1.22 702 to 0.59 to very roughly account for collisional ionization. 703 704 Parameters 705 ---------- 706 707 mass: array 708 Mass in Solar Mass units. 709 710 z: array 711 Redshift. 712 713 mu: array, optional 714 Mean mass per particle. 715 716 """ 717 718 omega_M_0 = cosmology['omega_M_0'] 719 omega = cden.omega_M_z(z, **cosmology) 720 d = omega - 1 721 deltac = 18. * math.pi**2. + 82. * d - 39. * d**2. 722 723 if mu is None: 724 mu_t = 1.22 725 else: 726 mu_t = mu 727 temp = (1.98e4 * 728 (mass * cosmology['h']/1.0e8)**(2.0/3.0) * 729 (omega_M_0 * deltac / (omega * 18. * math.pi**2.))**(1.0/3.0) * 730 ((1 + z) / 10.0) * 731 (mu_t / 0.6)) 732 733 if mu is None: 734 # Below is some magic to consistently use mu = 1.22 at temp < 1e4K 735 # and mu = 0.59 at temp >= 1e4K. 736 t_crit = 1e4 737 t_crit_large = 1.22 * t_crit / 0.6 738 t_crit_small = t_crit 739 mu = ((temp < t_crit_small) * 1.22 + 740 (temp > t_crit_large) * 0.59 + 741 (1e4 * 1.22/temp) * 742 (temp >= t_crit_small) * (temp <= t_crit_large)) 743 temp = temp * mu / 1.22 744 745 return temp

746

747 -def virial_mass(temp, z, mu=None, **cosmology):

748 r"""The mass of a halo of the given Virial temperature. 749 750 Uses equation 26 of Barkana & Loeb (2001PhR...349..125B), solved 751 for T_vir as a function of mass. 752 753 Parameters 754 ---------- 755 756 temp: array 757 Virial temperature of the halo in Kelvin. 758 759 z: array 760 Redshift. 761 762 Returns 763 ------- 764 765 mass: array 766 The mass of such a halo in Solar Masses. 767 768 Notes 769 ----- 770 771 At temp >= 10^4 k, the mean partical mass drops from 1.22 to 0.59 772 to very roughly account for collisional ionization. 773 774 Examples 775 -------- 776 777 >>> cosmo = {'omega_M_0' : 0.27, 778 ... 'omega_lambda_0' : 1-0.27, 779 ... 'omega_b_0' : 0.045, 780 ... 'omega_n_0' : 0.0, 781 ... 'N_nu' : 0, 782 ... 'h' : 0.72, 783 ... 'n' : 1.0, 784 ... 'sigma_8' : 0.9 785 ... } 786 >>> mass = virial_mass(1e4, 6.0, **cosmo) 787 >>> temp = virial_temp(mass, 6.0, **cosmo) 788 >>> print "Mass = %.3g M_sun" % mass 789 Mass = 1.68e+08 M_sun 790 >>> print round(temp, 4) 791 10000.0 792 793 """ 794 if mu is None: 795 t_crit = 1e4 796 mu = (temp < t_crit) * 1.22 + (temp >= t_crit) * 0.59 797 798 divisor = virial_temp(1.0e8 / cosmology['h'], z, mu=mu, **cosmology) 799 return 1.0e8 * (temp/divisor)**(3.0/2.0) / cosmology['h']

800

801 -def virial_temp_HB(mass, z):

802 """Virial temperature from halo mass according to Haiman & Bryan 803 (2006ApJ...650....7). 804 805 z is the redshift. 806 807 Units are Msun and kelvin. 808 809 """ 810 return 1800. * (mass/1e6)**(2./3.) * (1.+z)/21

811

812 -def virial_mass_HB(temp, z):

813 """Halo mass from Virial temperature according to Haiman & Bryan 814 (2006ApJ...650....7). 815 816 z is the redshift. 817 818 Units are Msun and kelvin. 819 820 """ 821 return 1e6 * (21. * temp / (1800 * (1+z)))**(3./2.)

822

823 -def sig_del(temp_min, z, mass=None, passed_min_mass = False, **cosmology):

824 """Convenience function to calculate collapse fraction inputs. 825 826 Parameters 827 ---------- 828 829 temp_min: 830 Minimum Virial temperature for a halo to be counted. Or minimum 831 mass, if passed_min_mass is True. 832 833 z: 834 Redshift. 835 836 mass: optional 837 The mass of the region under consideration. Defaults to 838 considering the entire universe. 839 840 passed_min_mass: boolean 841 Indicates that the first argument is actually the minimum mass, 842 not the minimum Virial temperature. 843 844 """ 845 846 if passed_min_mass: 847 mass_min = temp_min 848 else: 849 mass_min = virial_mass(temp_min, z, **cosmology) 850 r_min = mass_to_radius(mass_min, **cosmology) 851 sigma_min = sigma_r(r_min, 0., **cosmology) 852 sigma_min = sigma_min[0] 853 854 fg = fgrowth(z, cosmology['omega_M_0']) 855 delta_c = 1.686 / fg 856 857 if mass is None: 858 return sigma_min, delta_c 859 else: 860 r_mass = mass_to_radius(mass) 861 sigma_mass = sigma_r(r_mass, 0., **cosmology) 862 return sigma_min, delta_c, sigma_mass

863

864 -def collapse_fraction(sigma_min, delta_crit, sigma_mass=0, delta=0):

865 r"""Fraction of mass contained in collapsed objects. 866 867 Use sig_del to conveniently obtain sigma_min and delta_crit. See 868 Examples velow. 869 870 Parameters 871 ---------- 872 873 sigma_min: 874 The standard deviatiation of density fluctuations on the scale 875 corresponding to the minimum mass for a halo to be counted. 876 877 delta_crit: 878 The critical (over)density of collapse. 879 880 sigma_mass: 881 The standard deviation of density fluctuations on the scale 882 corresponding to the mass of the region under 883 consideration. Use zero to consider the entire universe. 884 885 delta: 886 The overdensity of the region under consideration. Zero 887 corresponds to the mean density of the universe. 888 889 Notes 890 ----- 891 892 The fraction of the mass in a region of mass m that has already 893 collapsed into halos above mass m_min is: 894 895 :: 896 897 f_\mathrm{col} = \mathrm{erfc} \left[ \frac{\delta_c - \delta(m)} 898 { \sqrt {2 [\sigma^2(m_\mathrm{min}) - \sigma^2(m)]}} \right] 899 900 901 The answer isn't real if sigma_mass > sigma_min. 902 903 Note that there is a slight inconsistency in the mass used to 904 calculate sigma in the above formula, since the region deviates 905 from the average density. 906 907 Examples 908 -------- 909 910 >>> import numpy 911 >>> import perturbation as cp 912 >>> cosmo = {'omega_M_0' : 0.27, 913 ... 'omega_lambda_0' : 1-0.27, 914 ... 'omega_b_0' : 0.045, 915 ... 'omega_n_0' : 0.0, 916 ... 'N_nu' : 0, 917 ... 'h' : 0.72, 918 ... 'n' : 1.0, 919 ... 'sigma_8' : 0.9, 920 ... 'baryonic_effects' : False 921 ... } 922 >>> fc = cp.collapse_fraction(*cp.sig_del(1e4, 0, **cosmo)) 923 >>> print round(fc, 4) 924 0.7328 925 >>> fc = cp.collapse_fraction(*cp.sig_del(1e2, 0, **cosmo)) 926 >>> print round(fc, 4) 927 0.8571 928 929 """ 930 931 fraction = scipy.special.erfc((delta_crit - delta) / 932 (numpy.sqrt(2. * 933 (sigma_min**2. - sigma_mass**2.) 934 ) 935 )) 936 return fraction

937 938 if __name__ == "__main__": 939 import doctest 940 doctest.testmod() 941

Source Code for Module cosmolopy.perturbation