aboutsummaryrefslogtreecommitdiffstats
path: root/fg21sim
diff options
context:
space:
mode:
Diffstat (limited to 'fg21sim')
-rw-r--r--fg21sim/extragalactic/clusters/emission.py7
-rw-r--r--fg21sim/extragalactic/clusters/formation.py19
-rw-r--r--fg21sim/extragalactic/clusters/halo.py46
-rw-r--r--fg21sim/extragalactic/clusters/main.py15
4 files changed, 34 insertions, 53 deletions
diff --git a/fg21sim/extragalactic/clusters/emission.py b/fg21sim/extragalactic/clusters/emission.py
index 93e04ff..4366025 100644
--- a/fg21sim/extragalactic/clusters/emission.py
+++ b/fg21sim/extragalactic/clusters/emission.py
@@ -18,11 +18,11 @@ import numpy as np
import scipy.integrate
import scipy.special
+from ...utils import cosmo
from ...utils.units import (Units as AU,
UnitConversions as AUC,
Constants as AC)
from ...utils.convert import Fnu_to_Tb_fast
-from ...utils.cosmology import Cosmology
logger = logging.getLogger(__name__)
@@ -57,7 +57,6 @@ class SynchrotronEmission:
self.n_e = n_e
self.z = z
self.radius = radius # [kpc]
- self.cosmo = Cosmology()
@property
def frequency_larmor(self):
@@ -165,7 +164,7 @@ class SynchrotronEmission:
Synchrotron flux at frequency ``nu``.
Unit: [Jy] = 1e-23 [erg/s/cm^2/Hz]
"""
- DL = self.cosmo.DL(self.z) * AUC.Mpc2cm # [cm]
+ DL = cosmo.DL(self.z) * AUC.Mpc2cm # [cm]
P_nu = self.power(nu)
F_nu = 1e23 * P_nu / (4*np.pi * DL*DL) # [Jy]
return F_nu
@@ -192,7 +191,7 @@ class SynchrotronEmission:
Synchrotron surface brightness at frequency ``nu``.
Unit: [K] <-> [Jy/pixel]
"""
- DA = self.cosmo.DL(self.z) * AUC.Mpc2cm # [cm]
+ DA = cosmo.DL(self.z) * AUC.Mpc2cm # [cm]
radius = self.radius * AUC.kpc2cm # [cm]
omega = (np.pi * radius**2 / DA**2) * AUC.rad2deg**2 # [deg^2]
pixelarea = (pixelsize * AUC.arcsec2deg) ** 2 # [deg^2]
diff --git a/fg21sim/extragalactic/clusters/formation.py b/fg21sim/extragalactic/clusters/formation.py
index 827c089..3bb44a9 100644
--- a/fg21sim/extragalactic/clusters/formation.py
+++ b/fg21sim/extragalactic/clusters/formation.py
@@ -21,7 +21,7 @@ import scipy.special
import scipy.optimize
from .mergertree import MergerTree
-from ...utils.cosmology import Cosmology
+from ...utils import cosmo
logger = logging.getLogger(__name__)
@@ -54,8 +54,6 @@ class ClusterFormation:
the merger is a major event or a minor one.
If ``M_main/M_sub < ratio_major``, then it is a major merger event.
Default: 3.0
- cosmo : `~Cosmology`, optional
- Adopted cosmological model with custom utility functions.
merger_mass_min : float, optional
Minimum mass change to be regarded as a merger event instead of
accretion.
@@ -70,12 +68,11 @@ class ClusterFormation:
major merger event, or None if not found.
"""
def __init__(self, M0, z0, zmax=3.0, ratio_major=3.0,
- cosmo=Cosmology(), merger_mass_min=1e12):
+ merger_mass_min=1e12):
self.M0 = M0 # [Msun]
self.z0 = z0
self.zmax = zmax
self.ratio_major = ratio_major
- self.cosmo = cosmo
self.merger_mass_min = merger_mass_min
@property
@@ -101,7 +98,7 @@ class ClusterFormation:
References: Ref.[1],Eq.(2)
"""
alpha = self.sigma_index
- sigma = self.cosmo.sigma8 * (mass / self.cosmo.M8) ** (-alpha)
+ sigma = cosmo.sigma8 * (mass / cosmo.M8) ** (-alpha)
return sigma
def f_delta_c(self, z):
@@ -113,7 +110,7 @@ class ClusterFormation:
References: Ref.[1],App.A,Eq.(A1)
"""
- return self.cosmo.overdensity_crit(z)
+ return cosmo.overdensity_crit(z)
def f_dw_max(self, mass):
"""
@@ -148,7 +145,7 @@ class ClusterFormation:
References: Ref.[1],Sec.(3)
"""
alpha = self.sigma_index
- mass = self.cosmo.M8 * (S / self.cosmo.sigma8**2)**(-1/(2*alpha))
+ mass = cosmo.M8 * (S / cosmo.sigma8**2)**(-1/(2*alpha))
return mass
@staticmethod
@@ -254,7 +251,7 @@ class ClusterFormation:
Mc = self.M0
mtree_root = MergerTree(data={"mass": Mc,
"z": zc,
- "age": self.cosmo.age(zc)})
+ "age": cosmo.age(zc)})
logger.debug("[main] z=%.4f : mass=%g [Msun]" % (zc, Mc))
mtree = mtree_root
@@ -275,7 +272,7 @@ class ClusterFormation:
dS = self.gen_dS(dw)
# Progenitor properties
z1 = self.calc_z(w2 + dw)
- age1 = self.cosmo.age(z1)
+ age1 = cosmo.age(z1)
S1 = S2 + dS
M1 = self.calc_mass(S1)
dM = Mc - M1
@@ -313,7 +310,7 @@ class ClusterFormation:
merger tree.
"""
z = 0.0 if _z is None else _z
- node_data = {"mass": M, "z": z, "age": self.cosmo.age(z)}
+ node_data = {"mass": M, "z": z, "age": cosmo.age(z)}
# Whether to stop the trace
if self.zmax is not None and z > self.zmax:
diff --git a/fg21sim/extragalactic/clusters/halo.py b/fg21sim/extragalactic/clusters/halo.py
index d7f351f..c78c2c3 100644
--- a/fg21sim/extragalactic/clusters/halo.py
+++ b/fg21sim/extragalactic/clusters/halo.py
@@ -26,7 +26,7 @@ import scipy.optimize
from .formation import ClusterFormation
from .solver import FokkerPlanckSolver
-from ...utils.cosmology import Cosmology
+from ...utils import cosmo
from ...utils.units import (Units as AU,
UnitConversions as AUC,
Constants as AC)
@@ -35,7 +35,7 @@ from ...utils.units import (Units as AU,
logger = logging.getLogger(__name__)
-class HaloSingle:
+class RadioHalo:
"""
Simulate a single (giant) radio halos following the "statistical
magneto-turbulent model" proposed by Cassano & Brunetti (2005).
@@ -63,8 +63,6 @@ class HaloSingle:
mec : float
Unit for electron momentum (p): mec = m_e * c, p = gamma * mec,
therefore value of p is the Lorentz factor.
- cosmo : `~Cosmology`
- Adopted cosmological model with custom utility functions.
mtree : `~MergerTree`
Merging history of this cluster.
"""
@@ -99,10 +97,6 @@ class HaloSingle:
self.buffer_np = self.configs.getn(comp+"/buffer_np")
self.time_step = self.configs.getn(comp+"/time_step")
self.injection_index = self.configs.getn(comp+"/injection_index")
- # Cosmology model
- self.H0 = self.configs.getn("cosmology/H0")
- self.OmegaM0 = self.configs.getn("cosmology/OmegaM0")
- self.cosmo = Cosmology(H0=self.H0, Om0=self.OmegaM0)
logger.info("Loaded and set up configurations")
@property
@@ -185,13 +179,13 @@ class HaloSingle:
Unit: [cm^-3 mec^-1]
"""
if zbegin is None:
- tstart = self.cosmo.age(self.zmax)
+ tstart = cosmo.age(self.zmax)
else:
- tstart = self.cosmo.age(zbegin)
+ tstart = cosmo.age(zbegin)
if zend is None:
- tstop = self.cosmo.age(self.z0)
+ tstop = cosmo.age(self.z0)
else:
- tstop = self.cosmo.age(zend)
+ tstop = cosmo.age(zend)
fpsolver = FokkerPlanckSolver(
xmin=self.pmin, xmax=self.pmax,
@@ -252,8 +246,8 @@ class HaloSingle:
Rvir : float
Virial radius (unit: kpc) of the cluster at given redshift
"""
- Dc = self.cosmo.overdensity_virial(z)
- rho = self.cosmo.rho_crit(z) # [g/cm^3]
+ Dc = cosmo.overdensity_virial(z)
+ rho = cosmo.rho_crit(z) # [g/cm^3]
R_vir = (3*mass*AUC.Msun2g / (4*np.pi * Dc * rho))**(1/3) # [cm]
R_vir *= AUC.cm2kpc # [kpc]
return R_vir
@@ -316,7 +310,7 @@ class HaloSingle:
http://adsabs.harvard.edu/abs/2005MNRAS.357.1313C
Eq.(12)
"""
- f_baryon = self.cosmo.Ob0 / self.cosmo.Om0
+ f_baryon = cosmo.Ob0 / cosmo.Om0
Rv = self._radius_virial(mass, z) * AUC.kpc2cm # [cm]
V = (4*np.pi / 3) * Rv**3 # [cm^3]
rho = f_baryon * mass*AUC.Msun2g / V # [g/cm^3]
@@ -346,7 +340,7 @@ class HaloSingle:
http://adsabs.harvard.edu/abs/2005MNRAS.357.1313C
Eq.(13)
"""
- f_baryon = self.cosmo.Ob0 / self.cosmo.Om0
+ f_baryon = cosmo.Ob0 / cosmo.Om0
M_ICM = mass * f_baryon * AUC.Msun2g # [g]
r *= AUC.kpc2cm # [cm]
Rv = self._radius_virial(mass, z) * AUC.kpc2cm # [cm]
@@ -477,12 +471,12 @@ class HaloSingle:
time : float
Elapsing time (unit: Gyr)
"""
- t_begin = self.cosmo.age(z_begin) # [Gyr]
+ t_begin = cosmo.age(z_begin) # [Gyr]
t_end = t_begin + time
- if t_end >= self.cosmo.age(0):
+ if t_end >= cosmo.age(0):
z_end = 0.0
else:
- z_end = self.cosmo.redshift(t_end)
+ z_end = cosmo.redshift(t_end)
return z_end
@property
@@ -542,7 +536,7 @@ class HaloSingle:
"""
zgrid, chigrid = self._chi_data
# XXX: force a minimal value instead of zero or too small
- chi_min = 1 / (10 * self.cosmo.age0)
+ chi_min = 1 / (10 * cosmo.age0)
try:
zi = np.where(z < zgrid)[0][0]
@@ -677,7 +671,7 @@ class HaloSingle:
"""
if not hasattr(self, "_electron_injection_rate"):
e_th = self.e_thermal # [erg/cm^3]
- age = self.cosmo.age(self.z0)
+ age = cosmo.age(self.z0)
term1 = (self.injection_index-2) * self.eta_e
term2 = e_th / (self.pmin * self.mec * AC.c) # [cm^-3]
term3 = 1.0 / (age * self.pmin) # [Gyr^-1 mec^-1]
@@ -714,7 +708,7 @@ class HaloSingle:
http://adsabs.harvard.edu/abs/2013AN....334..515D
Eq.(15)
"""
- z = self.cosmo.redshift(t)
+ z = cosmo.redshift(t)
chi = self._coef_acceleration(z) # [Gyr^-1]
# NOTE: Cassano & Brunetti's formula misses a factor of 2.
Dpp = chi * p**2 / 4 # [mec^2/Gyr]
@@ -760,7 +754,7 @@ class HaloSingle:
http://adsabs.harvard.edu/abs/2005MNRAS.357.1313C
Eq.(38)
"""
- z = self.cosmo.redshift(t)
+ z = cosmo.redshift(t)
n_th = self._n_thermal(self.M0, z)
coef = -3.3e-29 * AUC.Gyr2s / self.mec # [mec/Gyr]
dpdt = coef * n_th * (1 + np.log(p/n_th) / 75)
@@ -776,7 +770,7 @@ class HaloSingle:
http://adsabs.harvard.edu/abs/2005MNRAS.357.1313C
Eq.(39)
"""
- z = self.cosmo.redshift(t)
+ z = cosmo.redshift(t)
coef = -4.8e-4 * AUC.Gyr2s / self.mec # [mec/Gyr]
dpdt = (coef * (p*self.mec)**2 *
((self.magnetic_field/3.2)**2 + (1+z)**4))
@@ -793,7 +787,7 @@ class HaloSingle:
Energy density of the ICM (unit: erg/cm^3)
"""
mass = self.M0
- f_baryon = self.cosmo.Ob0 / self.cosmo.Om0
+ f_baryon = cosmo.Ob0 / cosmo.Om0
kT = self.kT_mass(mass) # [keV]
N = mass * AUC.Msun2g * f_baryon / (AC.mu * AC.u)
E_th = kT*AUC.keV2erg * N # [erg]
@@ -820,7 +814,7 @@ class HaloSingle:
Number density of the ICM
Unit: [cm^-3]
"""
- f_baryon = self.cosmo.Ob0 / self.cosmo.Om0
+ f_baryon = cosmo.Ob0 / cosmo.Om0
N = mass * AUC.Msun2g * f_baryon / (AC.mu * AC.u)
Rv = self._radius_virial(mass, z) * AUC.kpc2cm # [cm]
V = (4*np.pi / 3) * Rv**3 # [cm^3]
diff --git a/fg21sim/extragalactic/clusters/main.py b/fg21sim/extragalactic/clusters/main.py
index 438b860..b7c2d03 100644
--- a/fg21sim/extragalactic/clusters/main.py
+++ b/fg21sim/extragalactic/clusters/main.py
@@ -19,7 +19,7 @@ import pandas as pd
from .formation import ClusterFormation
from .halo import RadioHalo
from ...sky import get_sky
-from ...utils.cosmology import Cosmology
+from ...utils import cosmo
from ...utils.io import dataframe_to_csv
@@ -63,14 +63,6 @@ class GalaxyClusters:
self.checksum = self.configs.getn("output/checksum")
self.clobber = self.configs.getn("output/clobber")
- # Cosmology model
- self.H0 = self.configs.getn("cosmology/H0")
- self.OmegaM0 = self.configs.getn("cosmology/OmegaM0")
- self.Omegab0 = self.configs.getn("cosmology/Omegab0")
- self.sigma8 = self.configs.getn("cosmology/sigma8")
- self.cosmo = Cosmology(H0=self.H0, Om0=self.OmegaM0,
- Ob0=self.Omegab0, sigma8=self.sigma8)
-
# Sky and resolution
if self.sky.type_ == "patch":
self.resolution = self.sky.pixelsize # [arcsec]
@@ -184,11 +176,10 @@ class GalaxyClusters:
if ii % 50 == 0:
logger.info("[%d/%d] %.1f%% ..." % (ii, num, 100*ii/num))
z0, M0 = row.z, row.mass
- age0 = self.cosmo.age(z0)
- zmax = self.cosmo.redshift(age0 - self.tau_merger)
+ age0 = cosmo.age(z0)
+ zmax = cosmo.redshift(age0 - self.tau_merger)
clform = ClusterFormation(M0=M0, z0=z0, zmax=zmax,
ratio_major=self.ratio_major,
- cosmo=self.cosmo,
merger_mass_min=self.merger_mass_min)
clform.simulate_mergertree(main_only=True)
mmev = clform.last_major_merger