diff options
Diffstat (limited to 'fg21sim')
| -rw-r--r-- | fg21sim/configs/20-extragalactic.conf.spec | 2 | ||||
| -rw-r--r-- | fg21sim/extragalactic/__init__.py | 4 | ||||
| -rw-r--r-- | fg21sim/extragalactic/clusters/clusters.py | 732 |
3 files changed, 3 insertions, 735 deletions
diff --git a/fg21sim/configs/20-extragalactic.conf.spec b/fg21sim/configs/20-extragalactic.conf.spec index ddbabe6..ee58c8d 100644 --- a/fg21sim/configs/20-extragalactic.conf.spec +++ b/fg21sim/configs/20-extragalactic.conf.spec @@ -107,7 +107,7 @@ [[pointsources]] # Whether save this point source catelogue to disk save = boolean(default=True) - # Output directory to save the simulated catelogues + # Output directory to save the simulated catalog output_dir = string(default="PS_tables") # PS components to be simulated pscomponents = string_list(default=list()) diff --git a/fg21sim/extragalactic/__init__.py b/fg21sim/extragalactic/__init__.py index d4ba7a8..0935e03 100644 --- a/fg21sim/extragalactic/__init__.py +++ b/fg21sim/extragalactic/__init__.py @@ -1,5 +1,5 @@ -# Copyright (c) 2016 Weitian LI <liweitianux@live.com> +# Copyright (c) 2016-2017 Weitian LI <weitian@aaronly.me> # MIT license -from .clusters.clusters import GalaxyClusters +from .clusters import GalaxyClusters from .pointsources import PointSources diff --git a/fg21sim/extragalactic/clusters/clusters.py b/fg21sim/extragalactic/clusters/clusters.py deleted file mode 100644 index 94f00be..0000000 --- a/fg21sim/extragalactic/clusters/clusters.py +++ /dev/null @@ -1,732 +0,0 @@ -# Copyright (c) 2016-2017 Weitian LI <weitian@aaronly.me> -# MIT license - -""" -Simulation of the radio emissions from clusters of galaxies. - -XXX/TODO --------- -* Only consider radio *halos*, with many simplifications! -* Support radio *relics* simulations ... -""" - -import os -import logging -from datetime import datetime, timezone - -import numpy as np -import astropy.units as au -from astropy.io import fits -from astropy.cosmology import FlatLambdaCDM -import pandas as pd - -from ...sky import get_sky -from ...utils.wcs import make_wcs -from ...utils.random import spherical_uniform -from ...utils.convert import Fnu_to_Tb_fast -from ...utils.grid import make_ellipse -from ...utils.units import UnitConversions as AUC - - -logger = logging.getLogger(__name__) - - -class GalaxyClusters: - """ - Simulate the radio emissions from the clusters of galaxies, which - host radio halos and relics (currently not considered). - - The simulation follows the method adopted by [Jelic2008]_, which uses - the *ΛCDM deep wedge cluster catalog* derived from the *Hubble Volume - Project* [HVP]_, [Evard2002]_. - - Every radio cluster is simulated as an *ellipse* of *uniform brightness* - on a local coordinate grid with relatively higher resolution compared - to the output HEALPix map, which is then mapped to the output HEALPix - map by down-sampling, i.e., in a similar way as the simulations of SNRs. - - TODO: ??? - - Parameters - ---------- - configs : `ConfigManager` - A `ConfigManager` instance containing default and user configurations. - For more details, see the example configuration specifications. - - Attributes - ---------- - TODO: ??? - - NOTE - ---- - Currently, only radio *halos* are considered with many simplifications. - Radio *relics* simulations need more investigations ... - - References - ---------- - .. [Jelic2008] - Jelić, V. et al., - "Foreground simulations for the LOFAR-epoch of reionization experiment", - 2008, MNRAS, 389, 1319-1335, - http://adsabs.harvard.edu/abs/2008MNRAS.389.1319J - - .. [Evard2002] - Evard, A. E. et al., - "Galaxy Clusters in Hubble Volume Simulations: Cosmological Constraints - from Sky Survey Populations", - 2002, ApJ, 573, 7-36, - http://adsabs.harvard.edu/abs/2002ApJ...573....7E - - .. [EnBlin2002] - Enßlin, T. A. & Röttgering, H., - "The radio luminosity function of cluster radio halos", - 2002, A&A, 396, 83-89, - http://adsabs.harvard.edu/abs/2002A%26A...396...83E - - .. [Reiprich2002] - Reiprich, Thomas H. & Böhringer, Hans, - "The Mass Function of an X-Ray Flux-limited Sample of Galaxy Clusters", - 2002, ApJ, 567, 716-740, - http://adsabs.harvard.edu/abs/2002ApJ...567..716R - """ - # Component name - name = "extraglalactic clusters of galaxies" - - def __init__(self, configs): - self.configs = configs - self.sky = get_sky(configs) - self._set_configs() - - def _set_configs(self): - """Load the configs and set the corresponding class attributes.""" - comp = "extragalactic/clusters" - self.catalog_path = self.configs.get_path(comp+"/catalog") - self.catalog_outfile = self.configs.get_path(comp+"/catalog_outfile") - self.halo_fraction = self.configs.getn(comp+"/halo_fraction") - self.resolution = self.configs.getn(comp+"/resolution") # [ arcmin ] - self.prefix = self.configs.getn(comp+"/prefix") - self.save = self.configs.getn(comp+"/save") - self.output_dir = self.configs.get_path(comp+"/output_dir") - # - self.filename_pattern = self.configs.getn("output/filename_pattern") - self.use_float = self.configs.getn("output/use_float") - self.checksum = self.configs.getn("output/checksum") - self.clobber = self.configs.getn("output/clobber") - self.freq_unit = au.Unit(self.configs.getn("frequency/unit")) - # Cosmology model - self.H0 = self.configs.getn("cosmology/H0") - self.OmegaM0 = self.configs.getn("cosmology/OmegaM0") - self.cosmo = FlatLambdaCDM(H0=self.H0, Om0=self.OmegaM0) - # - logger.info("Loaded and set up configurations") - - def _load_catalog(self): - """Load the cluster catalog data and set up its properties.""" - self.catalog = pd.read_csv(self.catalog_path) - nrow, ncol = self.catalog.shape - logger.info("Loaded clusters catalog data from: {0}".format( - self.catalog_path)) - logger.info("Clusters catalog data: {0} objects, {1} columns".format( - nrow, ncol)) - # Set the properties for this catalog - self.catalog_prop = { - "omega_m": 0.3, - "omega_lambda": 0.7, - # Dimensionless Hubble constant - "h": 0.7, - "sigma8": 0.9, - # Number of particles - "n_particle": 1e9, - # Particle mass of the simulation [ h^-1 ] - "m_particle": 2.25e12 * au.solMass, - # Cube side length [ h^-1 ] - "l_side": 3000.0 * au.Mpc, - # Overdensity adopted to derive the clusters - "overdensity": 200, - # Sky coverage - "coverage": 10*au.deg * 10*au.deg - } - # Units for the catalog columns (also be populated by other methods) - self.units = {} - - def _save_catalog_inuse(self): - """Save the effective/inuse clusters catalog data to a CSV file.""" - if self.catalog_outfile is None: - logger.warning("Catalog output file not set; skip saving.") - return - # Create directory if necessary - dirname = os.path.dirname(self.catalog_outfile) - if not os.path.exists(dirname): - os.mkdir(dirname) - logger.info("Created directory: {0}".format(dirname)) - # Save catalog data - if os.path.exists(self.catalog_outfile): - if self.clobber: - logger.warning("Remove existing catalog file: {0}".format( - self.catalog_outfile)) - os.remove(self.catalog_outfile) - else: - raise OSError("Output file already exists: {0}".format( - self.catalog_outfile)) - self.catalog.to_csv(self.catalog_outfile, header=True, index=False) - logger.info("Save clusters catalog in use to: {0}".format( - self.catalog_outfile)) - - def _process_catalog(self): - """Process the catalog to prepare for the simulation.""" - # Dimensionless Hubble parameter adopted in THIS simulation - h = self.H0 / 100 - logger.info("Adopted dimensionless Hubble parameter: {0}".format(h)) - # Cluster masses, unit: solMass (NOTE: h dependence) - self.catalog["mass"] = (self.catalog["m"] * - self.catalog_prop["m_particle"].value / h) - self.units["mass"] = au.solMass - logger.info("Catalog: calculated cluster masses") - # Cluster distances from the observer, unit: Mpc - dist = ((self.catalog["x"]**2 + self.catalog["y"]**2 + - self.catalog["z"]**2) ** 0.5 * - self.catalog_prop["l_side"].value / h) - self.catalog["distance"] = dist - self.units["distance"] = au.Mpc - logger.info("Catalog: calculated cluster distances") - # Drop unnecessary columns to save memory - columns_drop = ["m", "sigma", "ip", "x", "y", "z", "vx", "vy", "vz"] - self.catalog.drop(columns_drop, axis=1, inplace=True) - logger.info("Catalog: dropped unnecessary columns: {0}".format( - ", ".join(columns_drop))) - - def _scale_catalog_coverage(self): - """ - Scale the catalog to match the coverage of the simulation sky - (patch or all sky), by assuming that clusters are uniformly - distributed. Also, the radio halo fraction is also considered - to determine the final number of clusters within the simulation - sky. - """ - skyarea = self.sky.area # [ deg^2 ] - logger.info("Simulation sky coverage: %s [deg^2]" % skyarea) - logger.info("Cluster catalog sky coverage: %s [deg^2]" % - self.catalog_prop["coverage"]) - factor = float(skyarea / self.catalog_prop["coverage"]) - n0_cluster = len(self.catalog) - logger.info("Radio halo fraction in clusters: {0}".format( - self.halo_fraction)) - # Total number of radio halos within the simulation sky - N_halo = int(n0_cluster * factor * self.halo_fraction) - logger.info("Total number of radio halos within the " + - "simulation sky: {0:,}".format(N_halo)) - logger.info("Scale the catalog to match the simulation sky ...") - idx = np.round(np.random.uniform(low=0, high=n0_cluster-1, - size=N_halo)).astype(np.int) - self.catalog = self.catalog.iloc[idx, :] - self.catalog.reset_index(inplace=True) - logger.info("DONE scale the catalog to match the simulation sky") - - def _add_random_position(self): - """ - Add random positions for each cluster as columns "glon" and - "glat" to the catalog data. - - Column "glon" is the Galactic longitudes, [0, 360) (degree). - Column "glat" is the Galactic latitudes, [-90, 90] (degree). - - The positions are uniformly distributed on the spherical surface. - """ - logger.info("Randomly generating positions for each cluster ...") - num = len(self.catalog) - theta, phi = spherical_uniform(num) - glon = np.degrees(phi) - glat = 90.0 - np.degrees(theta) - self.catalog["glon"] = glon - self.catalog["glat"] = glat - logger.info("Done add random positions for each cluster") - - def _add_random_eccentricity(self): - """ - Add random eccentricities for each cluster as column - "eccentricity" to the catalog data. - - The eccentricity of a ellipse is defined as: - e = sqrt((a^2 - b^2) / a^2) = f / a - where f is the distance from the center to either focus: - f = sqrt(a^2 - b^2) - - NOTE - ---- - The eccentricities are randomly generated from a *squared* - standard normalization distribution, and with an upper limit - at 0.9, i.e., the eccentricities are [0, 0.9]. - """ - logger.info("Adding random eccentricities for each cluster ...") - num = len(self.catalog) - eccentricity = np.random.normal(size=num) ** 2 - # Replace values beyond the upper limit by sampling from valid values - ulimit = 0.9 - idx_invalid = (eccentricity > ulimit) - num_invalid = idx_invalid.sum() - eccentricity[idx_invalid] = np.random.choice( - eccentricity[~idx_invalid], size=num_invalid) - self.catalog["eccentricity"] = eccentricity - logger.info("Done add random eccentricities to catalog") - - def _add_random_rotation(self): - """ - Add random rotation angles for each cluster as column "rotation" - to the catalog data. - - The rotation angles are uniformly distributed within [0, 360). - - The rotation happens on the spherical surface, i.e., not with respect - to the line of sight, but to the Galactic frame coordinate axes. - """ - logger.info("Adding random rotation angles for each cluster ...") - num = len(self.catalog) - rotation = np.random.uniform(low=0.0, high=360.0, size=num) - self.catalog["rotation"] = rotation - self.units["rotation"] = au.deg - logger.info("Done add random rotation angles to catalog") - - def _calc_sizes(self): - """ - Calculate the virial radii for each cluster from the masses, - and then calculate the elliptical angular sizes by considering - the added random eccentricities. - - Attributes - ---------- - catalog["r_vir"] : 1D `~numpy.ndarray` - The virial radii (unit: Mpc) calculated from the cluster masses - catalog["size_major"], catalog["size_minor] : 1D `~numpy.ndarray` - The major and minor axes (unit: arcmin) of the clusters calculated - from the above virial radii and the random eccentricities. - NOTE: These major and minor axes are corresponding to the - diameter values; NOT semi-major/semi-minor axes! - Unit: arcmin - - NOTE - ---- - The elliptical major and minor axes are calculated by assuming - the equal area between the ellipse and corresponding circle. - theta2 = r_vir / distance # half angular size - pi * a * b = pi * (theta2)^2 - e = sqrt((a^2 - b^2) / a^2) # eccentricity - thus, - a = theta2 / (1-e^2)^(1/4) # semi-major axis - b = theta2 * (1-e^2)^(1/4) # semi-minor axis - """ - logger.info("Calculating the virial radii ...") - overdensity = self.catalog_prop["overdensity"] - rho_crit = self.cosmo.critical_density(self.catalog["redshift"]) - mass = self.catalog["mass"].data * self.units["mass"] - r_vir = (3 * mass / (4*np.pi * overdensity * rho_crit)) ** (1.0/3.0) - self.catalog["r_vir"] = r_vir.to(au.Mpc).value - self.units["r_vir"] = au.Mpc - logger.info("Done calculate the virial radii") - # Calculate (elliptical) angular sizes, i.e., major and minor axes - logger.info("Calculating the elliptical angular sizes ...") - distance = self.catalog["distance"].data * self.units["distance"] - # Half angular size - theta2 = (r_vir / distance).decompose().value # [rad] - theta = theta2 * 2 # angular size - # Major and minor axes (corresponding to diameter values) - size_major = theta / (1 - self.catalog["eccentricity"]**2) ** 0.25 - size_minor = theta * (1 - self.catalog["eccentricity"]**2) ** 0.25 - self.catalog["size_major"] = size_major * au.rad.to(au.arcmin) - self.catalog["size_minor"] = size_minor * au.rad.to(au.arcmin) - self.units["size"] = au.arcmin - logger.info("Done calculate the elliptical angular sizes") - - def _calc_luminosity(self): - """ - Calculate the radio luminosity (at 1.4 GHz) using empirical - scaling relations. - - First, calculate the X-ray luminosity L_X using the empirical - scaling relation between mass and X-ray luminosity. - Then, derive the radio luminosity by employing the scaling - relation between X-ray and radio luminosity. - - Attributes - ---------- - catalog["luminosity"] : 1D `~numpy.ndarray` - The luminosity density (at 1.4 GHz) of each cluster. - catalog_prop["luminosity_freq"] : `~astropy.units.Quantity` - The frequency (as an ``astropy`` quantity) where the above - luminosity derived. - units["luminosity"] : `~astropy.units.Unit` - The unit used by the above luminosity. - - XXX/TODO - -------- - The scaling relations used here may be outdated, and some of the - quantities need trick conversions, which cause much confusion. - - Investigate for *more up-to-date scaling relations*, derived with - new observation constraints. - - NOTE - ---- - - The mass in the mass-X-ray luminosity scaling relation is NOT the - cluster real mass, since [Reiprich2002]_ refer to the - *critical density* ρ_c, while the scaling relation from - Jenkins et al. (2001) requires mass refer to the - *cosmic mean mass density ρ_m = Ω_m * ρ_c, - therefore, the mass needs following conversion (which is an - approximation): - M_{R&B} ≈ M * sqrt(OmegaM0) - - The derived X-ray luminosity is for the 0.1-2.4 keV energy band. - - The X-ray-radio luminosity scaling relation adopted here is - derived at 1.4 GHz. - - [EnBlin2002]_ assumes H0 = 50 h50 km/s/Mpc, so h50 = 1. - - References - ---------- - - [Jelic2008], Eq.(13,14) - - [EnBlin2002], Eq.(1,3) - """ - # Dimensionless Hubble parameter adopted here and in the literature - h_our = self.H0 / 100 - h_others = 50.0 / 100 - # - logger.info("Calculating the radio luminosity (at 1.4 GHz) ...") - # Calculate the X-ray luminosity from mass - # NOTE: mass conversion (see also the above notes) - mass_RB = (self.catalog["mass"].data * self.units["mass"] * - self.catalog_prop["omega_m"]**0.5) - a_X = 0.449 - b_X = 1.9 - # Hubble parameter conversion factor - h_conv1 = (h_our / h_others) ** (b_X-2) - # X-ray luminosity (0.1-2.4 keV) [ erg/s ] - L_X = ((a_X * 1e45 * - (mass_RB / (1e15*au.solMass)).decompose().value ** b_X) * - h_conv1) - # Calculate the radio luminosity from X-ray luminosity - a_r = 2.78 - b_r = 1.94 - # Hubble parameter conversion factor - h_conv2 = (h_our / h_others) ** (2*b_r-2) - # Radio luminosity density (at 1.4 GHz) [ W/Hz ] - L_r = (a_r * 1e24 * (L_X / 1e45)**b_r) * h_conv2 - self.catalog["luminosity"] = L_r - self.catalog_prop["luminosity_freq"] = 1400 * self.freq_unit - self.units["luminosity"] = au.W / au.Hz - logger.info("Done calculate the radio luminosity") - - def _calc_specindex(self): - """ - Calculate the radio spectral indexes for each cluster. - - Attributes - ---------- - catalog["specindex"] : 1D `~numpy.ndarray` - The radio spectral index of each cluster. - - XXX/TODO - -------- - Currently, a common/uniform spectral index (1.2) is assumed for all - clusters, which may be improved by investigating more recent results. - """ - specindex = 1.2 - logger.info("Use same spectral index for all clusters: " - "{0}".format(specindex)) - self.catalog["specindex"] = specindex - - def _calc_Tb(self, luminosity, distance, specindex, frequency, size): - """ - Calculate the brightness temperature at requested frequency - by assuming a power-law spectral shape. - - Parameters - ---------- - luminosity : float - The luminosity density (unit: [ W/Hz ]) at the reference - frequency (i.e., `self.catalog_prop["luminosity_freq"]`). - distance : float - The luminosity distance (unit: [ Mpc ]) to the object - specindex : float - The spectral index of the power-law spectrum. - Note the *negative* sign in the formula. - frequency : float - The frequency (unit: [ MHz ]) where the brightness - temperature requested. - size : 2-float tuple - The (major, minor) axes (unit: [ deg ]). - The order of major and minor can be arbitrary. - - Returns - ------- - Tb : float - Brightness temperature at the requested frequency, unit [ K ] - - NOTE - ---- - The power-law spectral shape is assumed for *flux density* other - than the *brightness temperature*. - Therefore, the flux density at the requested frequency should first - be calculated by extrapolating the spectrum, then convert the flux - density to derive the brightness temperature. - - XXX/NOTE - -------- - The *luminosity distance* is required to calculate the flux density - from the luminosity density. - Whether the distance (i.e., ``self.catalog["distance"]``) is the - *comoving distance* ?? - Whether a conversion is required to get the *luminosity distance* ?? - """ - freq = frequency # [ MHz ] - freq_ref = self.catalog_prop["luminosity_freq"].value - Lnu = luminosity * (freq / freq_ref) ** (-specindex) # [ W/Hz ] - # Conversion coefficient: [ W/Hz/Mpc^2 ] => [ Jy ] - coef = 1.0502650403056097e-19 - Fnu = coef * Lnu / (4*np.pi * distance**2) # [ Jy ] - omega = size[0] * size[1] * AUC.arcsec2deg**2 # [ arcsec^2 ] - Tb = Fnu_to_Tb_fast(Fnu, omega, freq) - return Tb - - def _simulate_templates(self): - """ - Simulate the template sky images for each cluster, and cache - these templates within the class. - - The template images/maps have values of (or approximate) ones for - these effective pixels, excluding the pixels corresponding to the - edges of original rotated ellipse, which may have values of - significantly less than 1 due to the rotation. - - Therefore, simulating the HEALPix map of one cluster at a requested - frequency is simply multiplying the cached template image by the - calculated brightness temperature (Tb) at that frequency. - - Furthermore, the total HEALPix map of all clusters are straightforward - additions of all the maps of each cluster. - - Attributes - ---------- - templates : list - A list containing the simulated templates for each cluster/halo. - Each element is a `(idx, val)` tuple with `hpidx` the indexes - of effective image pixels and `hpval` the values of the - corresponding pixels. - e.g., ``[ (idx1, val1), (idx2, val2), ... ]`` - """ - logger.info("Simulating sky templates for each cluster/halo ...") - templates = [] - # Make sure the index is reset, therefore, the *row indexes* can be - # simply used to identify the corresponding template image. - self.catalog.reset_index(inplace=True) - # XXX/TODO: be parallel - for row in self.catalog.itertuples(): - # TODO: progress bar - gcenter = (row.glon, row.glat) # [ deg ] - radii = (int(np.ceil(row.size_major * 0.5 / self.resolution)), - int(np.ceil(row.size_minor * 0.5 / self.resolution))) - rmax = max(radii) - pcenter = (rmax, rmax) - image = make_ellipse(pcenter, radii, row.rotation) - wcs = make_wcs(center=gcenter, size=image.shape, - pixelsize=self.resolution, - frame="Galactic", projection="CAR") - idx, val = self.sky.reproject_from(image, wcs, squeeze=True) - templates.append((idx, val)) - logger.info("Done simulate %d cluster templates" % len(templates)) - self.templates = templates - - def _simulate_single(self, data, frequency): - """ - Simulate one single cluster at the specified frequency, based - on the cached template image. - - Parameters - ---------- - data : namedtuple - The data of the SNR to be simulated, given in a ``namedtuple`` - object, from which can get the required properties by - ``data.key``. - e.g., elements of `self.catalog.itertuples()` - frequency : float - The simulation frequency (unit: `self.freq_unit`). - - Returns - ------- - idx : 1D `~numpy.ndarray` - The indexes of the effective map pixels for the single cluster. - val : 1D `~numpy.ndarray` - The values (i.e., brightness temperature) of each map pixel - with respect the above indexes. - - See Also - -------- - `self._simulate_template()` for more detailed description. - """ - idx, val = self.templates[data.Index] - # Calculate the brightness temperature - luminosity = data.luminosity - distance = data.distance - specindex = data.specindex - size = (data.size_major/60.0, data.size_minor/60.0) # [deg] - Tb = self._calc_Tb(luminosity, distance, specindex, frequency, size) - val = val * Tb - return (idx, val) - - def _make_filepath(self, **kwargs): - """ - Make the path of output file according to the filename pattern - and output directory loaded from configurations. - """ - data = { - "prefix": self.prefix, - } - data.update(kwargs) - filename = self.filename_pattern.format(**data) - filepath = os.path.join(self.output_dir, filename) - return filepath - - def _make_header(self): - """ - Make the header with detail information (e.g., parameters and - history) for the simulated products. - """ - header = fits.Header() - header["COMP"] = (self.name, "Emission component") - header["BUNIT"] = ("K", "data unit is Kelvin") - header["CREATOR"] = (__name__, "File creator") - # TODO: - history = [] - comments = [] - for hist in history: - header.add_history(hist) - for cmt in comments: - header.add_comment(cmt) - self.header = header - logger.info("Created FITS header") - - def output(self, skymap, frequency): - """ - Write the simulated free-free map to disk with proper header - keywords and history. - - Returns - ------- - outfile : str - The (absolute) path to the output sky map file. - """ - outfile = self._make_filepath(frequency=frequency) - if not hasattr(self, "header"): - self._make_header() - header = self.header.copy() - header["FREQ"] = (frequency, "Frequency [ MHz ]") - header["DATE"] = ( - datetime.now(timezone.utc).astimezone().isoformat(), - "File creation date" - ) - if self.use_float: - skymap = skymap.astype(np.float32) - sky = self.sky.copy() - sky.data = skymap - sky.header = header - sky.write(outfile, clobber=self.clobber, checksum=self.checksum) - return outfile - - def preprocess(self): - """ - Perform the preparation procedures for the final simulations. - - Attributes - ---------- - _preprocessed : bool - This attribute presents and is ``True`` after the preparation - procedures are performed, which indicates that it is ready to - do the final simulations. - """ - if hasattr(self, "_preprocessed") and self._preprocessed: - return - # - logger.info("{name}: preprocessing ...".format(name=self.name)) - self._load_catalog() - self._process_catalog() - # - self._scale_catalog_coverage() - self._add_random_position() - self._add_random_eccentricity() - self._add_random_rotation() - # - self._calc_sizes() - self._calc_luminosity() - self._calc_specindex() - # - self._simulate_templates() - # - self._preprocessed = True - - def simulate_frequency(self, frequency): - """ - Simulate the sky map of all extragalactic clusters of galaxies at - the specified frequency. - - Parameters - ---------- - frequency : float - The simulation frequency (unit: `self.freq_unit`). - - Returns - ------- - hpmap_f : 1D `~numpy.ndarray` - The sky map at the input frequency. - filepath : str - The (absolute) path to the output sky file if saved, - otherwise ``None``. - - See Also - -------- - `self._simulate_template()` for more detailed description. - """ - self.preprocess() - # - logger.info("Simulating {name} map at {freq} ({unit}) ...".format( - name=self.name, freq=frequency, unit=self.freq_unit)) - skymap_f = np.zeros(self.sky.shape) - # XXX/TODO: be parallel - for row in self.catalog.itertuples(): - # TODO: progress bar - index, value = self._simulate_single(row, frequency) - skymap_f[index] += value - # - if self.save: - filepath = self.output(skymap_f, frequency) - else: - filepath = None - return (skymap_f, filepath) - - def simulate(self, frequencies): - """ - Simulate the sky maps of all extragalactic clusters of galaxies - for every specified frequency. - - Parameters - ---------- - frequency : list[float] - List of frequencies (unit: `self.freq_unit`) where the - simulation performed. - - Returns - ------- - skymaps : list[1D `~numpy.ndarray`] - List of sky maps at each frequency. - paths : list[str] - List of (absolute) path to the output sky maps. - """ - skymaps = [] - paths = [] - for f in np.array(frequencies, ndmin=1): - skymap_f, outfile = self.simulate_frequency(f) - skymaps.append(skymap_f) - paths.append(outfile) - return (skymaps, paths) - - def postprocess(self): - """Perform the post-simulation operations before the end.""" - logger.info("{name}: postprocessing ...".format(name=self.name)) - # Save the catalog actually used in the simulation - self._save_catalog_inuse() |