diff options
Diffstat (limited to 'fg21sim')
| -rw-r--r-- | fg21sim/configs/20-extragalactic.conf.spec | 13 | ||||
| -rw-r--r-- | fg21sim/extragalactic/clusters/psformalism.py | 230 |
2 files changed, 241 insertions, 2 deletions
diff --git a/fg21sim/configs/20-extragalactic.conf.spec b/fg21sim/configs/20-extragalactic.conf.spec index a8f3510..515950e 100644 --- a/fg21sim/configs/20-extragalactic.conf.spec +++ b/fg21sim/configs/20-extragalactic.conf.spec @@ -20,11 +20,20 @@ # The configurations in this ``[[clusters]]`` section may also be # used by the following ``[[halos]]`` section. [[clusters]] - # The clusters catalog derived from the Hubble Volume Project (CSV file) - catalog = string(default=None) + # The Press-Schechter formalism predicted halo distribution densities. + ps_data = string(default=None) + # Output the effective/inuse clusters catalog data (CSV file) catalog_outfile = string(default=None) + # The fraction of the dark matter mass in galaxy clusters. + f_darkmatter = float(default=0.8, min=0.5, max=1.0) + + # The minimum mass for clusters when to determine the galaxy clusters + # total counts and their distributions. + # Unit: [Msun] + mass_min = float(default=2e14, min=1e12) + # Minimum mass change of the main cluster to be regarded as a merger # event instead of an accretion event. # Unit: [Msun] diff --git a/fg21sim/extragalactic/clusters/psformalism.py b/fg21sim/extragalactic/clusters/psformalism.py new file mode 100644 index 0000000..4a65b8f --- /dev/null +++ b/fg21sim/extragalactic/clusters/psformalism.py @@ -0,0 +1,230 @@ +# Copyright (c) 2017 Weitian LI <weitian@aaronly.me> +# MIT license + +""" +Press-Schechter (PS) formalism + +First determine the number of clusters within a sky patch (i.e., sky +coverage) according to the cluster distribution predicted by the PS +formalism; then sampling from the PS mass function to derive the mass +and redshift for each cluster. +""" + +import logging +import random + +import numpy as np +import pandas as pd + +from ...share import CONFIGS, COSMO +from ...utils.interpolate import bilinear +from ...utils.units import UnitConversions as AUC + + +logger = logging.getLogger(__name__) + + +class PSFormalism: + """ + Press-Schechter (PS) formalism + + Simulate the clusters number and their distribution (mass and z) + within a sky patch of certain coverage. + """ + def __init__(self, configs=CONFIGS): + self.configs = configs + self._set_configs() + self._load_data() + + def _set_configs(self): + """ + Load the required configurations and set them. + """ + comp = "extragalactic/clusters" + self.datafile = self.configs.get_path(comp+"/ps_data") + self.f_darkmatter = self.configs.getn(comp+"/f_darkmatter") + self.Mmin_cluster = self.configs.getn(comp+"/mass_min") # [Msun] + self.Mmin_halo = self.Mmin_cluster * self.f_darkmatter + + def _load_data(self, filepath=None): + """ + Load dndM data and reformat into a 2D density grid together with + redshifts and masses vectors. + + Data File Description + --------------------- + z1 mass1 density1 + z1 mass2 density2 + z1 .. density3 + z2 mass1 density4 + z2 mass2 density5 + z2 .. density6 + ... + + where, + * Redshifts: 0.0 -> 3.02, even-spacing, step 0.02 + * Mass: unit 1e12 -> 9.12e15 [Msun], log-even (dark matter) + * Density: [number]/dVc/dM + with, + - dVc: differential comvoing volume, [Mpc^3]/[sr]/[unit redshift] + """ + if filepath is None: + filepath = self.datafile + data = np.loadtxt(filepath) + redshifts = data[:, 0] + masses = data[:, 1] + densities = data[:, 2] + + redshifts = np.array(list(set(redshifts))) + redshifts.sort() + masses = np.array(list(set(masses))) + masses.sort() + densities = densities.reshape((len(redshifts), len(masses))) + + logger.info("Loaded PS data from file: %s" % filepath) + logger.info("Number of redshift bins: %d" % len(redshifts)) + logger.info("Number of mass bins: %d" % len(masses)) + self.redshifts = redshifts + self.masses = masses + self.densities = densities + + @staticmethod + def delta(x, logeven=False): + """ + Calculate the delta values for each element of a vector, + assuming they are evenly or log-evenly distributed, + with extrapolating. + """ + x = np.asarray(x) + if logeven: + x = np.log(x) + step = x[1] - x[0] + x1 = np.concatenate([[x[0]-step], x[:-1]]) + x2 = np.concatenate([x[1:], [x[-1]+step]]) + dx = (x2 - x1) * 0.5 + if logeven: + dx = np.exp(dx) + return dx + + @property + def number_grid(self): + """ + Calculate the number distribution w.r.t. redshift, mass, and + unit coverage [sr] from the density distribution. + """ + dz = self.delta(self.redshifts) + dM = self.delta(self.masses) + dMgrip, dzgrip = np.meshgrid(dM, dz) + Mgrip, zgrip = np.meshgrid(self.masses, self.redshifts) + dVcgrip = COSMO.differential_comoving_volume(zgrip).value # [Mpc^3/sr] + numgrid = self.densities * dVcgrip * dzgrip * dMgrip + return numgrid + + def calc_cluster_counts(self, coverage): + """ + Calculate the total number of clusters (>= minimum mass) within + the FoV coverage according to the number density distribution + (e.g., predicted by the Press-Schechter mass function) + + Parameters + ---------- + coverage : float + The coverage of the sky patch within which to determine the + total number of clusters. + Unit: [deg^2] + + Returns + ------- + counts : int + The total number of clusters within the sky patch. + + Attributes + ---------- + counts + """ + logger.info("Determine the total number of clusters within " + "sky patch of coverage %.1f [deg^2]" % coverage) + coverage *= AUC.deg2rad**2 # [deg^2] -> [rad^2] = [sr] + midx = (self.masses >= self.Mmin_halo) + numgrid = self.number_grid + counts = np.sum(numgrid[:, midx]) * coverage + self.counts = int(np.round(counts)) + logger.info("Total number of clusters: %d" % self.counts) + return self.counts + + def sample_z_m(self, counts=None): + """ + Randomly generate the requested number of pairs of (z, M) following + the specified number distribution. + + Parameters + ---------- + counts : int, optional + The number of (z, mass) pairs to be sampled. + If not specified, then default to ``self.counts`` + + Returns + ------- + df : `~pandas.DataFrame` + A Pandas data frame with 2 columns, i.e., ``z`` and ``mass``. + comment : list[str] + Comments to the above data frame. + + Attributes + ---------- + clusters : df + clusters_comment : comment + """ + if counts is None: + counts = self.counts + logger.info("Sampling (z, mass) pairs for %d clusters ..." % counts) + + redshifts = self.redshifts + masses = self.masses + zmin = redshifts.min() + zmax = redshifts.max() + Mmax = masses.max() + midx = (masses >= self.Mmin_halo) + numgrid = self.number_grid + numgrid2 = numgrid[:, midx] + NM = numgrid2.max() + z_list = [] + M_list = [] + i = 0 + while i < counts: + z = random.uniform(zmin, zmax) + M = random.uniform(self.mass_min, Mmax) + r = random.random() + zi1 = (self.redshifts < z).sum() + zi2 = zi1 - 1 + if zi2 < 0: + zi2 += 1 + zi1 += 1 + Mi1 = (self.masses < M).sum() + Mi2 = Mi1 - 1 + if Mi2 < 0: + Mi2 += 1 + Mi1 += 1 + N = bilinear( + z, np.log(M), + p11=(redshifts[zi1], np.log(masses[Mi1]), numgrid[zi1, Mi1]), + p12=(redshifts[zi1], np.log(masses[Mi2]), numgrid[zi1, Mi2]), + p21=(redshifts[zi2], np.log(masses[Mi1]), numgrid[zi2, Mi1]), + p22=(redshifts[zi2], np.log(masses[Mi2]), numgrid[zi2, Mi2])) + if r < N/NM: + z_list.append(z) + M_list.append(M) + i += 1 + logger.info("Sampled %d (z, mass) pairs for each cluster" % counts) + + df = pd.DataFrame(np.column_stack([z_list, M_list]), + columns=["z", "mass"]) + df["mass"] /= self.f_darkmatter + comment = [ + "cluster number counts : %d" % counts, + "z : redshift", + "mass : cluster total mass [Msun]", + ] + self.clusters = df + self.clusters_comment = comment + return (df, comment) |