From 797c3a59445f8e4e56748557812960c151b7f9bc Mon Sep 17 00:00:00 2001 From: Aaron LI Date: Sun, 8 Jan 2017 17:03:01 +0800 Subject: Add clusters/formation.py: Simulate cluster merging history --- fg21sim/configs/20-extragalactic.conf.spec | 4 + fg21sim/extragalactic/clusters/formation.py | 230 ++++++++++++++++++++++++++++ 2 files changed, 234 insertions(+) create mode 100644 fg21sim/extragalactic/clusters/formation.py diff --git a/fg21sim/configs/20-extragalactic.conf.spec b/fg21sim/configs/20-extragalactic.conf.spec index de6dbeb..f9d7c88 100644 --- a/fg21sim/configs/20-extragalactic.conf.spec +++ b/fg21sim/configs/20-extragalactic.conf.spec @@ -53,6 +53,10 @@ # merger. (unit: Msun) merger_mass_th = float(default=1e13, min=1e12) + # Minimum mass change of the main-cluster to be regarded as a merger + # event rather than accretion. (unit: Msun) + merger_mass_min = float(default=1e12, min=1e10) + # Radius of the giant radio halo in clusters (unit: kpc) # XXX: currently only support a constant radius of halos radius_halo = float(default=500, min=100) diff --git a/fg21sim/extragalactic/clusters/formation.py b/fg21sim/extragalactic/clusters/formation.py new file mode 100644 index 0000000..0b8e9d9 --- /dev/null +++ b/fg21sim/extragalactic/clusters/formation.py @@ -0,0 +1,230 @@ +# Copyright (c) 2017 Weitian LI +# MIT license + +""" +Simulate cluster formation (i.e., merging history) using the extended +Press-Schechter formalism. + +References +---------- +[1] Randall, Sarazin & Ricker 2002, ApJ, 577, 579 + http://adsabs.harvard.edu/abs/2002ApJ...577..579R +[2] Cassano & Brunetti 2005, MNRAS, 357, 1313 + http://adsabs.harvard.edu/abs/2005MNRAS.357.1313C +""" + +import logging + +import numpy as np +import scipy.integrate +import scipy.special +import scipy.optimize + +from .cosmology import Cosmology +from .mergertree import MergerTree + + +logger = logging.getLogger(__name__) + + +class ClusterFormation: + """ + Simulate the cluster formation (i.e., merging history) using the extended + Press-Schechter formalism by Monte Carlo methods. + + References + ---------- + [1] Randall, Sarazin & Ricker 2002, ApJ, 577, 579 + http://adsabs.harvard.edu/abs/2002ApJ...577..579R + [2] Cassano & Brunetti 2005, MNRAS, 357, 1313 + http://adsabs.harvard.edu/abs/2005MNRAS.357.1313C + + Parameters + ---------- + M0 : float + Present-day (z=0) mass (unit: Msun) of the cluster. + configs : `ConfigManager` + A `ConfigManager` instance containing default and user configurations. + For more details, see the example configuration specifications. + + Attributes + ---------- + cosmo : `~Cosmology` + Adopted cosmological model with custom utility functions. + mtree : `~MergerTree` + Merging history of this cluster. + """ + def __init__(self, M0, configs): + self.M0 = M0 # [Msun] + self.configs = configs + self._set_configs() + + def _set_configs(self): + """ + Set up the necessary class attributes according to the configs. + """ + comp = "extragalactic/halos" + # Minimum mass change (unit: Msun) of the main-cluster for a merger + self.merger_mass_min = self.configs.getn(comp+"/merger_mass_min") + # Cosmology model + self.H0 = self.configs.getn("cosmology/H0") + self.OmegaM0 = self.configs.getn("cosmology/OmegaM0") + self.sigma8 = self.configs.getn("cosmology/sigma8") + self.cosmo = Cosmology(H0=self.H0, Om0=self.OmegaM0, + sigma8=self.sigma8) + logger.info("Loaded and set up configurations") + + @property + def sigma_index(self): + """ + The power-law spectral index assumed for the following density + perturbations sigma(M). + + References: Ref.[1],Eq.(2) + """ + n = -7/5 + alpha = (n+3) / 6 + return alpha + + def f_sigma(self, mass): + """ + Current rms density fluctuations within a sphere of specified + mass (unit: Msun). + + It is generally sufficient to consider a power-law spectrum of + density perturbations, which is consistent with the CDM models. + + References: Ref.[1],Eq.(2) + """ + alpha = self.sigma_index + sigma = self.cosmo.sigma8 * (mass / self.cosmo.M8) ** (-alpha) + return sigma + + def f_delta_c(self, z): + """ + w = delta_c(z) is the critical linear overdensity for a region + to collapse at redshift z. + + This is a monotone decreasing function. + + References: Ref.[1],App.A,Eq.(A1) + """ + return self.cosmo.overdensity_crit(z) + + def f_dw_max(self, mass): + """ + Calculate the allowed maximum step size for tracing cluster + formation, therefore, the adopted step size is chosen to be half + of this maximum value. + + dw^2 ~< abs(d(ln(sigma(M)^2)) / d(ln(M))) * (dMc / M) * sigma(M)^2 + = 2 * alpha * sigma(M)^2 * dMc / M + + References: Ref.[1],Sec.(3.1),Para.(1) + """ + alpha = self.sigma_index + dMc = self.merger_mass_min + return np.sqrt(2 * alpha * self.f_sigma(mass)**2 * dMc / mass) + + def calc_z(self, delta_c): + """ + Solve the redshift from the specified delta_c (a.k.a. w). + """ + z = scipy.optimize.newton( + lambda x: self.f_delta_c(x) - delta_c, + x0=0, tol=1e-5) + return z + + def calc_mass(self, S): + """ + Calculate the mass corresponding to the given S. + + S = sigma(M)^2 + + References: Ref.[1],Sec.(3) + """ + alpha = self.sigma_index + mass = self.cosmo.M8 * (S / self.cosmo.sigma8**2)**(-1/(2*alpha)) + return mass + + @staticmethod + def cdf_K(dS, dw): + """ + The cumulative probability distribution function of sub-cluster + masses. + + References: Ref.[1],Eq.(5) + """ + p = scipy.special.erfc(dw / np.sqrt(2*dS)) + return p + + @staticmethod + def cdf_K_inv(p, dw): + """ + Inverse function of the above ``cdf_K()``. + """ + dS = 0.5 * (dw / scipy.special.erfcinv(p))**2 + return dS + + def gen_dS(self, dw, size=None): + """ + Randomly generate values of dS by sampling the CDF ``cdf_K()``. + """ + r = np.random.uniform(size=size) + dS = self.cdf_K_inv(r, dw) + return dS + + def simulate_mergertree(self): + """ + Simulate the merger tree of this cluster by tracing its formation + using the PS formalism. + + References: Ref.[1],Sec.(3.1) + """ + self.mtree = self._trace_formation(self.M0, dMc=self.merger_mass_min) + return self.mtree + + def _trace_formation(self, M, dMc, _z=None): + """ + Recursively trace the cluster formation and thus simulate its + merger tree. + """ + z = 0.0 if _z is None else _z + node_data = {"mass": M, "z": z, "age": self.cosmo.age(z)} + + if M <= dMc: + # Stop the trace + return MergerTree(data=node_data) + + # Trace the formation by simulate a merger/accretion event + # Notation: progenitor (*1) -> current (*2) + + # Current properties + w2 = self.f_delta_c(z=z) + S2 = self.f_sigma(M) ** 2 + dw = 0.5 * self.f_dw_max(M, dMc) + dS = self.gen_dS(dw) + # Progenitor properties + z1 = self.calc_z(w2 + dw) + S1 = S2 + dS + M1 = self.calc_mass(S1) + dM = M - M1 + + M_min = min(M1, dM) + if M_min <= dMc: + # Accretion + M_new = M - M_min + return MergerTree( + data=node_data, + main=self._trace_formation(M_new, dMc=dMc, _z=z1), + sub=None + ) + else: + # Merger event + M_main = max(M1, dM) + M_sub = M_min + return MergerTree( + data=node_data, + main=self._trace_formation(M_main, dMc=dMc, _z=z1), + sub=self._trace_formation(M_sub, dMc=dMc, _z=z1) + ) -- cgit v1.2.2