aboutsummaryrefslogtreecommitdiffstats
path: root/fg21sim/extragalactic/clusters
diff options
context:
space:
mode:
Diffstat (limited to 'fg21sim/extragalactic/clusters')
-rw-r--r--fg21sim/extragalactic/clusters/cosmology.py215
-rw-r--r--fg21sim/extragalactic/clusters/emission.py2
-rw-r--r--fg21sim/extragalactic/clusters/formation.py2
-rw-r--r--fg21sim/extragalactic/clusters/halo.py6
-rw-r--r--fg21sim/extragalactic/clusters/units.py72
5 files changed, 6 insertions, 291 deletions
diff --git a/fg21sim/extragalactic/clusters/cosmology.py b/fg21sim/extragalactic/clusters/cosmology.py
deleted file mode 100644
index 12e581a..0000000
--- a/fg21sim/extragalactic/clusters/cosmology.py
+++ /dev/null
@@ -1,215 +0,0 @@
-# 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/extragalactic/clusters/emission.py b/fg21sim/extragalactic/clusters/emission.py
index e71632b..1374b52 100644
--- a/fg21sim/extragalactic/clusters/emission.py
+++ b/fg21sim/extragalactic/clusters/emission.py
@@ -11,7 +11,7 @@ import numpy as np
import scipy.integrate
import scipy.special
-from .units import (Units as AU, Constants as AC)
+from ...utils.units import (Units as AU, Constants as AC)
logger = logging.getLogger(__name__)
diff --git a/fg21sim/extragalactic/clusters/formation.py b/fg21sim/extragalactic/clusters/formation.py
index f8c6d3a..1f48433 100644
--- a/fg21sim/extragalactic/clusters/formation.py
+++ b/fg21sim/extragalactic/clusters/formation.py
@@ -20,8 +20,8 @@ import scipy.integrate
import scipy.special
import scipy.optimize
-from .cosmology import Cosmology
from .mergertree import MergerTree
+from ...utils.cosmology import Cosmology
logger = logging.getLogger(__name__)
diff --git a/fg21sim/extragalactic/clusters/halo.py b/fg21sim/extragalactic/clusters/halo.py
index a6d5a20..02bd74b 100644
--- a/fg21sim/extragalactic/clusters/halo.py
+++ b/fg21sim/extragalactic/clusters/halo.py
@@ -18,10 +18,12 @@ import scipy.interpolate
import scipy.integrate
import scipy.optimize
-from .cosmology import Cosmology
from .formation import ClusterFormation
from .solver import FokkerPlanckSolver
-from .units import (Units as AU, UnitConversions as AUC, Constants as AC)
+from ...utils.cosmology import Cosmology
+from ...utils.units import (Units as AU,
+ UnitConversions as AUC,
+ Constants as AC)
logger = logging.getLogger(__name__)
diff --git a/fg21sim/extragalactic/clusters/units.py b/fg21sim/extragalactic/clusters/units.py
deleted file mode 100644
index a7e4ef5..0000000
--- a/fg21sim/extragalactic/clusters/units.py
+++ /dev/null
@@ -1,72 +0,0 @@
-# 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
Copyright © 2017-2018 Aaron LI. All Rights Reserved.