galactic: Add new component "galactic/snr": supernova remnants

NOTE: very preliminary, and missing many necessary functions!

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