path: root/fg21sim/galactic/freefree.py

# Copyright (c) 2016-2017 Weitian LI <weitian@aaronly.me>
# MIT license

"""
Diffuse Galactic free-free emission simulations.

References
----------
.. [dickinson2003]
   Dickinson, C.; Davies, R. D.; Davis, R. J.,
   "Towards a free-free template for CMB foregrounds",
   2003, MNRAS, 341, 369,
   http://adsabs.harvard.edu/abs/2003MNRAS.341..369D

.. [finkbeiner2003]
   Finkbeiner, Douglas P.,
   "A Full-Sky Hα Template for Microwave Foreground Prediction",
   2003, ApJS, 146, 407,
   http://adsabs.harvard.edu/abs/2003ApJS..146..407F

.. [schlegel1998]
   Schlegel, David J.; Finkbeiner, Douglas P.; Davis, Marc,
   "Maps of Dust Infrared Emission for Use in Estimation of Reddening
   and Cosmic Microwave Background Radiation Foregrounds",
   1998, ApJ, 500, 525,
   http://adsabs.harvard.edu/abs/1998ApJ...500..525S
"""

import os
import logging

import numpy as np

from ..sky import get_sky


logger = logging.getLogger(__name__)


class FreeFree:
    """
    Simulate the diffuse Galactic free-free emission.

    The [dickinson2003] method is followed to derive the free-free template.
    The all-sky Hα survey map [Finkbeiner2003] is first corrected for dust
    absorption using the infrared 100-μm dust map [Schlegel1998],
    and then converted to free-free emission map (brightness temperature).

    Parameters
    ----------
    configs : `~ConfigManager`
        An ``ConfigManager`` object contains default and user configurations.
        For more details, see the example config specification.

    Attributes
    ----------
    TODO
    """
    # Component name
    compID = "galactic/freefree"
    name = "Galactic free-free"

    def __init__(self, configs):
        self.configs = configs
        self._set_configs()

        self.sky = get_sky(configs)
        self.sky.add_header("CompID", self.compID, "Emission component ID")
        self.sky.add_header("CompName", self.name, "Emission component")
        self.sky.add_header("BUNIT", "K", "[Kelvin] Data unit")
        self.sky.creator = __name__

    def _set_configs(self):
        """
        Load the configs and set the corresponding class attributes.
        """
        comp = self.compID
        self.halphamap_path = self.configs.get_path(comp+"/halphamap")
        self.dustmap_path = self.configs.get_path(comp+"/dustmap")
        self.f_dust = self.configs.getn(comp+"/dust_fraction")
        self.halpha_abs_th = self.configs.getn(comp+"/halpha_abs_th")  # [mag]
        self.Te = self.configs.getn(comp+"/electron_temperature")  # [K]
        self.prefix = self.configs.getn(comp+"/prefix")
        self.output_dir = self.configs.get_path(comp+"/output_dir")
        #
        self.filename_pattern = self.configs.getn("output/filename_pattern")
        self.clobber = self.configs.getn("output/clobber")
        self.frequencies = self.configs.frequencies  # [MHz]
        #
        logger.info("Loaded and set up configurations")

    def _load_maps(self):
        """
        Load the Hα map, and 100-μm dust map.
        """
        logger.info("Loading H[alpha] map ...")
        self.halphamap = self.sky.open(self.halphamap_path)
        logger.info("Loading dust map ...")
        self.dustmap = self.sky.open(self.dustmap_path)

    def _correct_dust_absorption(self):
        """
        Correct the Hα map for dust absorption using the
        100-μm dust map.

        References: Ref.[dickinson2003],Eq.(1,3),Sec.(2.5)
        """
        if hasattr(self, "_dust_corrected") and self._dust_corrected:
            return

        logger.info("Correct H[alpha] map for dust absorption")
        logger.info("Effective dust fraction: {0}".format(self.f_dust))
        # Mask the regions where the true Halpha absorption is uncertain.
        # When the dust absorption goes rather large, the true Halpha
        # absorption can not well determined.
        # Corresponding dust absorption threshold, unit: [ MJy / sr ]
        dust_abs_th = self.halpha_abs_th / 0.0462 / self.f_dust
        logger.info("Dust absorption mask threshold: " +
                    "{0:.1f} MJy/sr ".format(dust_abs_th) +
                    "<-> H[alpha] absorption threshold: " +
                    "{0:.1f} mag".format(self.halpha_abs_th))
        mask = (self.dustmap.data > dust_abs_th)
        self.dustmap.data[mask] = np.nan
        fp_mask = 100 * mask.sum() / self.dustmap.data.size
        logger.warning("Dust map masked fraction: {0:.1f}%".format(fp_mask))
        #
        halphamap_corr = (self.halphamap.data *
                          10**(self.dustmap.data * 0.0185 * self.f_dust))
        self.halphamap.data = halphamap_corr
        self._dust_corrected = True
        logger.info("Done dust absorption correction")

    def _calc_factor_a(self, nu):
        """
        Calculate the ratio factor a(Te, ν), which will be used to
        convert the Halpha emission [Rayleigh] to free-free emission
        brightness temperature [K].

        Parameters
        ----------
        nu : float
            The frequency where to calculate the factor a(nu).
            Unit: [MHz]

        Returns
        -------
        a : float
            The factor for Hα to free-free conversion.

        References: [dickinson2003],Eq.(8)
        """
        term1 = 0.183 * nu**0.1 * self.Te**(-0.15)
        term2 = 3.91 - np.log(nu) + 1.5*np.log(self.Te)
        a = term1 * term2
        return a

    def _calc_halpha_to_freefree(self, nu):
        """
        Calculate the conversion factor between Hα emission [Rayleigh]
        to radio free-free emission [K] at frequency ν [MHz].

        Parameters
        ----------
        nu : float
            The frequency where to calculate the conversion factor.
            Unit: [MHz]

        Returns
        -------
        h2f : float
            The conversion factor between Hα emission and free-free emission.

        References: [dickinson2003],Eq.(11)
        NOTE: The above referred formula has a superfluous "10^3" term!
        """
        a = self._calc_factor_a(nu)
        h2f = 38.86 * a * nu**(-2.1) * 10**(290/self.Te) * self.Te**0.667
        return h2f

    def _outfilepath(self, frequency, **kwargs):
        """
        Generate the path/filename to the output file for writing
        the simulate sky images.

        Parameters
        ----------
        frequency : float
            The frequency of the output sky image.
            Unit: [MHz]

        Returns
        -------
        filepath : str
            The generated filepath for the output sky file.
        """
        filename = self.filename_pattern.format(
            prefix=self.prefix, frequency=frequency, **kwargs)
        filepath = os.path.join(self.output_dir, filename)
        return 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 subsequent simulations.
        """
        if hasattr(self, "_preprocessed") and self._preprocessed:
            return

        logger.info("{name}: preprocessing ...".format(name=self.name))
        self._load_maps()
        self._correct_dust_absorption()

        self._preprocessed = True

    def simulate_frequency(self, frequency):
        """
        Simulate the free-free map at the specified frequency.

        Parameters
        ----------
        frequency : float
            The frequency where to simulate the emission map.
            Unit: [MHz]

        Returns
        -------
        sky : `~SkyBase`
            The simulated sky image as a new sky instance.
        """
        logger.info("Simulating {name} map at {freq:.2f} [MHz] ...".format(
            name=self.name, freq=frequency))
        sky = self.sky.copy()
        sky.frequency = frequency
        ratio_K_R = self._calc_halpha_to_freefree(frequency)
        sky.data = self.halphamap * ratio_K_R
        logger.info("Done simulate map at %.2f [MHz]." % frequency)
        return sky

    def simulate(self, frequencies=None):
        """
        Simulate the emission maps.

        Parameters
        ----------
        frequencies : float, or list[float]
            The frequencies where to simulate the emission map.
            Unit: [MHz]
            Default: None (i.e., use ``self.frequencies``)

        Returns
        -------
        skyfiles : list[str]
            List of the filepath to the written sky files
        """
        if frequencies is None:
            frequencies = self.frequencies
        else:
            frequencies = np.array(frequencies, ndmin=1)

        logger.info("Simulating {name} ...".format(name=self.name))
        skyfiles = []
        for freq in frequencies:
            sky = self.simulate_frequency(freq)
            outfile = self._outfilepath(frequency=freq)
            sky.write(outfile)
            skyfiles.append(outfile)
        logger.info("Done simulate {name}!".format(name=self.name))
        return skyfiles

    def postprocess(self):
        """Perform the post-simulation operations before the end."""
        logger.info("{name}: postprocessing ...".format(name=self.name))
        logger.info("^_^ nothing to do :-)")