diff options
| -rw-r--r-- | fg21sim/galactic/__init__.py | 1 | ||||
| -rw-r--r-- | fg21sim/galactic/snr.py | 255 |
2 files changed, 256 insertions, 0 deletions
diff --git a/fg21sim/galactic/__init__.py b/fg21sim/galactic/__init__.py index ebf7999..5bdd659 100644 --- a/fg21sim/galactic/__init__.py +++ b/fg21sim/galactic/__init__.py @@ -3,3 +3,4 @@ from .synchrotron import Synchrotron from .freefree import FreeFree +from .snr import SuperNovaRemnants diff --git a/fg21sim/galactic/snr.py b/fg21sim/galactic/snr.py new file mode 100644 index 0000000..7709a4f --- /dev/null +++ b/fg21sim/galactic/snr.py @@ -0,0 +1,255 @@ +# Copyright (c) 2016 Weitian LI <liweitianux@live.com> +# MIT license + +""" +Galactic supernova remnants (SNRs) emission simulations. +""" + +import os +import logging +from datetime import datetime, timezone + +import numpy as np +from astropy.io import fits +import astropy.units as au +import healpy as hp +import pandas as pd + +from ..utils import write_fits_healpix + + +logger = logging.getLogger(__name__) + + +class SuperNovaRemnants: + """ + Simulate the Galactic supernova remnants emission. + + ??? + + Parameters + ---------- + configs : ConfigManager object + An `ConfigManager` object contains default and user configurations. + For more details, see the example config specification. + + Attributes + ---------- + ??? + + References + ---------- + .. [Green2014] + Green, D. A., + "A catalogue of 294 Galactic supernova remnants", + 2014, Bulletin of the Astronomical Society of India, 42, 47-58, + http://adsabs.harvard.edu/abs/2014BASI...42...47G + + .. [GreenSNRDataWeb] + A Catalogue of Galactic Supernova Remnants + http://www.mrao.cam.ac.uk/surveys/snrs/ + """ + def __init__(self, configs): + self.configs = configs + self._set_configs() + + def _set_configs(self): + """Load the configs and set the corresponding class attributes.""" + comp = "galactic/snr" + self.catalog_path = self.configs.get_path(comp+"/catalog") + self.catalog_outfile = self.configs.get_path(comp+"/catalog_outfile") + 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")) + # + logger.info("Loaded and set up configurations") + + def _load_catalog(self): + """Load the Galactic SNRs catalog data.""" + self.catalog = pd.read_csv(self.catalog_path) + nrow, ncol = self.catalog.shape + logger.info("Loaded SNRs catalog data from: {0}".format( + self.catalog_path)) + logger.info("SNRs catalog data: {0} objects, {1} columns".format( + nrow, ncol)) + # Set the units for columns + self.units = { + "glon": au.deg, + "glat": au.deg, + "size": au.arcmin, + "flux": au.Jy, + } + # The flux densities are given at 1 GHz + self.catalog_flux_freq = 1.0 * au.GHz + + def _save_catalog_inuse(self): + """Save the effective/inuse SNRs catalog data to a CSV file. + + NOTE + ---- + - Only the effective/inuse SNRs are saved (i.e., without the ones + that are filtered out). + - Also save the simulated rotation column. + - The unnecessary columns are striped. + """ + if self.catalog_outfile is None: + logger.warning("Catalog output file not set, so do NOT save.") + return + # TODO/XXX + pass + + def _filter_catalog(self): + """Filter the catalog data to remove the objects with incomplete + data. + + The following cases are filtered out: + - Missing angular size + - Missing flux density data + - Missing spectral index value + + NOTE + ---- + The objects with uncertain data are currently kept. + """ + cond1 = pd.isnull(self.catalog["size_major"]) + cond2 = pd.isnull(self.catalog["size_minor"]) + cond3 = pd.isnull(self.catalog["flux"]) + cond4 = pd.isnull(self.catalog["specindex"]) + cond_keep = ~(cond1 | cond2 | cond3 | cond4) + n_total = len(cond_keep) + n_delete = cond_keep.sum() + n_delete_p = n_delete / n_total * 100 + n_remain = n_total - n_delete + self.catalog = self.catalog[cond_keep] + # Reset index + self.catalog.reset_index(inplace=True) + logger.info("SNRs catalog: filtered out " + + "{0:d} ({1:.1f}) objects".format(n_delete, n_delete_p)) + logger.info("SNRs catalog: remaining {0} objects".format(n_remain)) + + def _add_random_rotation(self): + """Add random rotation angles for each SNR as column "rotation" + within the `catalog` data frame. + + 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. + """ + num = len(self.catalog) + angles = np.random.uniform(low=0.0, high=360.0, size=num) + rotation = pd.Series(data=angles, name="rotation") + self.catalog["rotation"] = rotation + logger.info("Added random rotation angles as the 'rotation' column") + + def _calc_Tb(self, flux, specindex, frequency): + """Calculate the brightness temperature at requested frequency + by assuming a power-law spectral shape. + + Parameters + ---------- + flux : float + The flux density (unit: [ Jy ]) at the reference frequency + (`self.catalog_flux_freq`). + specindex : float + The spectral index of the power-law spectrum + frequency : float + The frequency (unit: [ MHz ]) where the brightness temperature + requested. + + 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. + """ + pass + + def _simulate_frequency(self, frequency): + """Simulate the Galactic SNRs emission map at the specified frequency. + """ + pass + + def _make_filename(self, **kwargs): + """Make the path of output file according to the filename pattern + and output directory loaded from configurations. + """ + data = { + "prefix": self.prefix, + } + data.extend(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"] = ("Galactic supernova remnants (SNRs)", + "Emission component") + 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 simulate(self, frequencies): + """Simulate the free-free map at the specified frequencies.""" + hpmaps = [] + for f in np.array(frequencies, ndmin=1): + logger.info("Simulating free-free map at {0} ({1}) ...".format( + f, self.freq_unit)) + hpmap_f = self._simulate_frequency(f) + hpmaps.append(hpmap_f) + if self.save: + self.output(hpmap_f, f) + return hpmaps |