From c630afcff2a8529232dee68a3964e52cf8f080b4 Mon Sep 17 00:00:00 2001 From: Aaron LI Date: Sat, 27 May 2017 18:11:00 +0800 Subject: clusters: Move units.py and cosmology.py to utils --- fg21sim/extragalactic/clusters/cosmology.py | 215 ---------------------------- fg21sim/extragalactic/clusters/emission.py | 2 +- fg21sim/extragalactic/clusters/formation.py | 2 +- fg21sim/extragalactic/clusters/halo.py | 6 +- fg21sim/extragalactic/clusters/units.py | 72 ---------- fg21sim/utils/cosmology.py | 215 ++++++++++++++++++++++++++++ fg21sim/utils/units.py | 72 ++++++++++ 7 files changed, 293 insertions(+), 291 deletions(-) delete mode 100644 fg21sim/extragalactic/clusters/cosmology.py delete mode 100644 fg21sim/extragalactic/clusters/units.py create mode 100644 fg21sim/utils/cosmology.py create mode 100644 fg21sim/utils/units.py 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 -# 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 -# 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 diff --git a/fg21sim/utils/cosmology.py b/fg21sim/utils/cosmology.py new file mode 100644 index 0000000..12e581a --- /dev/null +++ b/fg21sim/utils/cosmology.py @@ -0,0 +1,215 @@ +# Copyright (c) 2016-2017 Weitian LI +# 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/utils/units.py b/fg21sim/utils/units.py new file mode 100644 index 0000000..a7e4ef5 --- /dev/null +++ b/fg21sim/utils/units.py @@ -0,0 +1,72 @@ +# Copyright (c) 2017 Weitian LI +# 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 -- cgit v1.2.2