diff options
Diffstat (limited to 'fg21sim/extragalactic/clusters')
| -rw-r--r-- | fg21sim/extragalactic/clusters/halo.py | 553 |
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 |