diff options
Diffstat (limited to 'fg21sim/extragalactic')
| -rw-r--r-- | fg21sim/extragalactic/clusters/cosmology.py | 118 |
1 files changed, 25 insertions, 93 deletions
diff --git a/fg21sim/extragalactic/clusters/cosmology.py b/fg21sim/extragalactic/clusters/cosmology.py index bedb17e..5f7b283 100644 --- a/fg21sim/extragalactic/clusters/cosmology.py +++ b/fg21sim/extragalactic/clusters/cosmology.py @@ -2,7 +2,7 @@ # MIT license """ -Cosmological models +Flat ΛCDM cosmological model. """ import logging @@ -18,30 +18,30 @@ logger = logging.getLogger(__name__) class Cosmology: """ - Cosmological model. + Flat ΛCDM cosmological model. Attributes ---------- H0 : float Hubble parameter at present day (z=0) Om0 : float - Density parameter of matter at present day + 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 - 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 + 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, Ode0=None, sigma8=None): + def __init__(self, H0=71.0, Om0=0.27, Ob0=0.046, Ode0=None, sigma8=0.834): 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!") @@ -49,22 +49,8 @@ class Cosmology: self.Om0 = Om0 self.Ob0 = Ob0 self.Ode0 = Ode0 - self._sigma8 = sigma8 + self.sigma8 = sigma8 self.cosmo = LambdaCDM(H0=H0, Om0=Om0, Ob0=Ob0, 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): @@ -74,32 +60,6 @@ class Cosmology: 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. @@ -176,17 +136,8 @@ class Cosmology: 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)) + omega_z = (1 / self.OmegaM(z)) - 1 + Delta_c = 18*np.pi**2 * (1 + 0.4093 * omega_z**0.9052) return Delta_c def overdensity_crit(self, z): @@ -200,21 +151,11 @@ class Cosmology: 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)) + 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)) return delta_c def growth_factor(self, z): @@ -227,19 +168,10 @@ class Cosmology: 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)) + 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)[0] + D = coef * integral return D |