aboutsummaryrefslogtreecommitdiffstats
path: root/fg21sim/extragalactic
diff options
context:
space:
mode:
authorAaron LI <aaronly.me@outlook.com>2017-01-08 17:03:01 +0800
committerAaron LI <aly@aaronly.me>2017-06-01 16:33:39 +0800
commit797c3a59445f8e4e56748557812960c151b7f9bc (patch)
tree5c9a1ef719c8dd0b2be9dbc1df98e7347a3e3e21 /fg21sim/extragalactic
parentadc328067f34bebb81b0ce3eaab471a69d4d44f8 (diff)
downloadfg21sim-797c3a59445f8e4e56748557812960c151b7f9bc.tar.bz2
Add clusters/formation.py: Simulate cluster merging history
Diffstat (limited to 'fg21sim/extragalactic')
-rw-r--r--fg21sim/extragalactic/clusters/formation.py230
1 files changed, 230 insertions, 0 deletions
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 <liweitianux@live.com>
+# 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)
+ )