Add clusters/cosmology.py with common cosmological models

Mainly the "flat LambdaCDM" and "EdS" cosmological models.

With some useful utility functions.

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