cosmology.py: Only support flat ΛCDM cosmological model

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