diff options
Diffstat (limited to 'fg21sim/utils')
| -rw-r--r-- | fg21sim/utils/cosmology.py | 215 | ||||
| -rw-r--r-- | fg21sim/utils/units.py | 72 |
2 files changed, 287 insertions, 0 deletions
diff --git a/fg21sim/utils/cosmology.py b/fg21sim/utils/cosmology.py new file mode 100644 index 0000000..12e581a --- /dev/null +++ b/fg21sim/utils/cosmology.py @@ -0,0 +1,215 @@ +# Copyright (c) 2016-2017 Weitian LI <liweitianux@live.com> +# MIT license + +""" +Flat ΛCDM cosmological model. +""" + +import numpy as np +from scipy import integrate +from astropy.cosmology import FlatLambdaCDM + +from .units import (UnitConversions as AUC, Constants as AC) + + +class Cosmology: + """ + Flat ΛCDM cosmological model. + + Attributes + ---------- + H0 : float + Hubble parameter at present day (z=0) + Om0 : float + Density parameter of (dark and baryon) matter at present day + Ob0 : float + Density parameter of baryon at present day + Ode0 : float + Density parameter of dark energy at present day + sigma8 : float + Present-day rms density fluctuation on a scale of 8 h^-1 Mpc. + + References + ---------- + [1] https://astro.uni-bonn.de/~pavel/WIKIPEDIA/Lambda-CDM_model.html + [2] https://en.wikipedia.org/wiki/Lambda-CDM_model + [3] Randall, Sarazin & Ricker 2002, ApJ, 577, 579 + http://adsabs.harvard.edu/abs/2002ApJ...577..579R + Sec.(2) + """ + def __init__(self, H0=71.0, Om0=0.27, Ob0=0.046, sigma8=0.834): + self.H0 = H0 # [km/s/Mpc] + self.Om0 = Om0 + self.Ob0 = Ob0 + self.Ode0 = 1.0 - Om0 + self.sigma8 = sigma8 + self._cosmo = FlatLambdaCDM(H0=H0, Om0=Om0, Ob0=Ob0) + + @property + def h(self): + """ + Dimensionless/reduced Hubble parameter + """ + return self.H0 / 100.0 + + @property + def M8(self): + """ + Mass contained in a sphere of radius of 8 h^-1 Mpc. + Unit: [Msun] + """ + r = 8 * AUC.Mpc2cm / self.h # [cm] + M8 = (4*np.pi/3) * r**3 * self.rho_crit(0) # [g] + M8 *= AUC.g2Msun # [Msun] + return M8 + + def E(self, z): + """ + Redshift evolution factor. + """ + return np.sqrt(self.Om0 * (1+z)**3 + self.Ode0) + + def H(self, z): + """ + Hubble parameter at redshift z. + """ + return self.H0 * self.E(z) + + @property + def hubble_time(self): + """ + Hubble time. + Unit: [Gyr] + """ + uconv = AUC.Mpc2km * AUC.s2Gyr + t_H = (1.0/self.H0) * uconv # [Gyr] + return t_H + + def age(self, z): + """ + Cosmic time (age) at redshift z. + + Parameters + ---------- + z : float + Redshift + + Returns + ------- + age : float + Age of the universe (cosmic time) at the given redshift. + Unit: [Gyr] + + References + ---------- + [1] Thomas & Kantowski 2000, Physical Review D, 62, 103507 + http://adsabs.harvard.edu/abs/2000PhRvD..62j3507T + Eq.(18) + """ + t_H = self.hubble_time + t = (t_H * (2/3/np.sqrt(1-self.Om0)) * + np.arcsinh(np.sqrt((1/self.Om0 - 1) / (1+z)**3))) + return t + + @property + def age0(self): + """ + Present age of the universe. + """ + return self.age(0) + + def redshift(self, age): + """ + Invert the above ``self.age(z)`` formula analytically, to calculate + the redshift corresponding to the given cosmic time (age). + + Parameters + ---------- + age : float + Age of the universe (cosmic time), unit [Gyr] + + Returns + ------- + z : float + Redshift corresponding to the specified age. + """ + t_H = self.hubble_time + term1 = (1/self.Om0) - 1 + term2 = np.sinh(3*age*np.sqrt(1-self.Om0) / (2*t_H)) ** 2 + z = (term1 / term2) ** (1/3) - 1 + return z + + def rho_crit(self, z): + """ + Critical density at redshift z. + Unit: [g/cm^3] + """ + rho = 3 * self.H(z)**2 / (8*np.pi * AC.G) + rho *= AUC.km2Mpc**2 + return rho + + def Om(self, z): + """ + Density parameter of matter at redshift z. + """ + return self.Om0 * (1+z)**3 / self.E(z)**2 + + def overdensity_virial(self, z): + """ + Calculate the virial overdensity, which generally used to + determine the virial radius of a cluster. + + References + ---------- + [1] Cassano & Brunetti 2005, MNRAS, 357, 1313 + http://adsabs.harvard.edu/abs/2005MNRAS.357.1313C + Eqs.(10,A4) + """ + omega_z = (1 / self.Om(z)) - 1 + Delta_c = 18*np.pi**2 * (1 + 0.4093 * omega_z**0.9052) + return Delta_c + + def overdensity_crit(self, z): + """ + Critical (linear) overdensity for a region to collapse + at a redshift z. + + References + ---------- + [1] Randall, Sarazin & Ricker 2002, ApJ, 577, 579 + http://adsabs.harvard.edu/abs/2002ApJ...577..579R + Appendix.A, Eq.(A1) + """ + coef = 3 * (12*np.pi) ** (2/3) / 20 + D0 = self.growth_factor0 + D_z = self.growth_factor(z) + Om_z = self.Om(z) + delta_c = coef * (D0 / D_z) * (1 + 0.0123*np.log10(Om_z)) + return delta_c + + def growth_factor(self, z): + """ + Growth factor at redshift z. + + References + ---------- + [1] Randall, Sarazin & Ricker 2002, ApJ, 577, 579 + http://adsabs.harvard.edu/abs/2002ApJ...577..579R + Appendix.A, Eq.(A7) + """ + x0 = (2 * self.Ode0 / self.Om0) ** (1/3) + x = x0 / (1 + z) + coef = np.sqrt(x**3 + 2) / (x**1.5) + integral = integrate.quad(lambda y: y**1.5 * (y**3+2)**(-1.5), + a=0, b=x, epsabs=1e-5, epsrel=1e-5)[0] + D = coef * integral + return D + + @property + def growth_factor0(self): + """ + Present-day (z=0) growth factor. + """ + if not hasattr(self, "_growth_factor0"): + self._growth_factor0 = self.growth_factor(0) + return self._growth_factor0 diff --git a/fg21sim/utils/units.py b/fg21sim/utils/units.py new file mode 100644 index 0000000..a7e4ef5 --- /dev/null +++ b/fg21sim/utils/units.py @@ -0,0 +1,72 @@ +# Copyright (c) 2017 Weitian LI <liweitianux@live.com> +# MIT license + +""" +Commonly used units and their conversions relations, as well as constants. + +Astropy's units system is very powerful, but also very slow, +and may even be the speed bottleneck of the program. + +This module provides commonly used units conversions by holding +them directly in a class, thus avoid repeated/unnecessary calculations. +""" + +import astropy.units as au +import astropy.constants as ac + + +class Units: + """ + Commonly used units, especially in the CGS unit system. + """ + # Unit for electron momentum (p), thus its value is the Lorentz factor + mec = ac.m_e.cgs.value*ac.c.cgs.value # [g cm / s] + + +class UnitConversions: + """ + Commonly used units conversion relations. + + Hold the conversion relations directly to avoid repeated/unnecessary + calculations. + """ + # Mass + Msun2g = au.solMass.to(au.g) + g2Msun = au.g.to(au.solMass) + # Time + Gyr2s = au.Gyr.to(au.s) + s2Gyr = au.s.to(au.Gyr) + # Length + kpc2cm = au.kpc.to(au.cm) + cm2kpc = au.cm.to(au.kpc) + Mpc2cm = au.Mpc.to(au.cm) + cm2Mpc = au.cm.to(au.Mpc) + Mpc2km = au.Mpc.to(au.km) + km2Mpc = au.km.to(au.Mpc) + kpc2km = au.kpc.to(au.km) + km2kpc = au.km.to(au.kpc) + km2cm = au.km.to(au.cm) + # Energy + keV2erg = au.keV.to(au.erg) + + +class Constants: + """ + Commonly used constants, especially in the CGS unit system. + + Astropy's constants are stored in SI units by default. + When request a constant in CGS unit system, additional (and slow) + conversions required. + """ + # Speed of light + c = ac.c.cgs.value # [cm/s] + # Atomic mass unit (i.e., a.m.u.) + u = ac.u.cgs.value # [g] + # Gravitational constant + G = ac.G.cgs.value # [cm^3/g/s^2] + # Electron charge + e = ac.e.gauss.value # [Fr] = [esu] + + # Mean molecular weight + # Ref.: Ettori et al, 2013, Space Science Review, 177, 119-154, Eq.(6) + mu = 0.6 |