From e252342e4fb0c106734e5659203d36985242bc7c Mon Sep 17 00:00:00 2001
From: Aaron LI <aaronly.me@outlook.com>
Date: Tue, 11 Oct 2016 10:00:55 +0800
Subject: galactic: Add free-free component simulation
* New class "galactic.FreeFree" to simulate free-free emission
* Add new config section "galactic/freefree"
NOTE: current untested
---
fg21sim/galactic/__init__.py | 1 +
fg21sim/galactic/freefree.py | 227 +++++++++++++++++++++++++++++++++++++++++++
2 files changed, 228 insertions(+)
create mode 100644 fg21sim/galactic/freefree.py
(limited to 'fg21sim/galactic')
diff --git a/fg21sim/galactic/__init__.py b/fg21sim/galactic/__init__.py
index 4fe2659..ebf7999 100644
--- a/fg21sim/galactic/__init__.py
+++ b/fg21sim/galactic/__init__.py
@@ -2,3 +2,4 @@
# MIT license
from .synchrotron import Synchrotron
+from .freefree import FreeFree
diff --git a/fg21sim/galactic/freefree.py b/fg21sim/galactic/freefree.py
new file mode 100644
index 0000000..9a51777
--- /dev/null
+++ b/fg21sim/galactic/freefree.py
@@ -0,0 +1,227 @@
+# Copyright (c) 2016 Weitian LI <liweitianux@live.com>
+# MIT license
+
+"""
+Diffuse Galactic free-free 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
+
+from ..utils import read_fits_healpix, write_fits_healpix
+
+
+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 H\alpha survey map [Finkbeiner2003] is first corrected for dust
+ absorption using the infrared 100-\mu{}m dust map [Schlegel1998],
+ and then converted to free-free emission map (brightness temperature).
+
+ Parameters
+ ----------
+ configs : ConfigManager object
+ An `ConfigManager` object contains default and user configurations.
+ For more details, see the example config specification.
+
+ Attributes
+ ----------
+ ???
+
+ 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
+ """
+ def __init__(self, configs):
+ self.configs = configs
+ self._set_configs()
+ self._load_halphamap()
+ self._load_dustmap()
+
+ def _set_configs(self):
+ """Load the configs and set the corresponding class attributes."""
+ self.halphamap_path = self.configs.get_path(
+ "galactic/freefree/halphamap")
+ self.halphamap_unit = au.Unit(
+ self.configs.getn("galactic/freefree/halphamap_unit"))
+ self.dustmap_path = self.configs.get_path(
+ "galactic/freefree/dustmap")
+ self.dustmap_unit = au.Unit(
+ self.configs.getn("galactic/freefree/dustmap_unit"))
+ # output
+ self.prefix = self.configs.getn("galactic/freefree/prefix")
+ self.save = self.configs.getn("galactic/freefree/save")
+ self.output_dir = self.configs.get_path(
+ "galactic/freefree/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")
+ logger.info("Loaded and set up configurations")
+
+ def _load_halphamap(self):
+ """Load the H{\alpha} map, and upgrade/downgrade the resolution
+ to match the output Nside."""
+ self.halphamap, self.halphamap_header = read_fits_healpix(
+ self.halphamap_path)
+ halphamap_nside = self.halphamap_header["NSIDE"]
+ logger.info("Loaded H[alpha] map from {0} (Nside={1})".format(
+ self.halphamap_path, halphamap_nside))
+ # TODO: Validate & convert unit
+ if self.halphamap_unit != au.Unit("Rayleigh"):
+ raise ValueError("unsupported Halpha map unit: {0}".format(
+ self.halphamap_unit))
+ # Upgrade/downgrade resolution
+ if halphamap_nside != self.nside:
+ self.halphamap = hp.ud_grade(self.halphamap, nside_out=self.nside)
+ logger.info("Upgrade/downgrade H[alpha] map from Nside "
+ "{0} to {1}".format(halphamap_nside, self.nside))
+
+ def _load_dustmap(self):
+ """Load the dust map, and upgrade/downgrade the resolution
+ to match the output Nside."""
+ self.dustmap, self.dustmap_header = read_fits_healpix(
+ self.dustmap_path)
+ dustmap_nside = self.dustmap_header["NSIDE"]
+ logger.info("Loaded dust map from {0} (Nside={1})".format(
+ self.dustmap_path, dustmap_nside))
+ # TODO: Validate & convert unit
+ if self.dustmap_unit != au.Unit("MJy / sr"):
+ raise ValueError("unsupported dust map unit: {0}".format(
+ self.dustmap_unit))
+ # Upgrade/downgrade resolution
+ if dustmap_nside != self.nside:
+ self.dustmap = hp.ud_grade(self.dustmap, nside_out=self.nside)
+ logger.info("Upgrade/downgrade dust map from Nside "
+ "{0} to {1}".format(dustmap_nside, self.nside))
+
+ def _correct_dust_absorption(self):
+ """Correct the H{\alpha} map for dust absorption using the
+ 100-{\mu}m dust map.
+
+ References: [Dickinson2003], Eq.(3)
+ """
+ if hasattr(self, "dust_corrected") and self.dust_corrected:
+ return
+ #
+ # Effective dust fraction in the LoS actually absorbing Halpha
+ f_dust = 0.33
+ halphamap_corr = self.halphamap * 10**(self.dustmap * 0.0185 * f_dust)
+ self.halphamap = halphamap_corr
+ self.dust_corrected = True
+
+ def _calc_ratio_a(self, Te, nu_GHz):
+ """Calculate the ratio factor a(T, nu), which will be used to
+ transform the Halpha emission (Rayleigh) to free-free emission
+ brightness temperature (mK).
+
+ References: [Dickinson2003], Eq.(8)
+ """
+ term1 = 0.366 * nu_GHz**0.1 * Te**(-0.15)
+ term2 = np.log(4.995e-2 / nu_GHz) + 1.5*np.log(Te)
+ a = term1 * term2
+ return a
+
+ def _simulate_frequency(self, frequency):
+ """Simulate the free-free map at the specified frequency.
+
+ References: [Dickinson2003], Eq.(11)
+
+ NOTE: [Dickinson2003], Eq.(11) may wrongly have the "10^3" term.
+ """
+ # Assumed electron temperature [ K ]
+ Te = 7000.0
+ T4 = Te / 1e4
+ nu = frequency * self.freq_unit.to(au.GHz) # frequency [ GHz ]
+ ratio_a = self._calc_ratio_a(T4, nu)
+ # NOTE: ignored the "10^3" term in the referred equation
+ ratio_mK_R = (8.396 * ratio_a * nu**(-2.1) *
+ T4**0.667 * 10**(0.029/T4) * (1+0.08))
+ # Use "Kelvin" as the brightness temperature unit
+ ratio_K_R = ratio_mK_R * au.mK.to(au.K)
+ hpmap_f = self.halphamap * ratio_K_R
+ return hpmap_f
+
+ 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 free-free emission",
+ "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))
+ #
+ filename = self.filename_pattern.format(prefix=self.prefix,
+ frequency=frequency)
+ filename += ".fits"
+ filepath = os.path.join(self.output_dir, filename)
+ 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."""
+ self._correct_dust_absorption()
+ #
+ 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
--
cgit v1.2.2