aboutsummaryrefslogtreecommitdiffstats
diff options
context:
space:
mode:
authorAaron LI <aly@aaronly.me>2017-07-22 14:38:37 +0800
committerAaron LI <aly@aaronly.me>2017-07-22 14:38:37 +0800
commit86d29755965a348c7c9f8bece775e8e110376693 (patch)
tree2266eebd6c1d1f3866a670575f07efe85d62430b
parentb72fedcf0a689aefa8d0659ac195fe61f3197201 (diff)
downloadfg21sim-86d29755965a348c7c9f8bece775e8e110376693.tar.bz2
clusters/halo.py: Significant cleanups
Several methods/functions have been migrated into "helper.py", while other methods/functions are obsolete. Signed-off-by: Aaron LI <aly@aaronly.me>
-rw-r--r--fg21sim/extragalactic/clusters/halo.py553
1 files changed, 2 insertions, 551 deletions
diff --git a/fg21sim/extragalactic/clusters/halo.py b/fg21sim/extragalactic/clusters/halo.py
index c78c2c3..958dd6c 100644
--- a/fg21sim/extragalactic/clusters/halo.py
+++ b/fg21sim/extragalactic/clusters/halo.py
@@ -20,11 +20,7 @@ References
import logging
import numpy as np
-import scipy.interpolate
-import scipy.integrate
-import scipy.optimize
-from .formation import ClusterFormation
from .solver import FokkerPlanckSolver
from ...utils import cosmo
from ...utils.units import (Units as AU,
@@ -54,101 +50,10 @@ class RadioHalo:
Unit: [Msun]
z0 : float
Redshift from where to simulate former merging history.
- configs : `ConfigManager`
- A `ConfigManager` instance containing default and user configurations.
- For more details, see the example configuration specifications.
-
- Attributes
- ----------
- mec : float
- Unit for electron momentum (p): mec = m_e * c, p = gamma * mec,
- therefore value of p is the Lorentz factor.
- mtree : `~MergerTree`
- Merging history of this cluster.
"""
- # Unit for electron momentum (p), thus its value is the Lorentz factor
- mec = AU.mec # [g cm / s]
- # Merger tree (i.e., merging history) of this cluster
- mtree = None
-
- def __init__(self, M0, z0, configs):
- self.M0 = M0 # [Msun]
+ def __init__(self, M0, z0):
+ self.M0 = M0
self.z0 = z0
- self.configs = configs
- self._set_configs()
-
- def _set_configs(self):
- """
- Set up the necessary class attributes according to the configs.
- """
- comp = "extragalactic/halos"
- self.zmax = self.configs.getn(comp+"/zmax")
- self.zbinsize = self.configs.getn(comp+"/zbinsize")
- # Mass threshold of the sub-cluster for a significant merger
- self.merger_mass_th = self.configs.getn(comp+"/merger_mass_th")
- self.radius = self.configs.getn(comp+"/radius")
- self.b_mean = self.configs.getn(comp+"/b_mean")
- self.b_index = self.configs.getn(comp+"/b_index")
- self.eta_t = self.configs.getn(comp+"/eta_t")
- self.eta_e = self.configs.getn(comp+"/eta_e")
- self.pmin = self.configs.getn(comp+"/pmin")
- self.pmax = self.configs.getn(comp+"/pmax")
- self.pgrid_num = self.configs.getn(comp+"/pgrid_num")
- 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")
- logger.info("Loaded and set up configurations")
-
- @property
- def zgrid(self):
- """
- Redshift grid between 0 and ``self.zmax``.
- """
- return np.arange(start=0.0, stop=self.zmax+self.zbinsize,
- step=self.zbinsize)
-
- @property
- def pgrid(self):
- """
- Electron momentum values of the adopted logarithmic grid
- used for solving the Fokker-Planck equation.
- """
- grid = np.logspace(np.log10(self.pmin), np.log10(self.pmax),
- num=self.pgrid_num)
- return grid
-
- @property
- def magnetic_field(self):
- """
- Calculate the mean magnetic field of this cluster from the
- scaling relation between magnetic field and cluster mass.
-
- Unit: [uG]
-
- Reference: Ref.[3],Eq.(1)
- """
- M_mean = 1.6e15 # [Msun]
- B = self.b_mean * (self.M0/M_mean) ** self.b_index
- return B
-
- def simulate_mergertree(self):
- """
- Simulate the merging history of the cluster using the extended
- Press-Schechter formalism.
-
- Attributes
- ----------
- mtree : `~MergerTree`
- Generated merger tree of this cluster.
-
- Returns
- -------
- mtree : `~MergerTree`
- Generated merger tree of this cluster.
- """
- self.formation = ClusterFormation(self.M0, self.z0, self.configs)
- self.mtree = self.formation.simulate_mergertree()
- return self.mtree
def calc_electron_spectrum(self, zbegin=None, zend=None, n0_e=None):
"""
@@ -204,261 +109,6 @@ class RadioHalo:
n_e = fpsolver.solve(u0=n0_e, tstart=tstart, tstop=tstop)
return (p, n_e)
- def kT_mass(self, mass):
- """
- Estimate the cluster ICM temperature from its mass by assuming
- an (observed) temperature-mass relation.
-
- TODO: upgrade this M-T relation.
-
- Parameters
- ----------
- mass : float
- Mass (unit: Msun) of the cluster
-
- Returns
- -------
- kT : float
- Temperature of the ICM (unit: keV)
-
- References
- ----------
- [1] Nevalainen et al. 2000, ApJ, 532, 694
- Ettori et al, 2013, Space Science Review, 177, 119-154
- NOTE: H0 = 50 * h50 [km/s/Mpc]
- """
- kT = 10 * (mass/1.23e15) ** (1/1.79) # [keV]
- return kT
-
- def _radius_virial(self, mass, z=0.0):
- """
- Calculate the virial radius of a cluster.
-
- Parameters
- ----------
- mass : float
- Mass (unit: Msun) of the cluster
- z : float
- Redshift
-
- Returns
- -------
- Rvir : float
- Virial radius (unit: kpc) of the cluster at given redshift
- """
- 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
-
- def _radius_stripping(self, mass, M_main, z):
- """
- Calculate the stripping radius of the sub-cluster at which
- equipartition between static and ram pressure is established,
- and the stripping is efficient outside this stripping radius.
-
- Note that the value of the stripping radius obtained would
- give the *mean value* of the actual stripping radius during
- a merger.
-
- Parameters
- ----------
- mass : float
- The mass (unit: Msun) of the sub-cluster.
- M_main : float
- The mass (unit: Msun) of the main cluster.
- z : float
- Redshift
-
- Returns
- -------
- rs : float
- The stripping radius of the sub-cluster.
- Unit: kpc
-
- References
- ----------
- [1] Cassano & Brunetti 2005, MNRAS, 357, 1313
- http://adsabs.harvard.edu/abs/2005MNRAS.357.1313C
- Eq.(11)
- """
- vi = self._velocity_impact(M_main, mass, z) * 1e5 # [cm/s]
- kT = self.kT_mass(mass) * AUC.keV2erg # [erg]
- coef = kT / (AC.mu * AC.u * vi**2) # dimensionless
- rho_avg = self._density_average(M_main, z) # [g/cm^3]
-
- def equation(r):
- return coef * self.density_profile(r, mass, z) / rho_avg - 1
-
- r_vir = self._radius_virial(mass, z) # [kpc]
- rs = scipy.optimize.brentq(equation, a=0, b=r_vir) # [kpc]
- return rs
-
- def _density_average(self, mass, z=0.0):
- """
- Average density of the cluster ICM.
-
- Returns
- -------
- rho : float
- Average ICM density (unit: g/cm^3)
-
- References
- ----------
- [1] Cassano & Brunetti 2005, MNRAS, 357, 1313
- http://adsabs.harvard.edu/abs/2005MNRAS.357.1313C
- Eq.(12)
- """
- 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]
- return rho
-
- def density_profile(self, r, mass, z):
- """
- ICM (baryon) density profile, assuming the beta model.
-
- Parameters
- ----------
- r : float
- Radius (unit: kpc) where to calculate the density
- mass : float
- Cluster mass (unit: Msun)
- z : float
- Redshift
-
- Returns
- -------
- rho_r : float
- Density at the specified radius (unit: g/cm^3)
-
- References
- ----------
- [1] Cassano & Brunetti 2005, MNRAS, 357, 1313
- http://adsabs.harvard.edu/abs/2005MNRAS.357.1313C
- Eq.(13)
- """
- 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]
- rc = self._beta_rc(Rv)
- beta = self._beta_beta()
- norm = self._beta_norm(M_ICM, beta, rc, Rv) # [g/cm^3]
- rho_r = norm * (1 + (r/rc)**2) ** (-3*beta/2) # [g/cm^3]
- return rho_r
-
- @staticmethod
- def _beta_rc(r_vir):
- """
- Core radius of the beta model for the ICM density profile.
-
- TODO: upgrade this!
- """
- return 0.1*r_vir
-
- @staticmethod
- def _beta_beta():
- """
- Beta value of the beta model for the ICM density profile.
-
- TODO: upgrade this!
- """
- return 0.8
-
- @staticmethod
- def _beta_norm(mass, beta, rc, r_vir):
- """
- Calculate the normalization of the beta model for the ICM
- density profile.
-
- Parameters
- ----------
- mass : float
- The mass (unit: g) of ICM
- beta : float
- Beta value of the assumed beta profile
- rc : float
- Core radius (unit: cm) of the assumed beta profile
- r_vir : float
- The virial radius (unit: cm) of the cluster
-
- Returns
- -------
- norm : float
- Normalization of the beta model (unit: g/cm^3)
-
- References
- ----------
- [1] Cassano & Brunetti 2005, MNRAS, 357, 1313
- http://adsabs.harvard.edu/abs/2005MNRAS.357.1313C
- Eq.(14)
- """
- integration = scipy.integrate.quad(
- lambda r: r*r * (1+(r/rc)**2) ** (-3*beta/2),
- 0, r_vir)[0]
- norm = mass / (4*np.pi * integration) # [g/cm^3]
- return norm
-
- def _velocity_impact(self, M_main, M_sub, z=0.0):
- """
- Calculate the relative impact velocity between the two merging
- clusters when they are at a distance of virial radius.
-
- Parameters
- ----------
- M_main : float
- Mass of the main cluster (unit: Msun)
- M_sub : float
- Mass of the sub cluster (unit: Msun)
- z : float
- Redshift
-
- Returns
- -------
- vi : float
- Relative impact velocity (unit: km/s)
-
- References
- ----------
- [1] Cassano & Brunetti 2005, MNRAS, 357, 1313
- http://adsabs.harvard.edu/abs/2005MNRAS.357.1313C
- Eq.(9)
- """
- eta_v = 4 * (1 + M_main/M_sub) ** (1/3)
- R_vir = self._radius_virial(M_main, z) * AUC.kpc2cm # [cm]
- vi = np.sqrt(2*AC.G * (1-1/eta_v) *
- (M_main+M_sub)*AUC.Msun2g / R_vir) # [cm/s]
- vi /= AUC.km2cm # [km/s]
- return vi
-
- def _time_crossing(self, M_main, M_sub, z):
- """
- Calculate the crossing time of the sub-cluster during a merger.
-
- Parameters
- ----------
- M_main : float
- Mass of the main cluster (unit: Msun)
- M_sub : float
- Mass of the sub cluster (unit: Msun)
- z : float
- Redshift where the merger occurs.
-
- Returns
- -------
- time : float
- Crossing time (unit: Gyr)
- """
- R_vir = self._radius_virial(M_main, z) # [kpc]
- vi = self._velocity_impact(M_main, M_sub, z) # [km/s]
- # Unit conversion coefficient: [s kpc/km] => [Gyr]
- uconv = AUC.kpc2km * AUC.s2Gyr
- time = uconv * R_vir / vi # [Gyr]
- return time
-
def _z_end(self, z_begin, time):
"""
Calculate the ending redshift from ``z_begin`` after elapsing
@@ -479,165 +129,6 @@ class RadioHalo:
z_end = cosmo.redshift(t_end)
return z_end
- @property
- def merger_events(self):
- """
- Trace only the main cluster, and filter out the significant
- merger events.
-
- Returns
- -------
- mevents : list[dict]
- List of dictionaries that records all the merger events
- of the main cluster.
- NOTE:
- The merger events are ordered by increasing redshifts.
- """
- events = []
- tree = self.mtree
- while tree:
- if (tree.main and tree.sub and
- tree.sub.data["mass"] >= self.merger_mass_th and
- tree.main.data["z"] <= self.zmax):
- events.append({
- "z": tree.main.data["z"],
- "age": tree.main.data["age"],
- "M_main": tree.main.data["mass"],
- "M_sub": tree.sub.data["mass"],
- })
- tree = tree.main
- return events
-
- def _coef_acceleration(self, z):
- """
- Calculate the electron-acceleration coefficient at arbitrary
- redshift.
-
- Parameters
- ----------
- z : float
- Redshift where to calculate the acceleration coefficient.
-
- Returns
- -------
- chi : float
- The calculated electron-acceleration coefficient.
- (unit: Gyr^-1)
-
- XXX/NOTE
- --------
- This coefficient may be very small and even zero, then the
- diffusion coefficient of the Fokker-Planck equation is thus
- very small and even zero, which cause problems for calculating
- some quantities (e.g., w(x), C(x)) and wrong/invalid results.
- To avoid these problems, force the minimal value of this
- coefficient to be 1/(10*t0), which t0 is the present-day age
- of the universe.
- """
- zgrid, chigrid = self._chi_data
- # XXX: force a minimal value instead of zero or too small
- chi_min = 1 / (10 * cosmo.age0)
-
- try:
- zi = np.where(z < zgrid)[0][0]
- slope = (chigrid[zi]-chigrid[zi-1]) / (zgrid[zi]-zgrid[zi-1])
- chi = (z-zgrid[zi]) * slope + chigrid[zi]
- except IndexError:
- # Given z >= zgrid[-1]
- chi = 0.0
-
- if chi > chi_min:
- return chi
- else:
- return chi_min
-
- @property
- def _chi_data(self):
- """
- Returns
- -------
- zgrid : 1D `~numpy.ndarray`
- Redshift positions where the chi values are calculated.
- Same as the attribute ``self.zgrid``
- chigrid : 1D `~numpy.ndarray`
- Values of acceleration coefficients at each redshift.
- """
- if hasattr(self, "_chi_data_"):
- zgrid, chigrid = self._chi_data_
- else:
- # Order the merger events by decreasing redshifts
- mevents = list(reversed(self.merger_events))
- chi_z = [self._chi_at_zidx(zidx, mevents)
- for zidx in range(len(mevents))]
- zgrid = self.zgrid
- chigrid = np.zeros(zgrid.shape)
- for (chi_, zbegin, zend) in chi_z:
- # NOTE: zbegin > zend
- mask = (zgrid <= zbegin) & (zgrid >= zend)
- chigrid[mask] += chi_
- self._chi_data_ = (zgrid, chigrid)
- return (zgrid, chigrid)
-
- def _chi_at_zidx(self, zidx, mevents):
- """
- Calculate electron-acceleration coefficient at the specified
- merger event which is specified with a redshift index.
-
- Parameters
- ----------
- zidx : int
- Index of the redshift where to calculate the coefficient.
- mevents : list[dict]
- List of dictionaries that records all the merger events
- of the main cluster.
- NOTE:
- The merger events should be ordered by increasing time
- (or decreasing redshifts).
-
- Returns
- -------
- chi : float
- The calculated electron-acceleration coefficient.
- Unit: [Gyr^-1]
- zbegin : float
- The redshift when the merger begins
- zend : float
- The redshift when the merger ends
- NOTE: zbegin > zend
-
- References: Ref.[1],Eq.(40)
- """
- redshifts = np.array([ev["z"] for ev in mevents])
- zbegin = mevents[zidx]["z"]
- M_main = mevents[zidx]["M_main"]
- M_sub = mevents[zidx]["M_sub"]
- t_crossing = self._time_crossing(M_main, M_sub, zbegin) # [Gyr]
- zend = self._z_end(zbegin, t_crossing)
- try:
- zend_idx = np.where(redshifts < zend)[0][0]
- except IndexError:
- # Specified redshift already the last/smallest one
- zend_idx = zidx + 1
- #
- coef = 2.23e-16 * self.eta_t / (self.radius/500)**3 # [s^-1]
- coef *= AUC.Gyr2s # [Gyr^-1]
- chi = 0.0
- for ev in mevents[zidx:zend_idx]:
- M_main = ev["M_main"]
- M_sub = ev["M_sub"]
- z = ev["z"]
- R_vir = self._radius_virial(M_main, z)
- rs = self._radius_stripping(M_sub, M_main, z)
- kT = self.kT_mass(M_main)
- term1 = ((M_main+M_sub)/2e15 * (2.6e3/R_vir)) ** (3/2)
- term2 = (rs/500)**2 / np.sqrt(kT/7)
- if rs <= self.radius:
- term3 = 1.0
- else:
- term3 = (self.radius/rs) ** 2
- chi += coef * term1 * term2 * term3
- return (chi, zbegin, zend)
-
def fp_injection(self, p, t=None):
"""
Electron injection term for the Fokker-Planck equation.
@@ -776,47 +267,7 @@ class RadioHalo:
((self.magnetic_field/3.2)**2 + (1+z)**4))
return dpdt
- @property
- def e_thermal(self):
- """
- Calculate the thermal energy density of the ICM.
-
- Returns
- -------
- e_th : float
- Energy density of the ICM (unit: erg/cm^3)
"""
- mass = self.M0
- 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]
- Rv = self._radius_virial(mass, self.z0) * AUC.kpc2cm # [cm]
- V = (4*np.pi / 3) * Rv**3 # [cm^3]
- e_th = E_th / V # [erg/cm^3]
- return e_th
-
- def _n_thermal(self, mass, z=0.0):
- """
- Calculate the number density of the ICM thermal plasma.
- Parameters
----------
- mass : float
- Mass of the cluster at the redshift
- Unit: [Msun]
- z : float
- Redshift
-
- Returns
- -------
- n_th : float
- Number density of the ICM
- Unit: [cm^-3]
"""
- 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]
- n_th = N / V # [cm^-3]
- return n_th