aboutsummaryrefslogtreecommitdiffstats
diff options
context:
space:
mode:
-rw-r--r--fg21sim/galactic/synchrotron.py204
1 files changed, 204 insertions, 0 deletions
diff --git a/fg21sim/galactic/synchrotron.py b/fg21sim/galactic/synchrotron.py
new file mode 100644
index 0000000..14b4065
--- /dev/null
+++ b/fg21sim/galactic/synchrotron.py
@@ -0,0 +1,204 @@
+# Copyright (c) 2016 Weitian LI <liweitianux@live.com>
+# MIT license
+
+"""
+Diffuse Galactic synchrotron emission (unpolarized) simulations.
+"""
+
+import os
+from datetime import datetime, timezone
+
+import numpy as np
+from astropy.io import fits
+import astropy.units as au
+import healpy as hp
+
+
+class Synchrotron:
+ """Simulate the diffuse Galactic synchrotron emission based on an
+ existing template.
+
+ Parameters
+ ----------
+ configs : ConfigManager object
+ An `ConfigManager` object contains default and user configurations.
+ For more details, see the example config specification.
+
+ Attributes
+ ----------
+ ???
+
+ References
+ ----------
+ ???
+ """
+ def __init__(self, configs):
+ self.configs = configs
+ self._set_configs()
+ self._load_template()
+ self._load_indexmap()
+
+ def _set_configs(self):
+ """Load the configs and set the corresponding class attributes."""
+ data_dir = self.configs.getn("common/data_dir")
+ self.template_path = os.path.join(
+ data_dir, self.configs.getn("galactic/synchrotron/template"))
+ self.template_freq = self.configs.getn(
+ "galactic/synchrotron/template_freq")
+ self.template_unit = au.Unit(
+ self.configs.getn("galactic/synchrotron/template_unit"))
+ self.indexmap_path = os.path.join(
+ data_dir, self.configs.getn("galactic/synchrotron/indexmap"))
+ self.smallscales = self.configs.getn(
+ "galactic/synchrotron/add_smallscales")
+ # output
+ self.prefix = self.configs.getn("galactic/synchrotron/prefix")
+ self.save = self.configs.getn("galactic/synchrotron/save")
+ self.output_dir = self.configs.getn("galactic/synchrotron/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")
+ # common
+ self.nside = self.configs.getn("common/nside")
+ self.lmin = self.configs.getn("common/lmin")
+ self.lmax = self.configs.getn("common/lmax")
+ # unit of the frequency
+ self.freq_unit = au.Unit(self.configs.getn("frequency/unit"))
+
+ def _load_template(self):
+ """Load the template map"""
+ self.template, header = hp.read_map(self.template_path,
+ nest=False, h=True, verbose=False)
+ self.template_header = fits.header.Header(header)
+ self.template_nside = self.template_header["NSIDE"]
+ self.template_ordering = self.template_header["ORDERING"]
+
+ def _load_indexmap(self):
+ """Load the spectral index map"""
+ self.indexmap, header = hp.read_map(self.indexmap_path,
+ nest=False, h=True, verbose=False)
+ self.indexmap_header = fits.header.Header(header)
+ self.indexmap_nside = self.indexmap_header["NSIDE"]
+ self.indexmap_ordering = self.indexmap_header["ORDERING"]
+
+ def _process_template(self):
+ """Upgrade/downgrade the template to have the same Nside as
+ requested by the output."""
+ if self.nside == self.template_nside:
+ self.hpmap = self.template.copy()
+ else:
+ # upgrade/downgrade the resolution
+ self.hpmap = hp.ud_grade(self.template, nside_out=self.nside)
+
+ def _process_indexmap(self):
+ """Upgrade/downgrade the spectral index map to have the same
+ Nside as requested by the output."""
+ if self.nside == self.indexmap_nside:
+ self.hpmap_index = self.indexmap.copy()
+ else:
+ # upgrade/downgrade the resolution
+ self.hpmap_index = hp.ud_grade(self.indexmap, nside_out=self.nside)
+
+ def _add_smallscales(self):
+ """Add fluctuations on small scales to the template map.
+
+ XXX/TODO:
+ * Support using different models.
+ * This should be extensible/plug-able, e.g., a separate module
+ and allow easily add new models for use.
+
+ References
+ ----------
+ [1] M. Remazeilles et al. 2015, MNRAS, 451, 4311-4327
+ "An improved source-subtracted and destriped 408-MHz all-sky map"
+ Sec. 4.2: Small-scale fluctuations
+ """
+ if not self.smallscales:
+ return
+ # To add small scale fluctuations
+ # model: Remazeilles15
+ gamma = -2.703 # index of the power spectrum between l [30, 90]
+ sigma_temp = 56 # original beam resolution of the template [ arcmin ]
+ alpha = 0.0599
+ beta = 0.782
+ # angular power spectrum of the Gaussian random field
+ ell = np.arange(self.lmax+1).astype(np.int)
+ cl = np.zeros(ell.shape)
+ ell_idx = ell >= self.lmin
+ cl[ell_idx] = (ell[ell_idx] ** gamma *
+ 1.0 - np.exp(-ell[ell_idx]**2 * sigma_temp**2))
+ cl[ell < self.lmin] = cl[self.lmin]
+ # generate a realization of the Gaussian random field
+ gss = hp.synfast(cls=cl, nside=self.nside)
+ # whiten the Gaussian random field
+ gss = (gss - gss.mean()) / gss.std()
+ self.hpmap_smallscales = alpha * gss * self.hpmap**beta
+ self.hpmap += self.hpmap_smallscales
+
+ def _transform_frequency(self, frequency):
+ """Transform the template map to the requested frequency,
+ according to the spectral model and using an spectral index map.
+ """
+ hpmap_f = (self.hpmap *
+ (frequency / self.template_freq) ** self.hpmap_index)
+ 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 synchrotron (unpolarized)",
+ "Emission component")
+ # TODO:
+ history = []
+ comments = []
+ for hist in history:
+ header.add_history(hist)
+ for cmt in comments:
+ header.add_comment(cmt)
+ self.header = header
+
+ def output(self, hpmap, frequency):
+ """Write the simulated synchrotron map to disk with proper
+ header keywords and history.
+ """
+ FITS_COLUMN_FORMATS = {
+ np.dtype("float32"): "E",
+ np.dtype("float64"): "D",
+ }
+ #
+ filename = self.filename_pattern.format(prefix=self.prefix,
+ frequency=frequency)
+ 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)
+ hdu = fits.BinTableHDU.from_columns([
+ fits.Column(name="I", array=hpmap,
+ format=FITS_COLUMN_FORMATS.get(hpmap.dtype))
+ ], header=header)
+ hdu.writeto(filepath, clobber=self.clobber, checksum=True)
+
+ def simulate(self, frequencies):
+ """Simulate the synchrotron map at the specified frequencies."""
+ if not hasattr(self, "hpmap"):
+ self._process_template()
+ self._add_smallscales()
+ if not hasattr(self, "hpmap_index"):
+ self._process_indexmap()
+ #
+ hpmaps = []
+ for f in np.array(frequencies, ndmin=1):
+ hpmap_f = self._transform_frequency(f)
+ hpmaps.append(hpmap_f)
+ if self.save:
+ self.output(hpmap_f, f)
+ return hpmaps