aboutsummaryrefslogtreecommitdiffstats
path: root/fg21sim
diff options
context:
space:
mode:
Diffstat (limited to 'fg21sim')
-rw-r--r--fg21sim/configs/20-extragalactic.conf.spec2
-rw-r--r--fg21sim/extragalactic/__init__.py4
-rw-r--r--fg21sim/extragalactic/clusters/clusters.py732
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()