diff options
Diffstat (limited to 'fg21sim')
| -rw-r--r-- | fg21sim/extragalactic/clusters/cosmology.py | 230 |
1 files changed, 230 insertions, 0 deletions
diff --git a/fg21sim/extragalactic/clusters/cosmology.py b/fg21sim/extragalactic/clusters/cosmology.py new file mode 100644 index 0000000..a66ff66 --- /dev/null +++ b/fg21sim/extragalactic/clusters/cosmology.py @@ -0,0 +1,230 @@ +# Copyright (c) 2016-2017 Weitian LI <liweitianux@live.com> +# MIT license + +""" +Cosmological models +""" + +import logging + +import numpy as np +from scipy import integrate +import astropy.units as au +from astropy.cosmology import LambdaCDM, z_at_value + + +logger = logging.getLogger(__name__) + + +class Cosmo: + """ + Cosmological model. + + Attributes + ---------- + H0 : float + Hubble parameter at present day (z=0) + Om0 : float + Density parameter of matter at present day + Ode0 : float + Density parameter of dark energy at present day + model : str + Type of the current cosmological model: + * open : Om0 < 1, Ode0 = 0 + * closed : Om0 > 1, Ode0 = 0 + * EdS (Einstein-de Sitter) : Om0 = 1, Ode0 = 0 + * flatLambdaCDM : Om0 + Ode0 = 1, Ode0 > 0 + """ + + def __init__(self, H0=71.0, Om0=0.27, Ode0=None, sigma8=None): + Ode0 = 1.0 - Om0 if Ode0 is None else Ode0 + if (Ode0 > 0) and abs(Om0 + Ode0 - 1) > 1e-5: + raise ValueError("non-flat LambdaCDM model not supported!") + self.H0 = H0 # [km/s/Mpc] + self.Om0 = Om0 + self.Ode0 = Ode0 + self._sigma8 = sigma8 + self.cosmo = LambdaCDM(H0=H0, Om0=Om0, Ode0=Ode0) + logger.info("Cosmological model: {0}".format(self.model)) + + @property + def model(self): + if self.Ode0 < 1e-5: + if self.Om0 < 1: + model = "open" + elif self.Om0 == 1: + model = "EdS" + else: + model = "closed" + else: + model = "flatLambdaCDM" + return model + + @property + def h(self): + """ + Dimensionless/reduced Hubble parameter + """ + return self.H0 / 100.0 + + @property + def sigma8(self): + """ + Present-day rms density fluctuation on a scale of 8 h^-1 Mpc. + + References + ---------- + [1] Randall, Sarazin & Ricker 2002, ApJ, 577, 579 + http://adsabs.harvard.edu/abs/2002ApJ...577..579R + Sec.2 + """ + if hasattr(self, "_sigma8"): + return self._sigma8 + # + if self.model == "open": + sigma8 = 0.827 + elif self.model == "closed": + raise NotImplementedError + elif self.model == "EdS": + sigma8 = 0.514 + elif self.model == "flatLambdaCDM": + sigma8 = 0.834 + else: + raise ValueError("unknown model: {0}".format(self.model)) + return sigma8 + + @property + def M8(self): + """ + Mass contained in a sphere of radius of 8 h^-1 Mpc. + Unit: [Msun] + """ + r = 8 * au.Mpc.to(au.cm) / self.h # [cm] + M8 = (4*np.pi/3) * r**3 * self.rho_crit(0) + return (M8 * au.g.to(au.solMass)) + + def age(self, z): + """ + Cosmic time at redshift z. + + Parameters + ---------- + z : float + Redshift + + Returns + ------- + age : float + Age of the universe (cosmic time) at the given redshift. + Unit: [Gyr] + """ + return self.cosmo.age(z).value + + def redshift(self, age): + """ + Invert the above age calculation, to return the redshift + corresponding to the given cosmic time. + + Parameters + ---------- + age : float + Age of the universe (cosmic time), unit [Gyr] + + Returns + ------- + z : float + Redshift corresponding to the input age. + """ + return z_at_value(self.age, age) + + def rho_crit(self, z): + """ + Critical density at redshift z. + Unit: [g/cm^3] + """ + return self.cosmo.critical_density(0).value + + def OmegaM(self, z): + """ + Density parameter of matter at redshift z. + """ + return self.Om0 * (1+z)**3 / self.cosmo.efunc(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) + """ + if self.model == "open": + raise NotImplementedError + elif self.model == "closed": + raise NotImplementedError + elif self.model == "EdS": + Delta_c = 18 * np.pi**2 + elif self.model == "flatLambdaCDM": + omega_z = (1 / self.OmegaM(z)) - 1 + Delta_c = 18*np.pi**2 * (1 + 0.4093 * omega_z**0.9052) + else: + raise ValueError("unknown model: {0}".format(self.model)) + 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) + """ + if self.model == "open": + raise NotImplementedError + elif self.model == "closed": + raise NotImplementedError + elif self.model == "EdS": + coef = 3 * (12*np.pi) ** (2/3) / 20 + delta_c = coef * (self.age(0) / self.age(z)) ** (2/3) + elif self.model == "flatLambdaCDM": + coef = 3 * (12*np.pi) ** (2/3) / 20 + D0 = self.growth_factor(0) + D_z = self.growth_factor(z) + Om_z = self.OmegaM(z) + delta_c = coef * (D0 / D_z) * (1 + 0.0123*np.log10(Om_z)) + else: + raise ValueError("unknown model: {0}".format(self.model)) + 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) + """ + if self.model == "open": + raise NotImplementedError + elif self.model == "closed": + raise NotImplementedError + elif self.model == "EdS": + raise NotImplementedError + elif self.model == "flatLambdaCDM": + 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), + 0, x)[0] + D = coef * integral + else: + raise ValueError("unknown model: {0}".format(self.model)) + return D |