aboutsummaryrefslogtreecommitdiffstats
path: root/fg21sim/extragalactic/clusters.py
diff options
context:
space:
mode:
Diffstat (limited to 'fg21sim/extragalactic/clusters.py')
-rw-r--r--fg21sim/extragalactic/clusters.py650
1 files changed, 650 insertions, 0 deletions
diff --git a/fg21sim/extragalactic/clusters.py b/fg21sim/extragalactic/clusters.py
new file mode 100644
index 0000000..b20d261
--- /dev/null
+++ b/fg21sim/extragalactic/clusters.py
@@ -0,0 +1,650 @@
+# Copyright (c) 2016 Weitian LI <liweitianux@live.com>
+# MIT license
+
+"""
+Simulation of the radio emissions from clusters of galaxies.
+
+NOTE
+----
+Currently, only radio *halos* are considered with many simplifications.
+Radio *relics* simulations need more investigations ...
+"""
+
+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 healpy as hp
+import pandas as pd
+
+from ..utils import write_fits_healpix
+from ..utils.random import spherical_uniform
+from ..utils.convert import Fnu_to_Tb
+from ..utils.grid import make_grid_ellipse, map_grid_to_healpix
+
+
+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 = "clusters of galaxies"
+
+ def __init__(self, configs):
+ self.configs = 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") * au.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.clobber = self.configs.getn("output/clobber")
+ self.nside = self.configs.getn("common/nside")
+ 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, so do NOT save.")
+ 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:
+ 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_particles"].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 _expand_catalog_fullsky(self):
+ """Expand the catalog to be full sky, by assuming that clusters are
+ uniformly distributed. Also, the radio halo fraction is also
+ considered to determine the final number of clusters on the full sky.
+ """
+ fullsky = 4*np.pi * au.sr
+ factor = float(fullsky / self.catalog_prop["coverage"])
+ n0_cluster = len(self.catalog)
+ logger.info("Halo fraction in clusters: %f" % self.halo_fraction)
+ # Total number of clusters on the full sky
+ N_cluster = int(n0_cluster * factor * self.halo_fraction)
+ logger.info("Total number of clusters on the full sky: {0:,}".format(
+ N_cluster))
+ logger.info("Expanding the catalog to be full sky ...")
+ idx = np.round(np.random.uniform(low=0, high=n0_cluster-1,
+ size=N_cluster)).astype(np.int)
+ self.catalog = self.catalog.iloc[idx, :]
+ self.catalog.reset_index(inplace=True)
+ logger.info("Done expand the catalog to be full 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
+ self.units["coord"] = au.deg
+ 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: degree) of the clusters calculated
+ from the above virial radii and the random eccentricities.
+
+ NOTE
+ ----
+ The elliptical major and minor axes are calculated by assuming
+ the equal area between the ellipse and corresponding circle.
+ theta2 = r_vir / distance
+ pi * a * b = pi * (theta2)^2
+ e = sqrt((a^2 - b^2) / a^2)
+ thus,
+ a = theta2 / (1-e^2)^(1/4)
+ b = theta2 * (1-e^2)^(1/4)
+ """
+ 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"]
+ theta2 = float(2 * r_vir / distance) # [ rad ]
+ size_major = theta2 / (1 - self.catalog["eccentricity"]**2) ** 0.25
+ size_minor = theta2 * (1 - self.catalog["eccentricity"]**2) ** 0.25
+ self.catalog["size_major"] = size_major * au.rad.to(au.deg)
+ self.catalog["size_minor"] = size_minor * au.rad.to(au.deg)
+ self.units["size"] = au.deg
+ 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.
+
+ 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
+ # XXX: why the mass conversion ??
+ 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 * float(mass_RB / (1e15*au.solMass))**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"] = 1.4 * au.GHz
+ self.units["luminosity"] = au.W / au.Hz
+ logger.info("Done Calculate the radio luminosity")
+
+ 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: `self.units["luminosity"]`) at
+ the reference frequency (i.e.,
+ `self.catalog_prop["luminosity_freq"]`).
+ distance : float
+ The luminosity distance (unit: `self.units["distance"]`) 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: `self.freq_unit`) where the brightness
+ temperature requested.
+ size : 2-float tuple
+ The (major, minor) axes (unit: `self.units["size"]`).
+ 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.
+ """
+ freq = frequency * self.freq_unit
+ freq_ref = self.catalog_prop["luminosity_freq"]
+ luminosity = luminosity * self.units["luminosity"]
+ Lnu = luminosity * float(freq / freq_ref) ** (-specindex)
+ Fnu = Lnu / (4*np.pi * (distance*self.units["distance"])**2)
+ omega = size[0]*self.units["size"] * size[1]*self.units["size"]
+ Tb = Fnu_to_Tb(Fnu, omega, freq)
+ return Tb.value
+
+ def _simulate_templates(self):
+ """Simulate the template (HEALPix) 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.
+ Each element is a `(hpidx, hpval)` tuple with `hpidx` the
+ indexes of effective HEALPix pixels (RING ordering) and `hpval`
+ the values of the corresponding pixels.
+ e.g., ``[ (hpidx1, hpval1), (hpidx2, hpval2), ... ]``
+ """
+ logger.info("Simulating HEALPix templates for each cluster ...")
+ templates = []
+ resolution = self.resolution.to(au.deg).value
+ # XXX/TODO: be parallel
+ for row in self.catalog.itertuples():
+ # TODO: progress bar
+ center = (row.glon, row.glat)
+ size = (row.size_major, row.size_minor) # already [ degree ]
+ rotation = row.rotation # already [ degree ]
+ grid = make_grid_ellipse(center, size, resolution, rotation)
+ hpidx, hpval = map_grid_to_healpix(grid, self.nside)
+ templates.append((hpidx, hpval))
+ 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
+ -------
+ hpidx : 1D `~numpy.ndarray`
+ The indexes (in RING ordering) of the effective HEALPix
+ pixels for the SNR.
+ hpval : 1D `~numpy.ndarray`
+ The values (i.e., brightness temperature) of each HEALPix
+ pixel with respect the above indexes.
+
+ See Also
+ --------
+ `self._simulate_template()` for more detailed description.
+ """
+ name = data.name
+ hpidx, hpval = self.templates[name]
+ # Calculate the brightness temperature
+ luminosity = data.luminosity
+ distance = data.distance
+ specindex = data.specindex
+ size = (data.size_major, data.size_minor)
+ Tb = self._calc_Tb(luminosity, distance, specindex, frequency, size)
+ hpval = hpval * Tb
+ return (hpidx, hpval)
+
+ 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)
+ filetype = self.configs.getn("output/filetype")
+ if filetype == "fits":
+ filename += ".fits"
+ else:
+ raise NotImplementedError("unsupported filetype: %s" % filetype)
+ 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"] = ("extragalactic clusters of galaxies (radio halos)",
+ "Emission component")
+ header["UNIT"] = ("Kelvin", "Map unit")
+ 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, hpmap, frequency):
+ """Write the simulated free-free map to disk with proper header
+ keywords and history.
+ """
+ if not os.path.exists(self.output_dir):
+ os.mkdir(self.output_dir)
+ logger.info("Created output dir: {0}".format(self.output_dir))
+ #
+ filepath = 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:
+ hpmap = hpmap.astype(np.float32)
+ write_fits_healpix(filepath, hpmap, header=header,
+ clobber=self.clobber)
+ logger.info("Write simulated map to file: {0}".format(filepath))
+
+ 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._expand_catalog_fullsky()
+ self._add_random_position()
+ self._add_random_eccentricity()
+ self._add_random_rotation()
+ #
+ self._calc_sizes()
+ self._calc_luminosity()
+ #
+ self._simulate_templates()
+ #
+ self._preprocessed = True
+
+ def simulate_frequency(self, frequency):
+ """Simulate the emission (HEALPix) map of all Galactic SNRs at
+ the specified frequency.
+
+ Parameters
+ ----------
+ frequency : float
+ The simulation frequency (unit: `self.freq_unit`).
+
+ Returns
+ -------
+ hpmap_f : 1D `~numpy.ndarray`
+ HEALPix map data in RING ordering
+
+ 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))
+ hpmap_f = np.zeros(hp.nside2npix(self.nside))
+ for row in self.catalog.itertuples():
+ hpidx, hpval = self._simulate_single(row, frequency)
+ hpmap_f[hpidx] += hpval
+ #
+ if self.save:
+ self.output(hpmap_f, frequency)
+ return hpmap_f
+
+ def simulate(self, frequencies):
+ """Simulate the emission (HEALPix) maps of all Galactic SNRs for
+ every specified frequency.
+
+ Parameters
+ ----------
+ frequency : list[float]
+ List of frequencies (unit: `self.freq_unit`) where the
+ simulation performed.
+
+ Returns
+ -------
+ hpmaps : list[1D `~numpy.ndarray`]
+ List of HEALPix maps (in RING ordering) at each frequency.
+ """
+ hpmaps = []
+ for f in np.array(frequencies, ndmin=1):
+ hpmap_f = self.simulate_frequency(f)
+ hpmaps.append(hpmap_f)
+ return hpmaps
+
+ 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()