diff options
| -rwxr-xr-x | bin/fg21sim | 67 | ||||
| -rw-r--r-- | fg21sim/galactic/__init__.py | 4 | ||||
| -rw-r--r-- | fg21sim/galactic/synchrotron.py | 34 | ||||
| -rwxr-xr-x | setup.py | 1 |
4 files changed, 102 insertions, 4 deletions
diff --git a/bin/fg21sim b/bin/fg21sim new file mode 100755 index 0000000..8098ac3 --- /dev/null +++ b/bin/fg21sim @@ -0,0 +1,67 @@ +#!/usr/bin/env python3 +# -*- mode: python -*- +# +# Copyright (c) 2016 Weitian LI <liweitianux@live.com> +# MIT license + +""" +Simulate the low-frequency radio foregrounds for the 21cm EoR signal. +""" + +import os +import sys +import argparse +import logging + +import numpy as np + +from fg21sim.configs import configs, validate_configs +from fg21sim.utils import setup_logging +from fg21sim.galactic import Synchrotron as GalacticSynchrotron + + +def main(): + parser = argparse.ArgumentParser( + description="Simulate the radio foregrounds for 21cm EoR signal") + parser.add_argument("config", help="user configuration file") + parser.add_argument("-l", "--log", dest="loglevel", default=None, + choices=["DEBUG", "INFO", "WARNING", + "ERROR", "CRITICAL"], + help="set the log level") + parser.add_argument("-L", "--logfile", default=None, + help="filename where to save the log messages") + parser.add_argument("-Q", "--quiet", action="store_true", + help="be quiet so do not log messages to screen") + args = parser.parse_args() + + configs.read_userconfig(args.config) + validate_configs(configs) + + if args.quiet: + log_stream = "" + else: + log_stream = None + + setup_logging(dict_config=configs.logging, + level=args.loglevel, + stream=log_stream, + logfile=args.logfile) + tool = os.path.basename(sys.argv[0]) + logger = logging.getLogger(tool) + logger.info("COMMAND: {0}".format(" ".join(sys.argv))) + + frequencies = np.array(configs.frequencies, ndmin=1) + freq_unit = configs.getn("frequency/unit") + logger.info("Simulation frequencies: " + "{min:.2f} - {max:.2f} {unit} (#{num:d})".format( + min=frequencies.min(), max=frequencies.max(), + num=len(frequencies), unit=freq_unit)) + + logger.info("Simulating the Galactic synchrotron foreground ...") + gsynchrotron = GalacticSynchrotron(configs) + map_gsync = gsynchrotron.simulate(frequencies) + logger.info("Done simulate Galactic synchrotron foreground!") + + +if __name__ == "__main__": + main() diff --git a/fg21sim/galactic/__init__.py b/fg21sim/galactic/__init__.py new file mode 100644 index 0000000..4fe2659 --- /dev/null +++ b/fg21sim/galactic/__init__.py @@ -0,0 +1,4 @@ +# Copyright (c) 2016 Weitian LI <liweitianux@live.com> +# MIT license + +from .synchrotron import Synchrotron diff --git a/fg21sim/galactic/synchrotron.py b/fg21sim/galactic/synchrotron.py index d05292f..614aae9 100644 --- a/fg21sim/galactic/synchrotron.py +++ b/fg21sim/galactic/synchrotron.py @@ -6,6 +6,7 @@ Diffuse Galactic synchrotron emission (unpolarized) simulations. """ import os +import logging from datetime import datetime, timezone import numpy as np @@ -14,6 +15,9 @@ import astropy.units as au import healpy as hp +logger = logging.getLogger(__name__) + + class Synchrotron: """Simulate the diffuse Galactic synchrotron emission based on an existing template. @@ -40,20 +44,21 @@ class Synchrotron: def _set_configs(self): """Load the configs and set the corresponding class attributes.""" - self.template_path = self.configs.getn( + self.template_path = self.configs.get_path( "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 = self.configs.getn( + self.indexmap_path = self.configs.get_path( "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.output_dir = self.configs.get_path( + "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") @@ -63,6 +68,8 @@ class Synchrotron: self.lmax = self.configs.getn("common/lmax") # unit of the frequency self.freq_unit = au.Unit(self.configs.getn("frequency/unit")) + # + logger.info("Loaded and setup configurations") def _load_template(self): """Load the template map""" @@ -71,6 +78,8 @@ class Synchrotron: self.template_header = fits.header.Header(header) self.template_nside = self.template_header["NSIDE"] self.template_ordering = self.template_header["ORDERING"] + logger.info("Loaded template map from {0} (Nside={1})".format( + self.template_path, self.template_nside)) def _load_indexmap(self): """Load the spectral index map""" @@ -79,6 +88,8 @@ class Synchrotron: self.indexmap_header = fits.header.Header(header) self.indexmap_nside = self.indexmap_header["NSIDE"] self.indexmap_ordering = self.indexmap_header["ORDERING"] + logger.info("Loaded spectral index map from {0} (Nside={1})".format( + self.indexmap_path, self.indexmap_nside)) def _process_template(self): """Upgrade/downgrade the template to have the same Nside as @@ -88,6 +99,8 @@ class Synchrotron: else: # upgrade/downgrade the resolution self.hpmap = hp.ud_grade(self.template, nside_out=self.nside) + logger.info("Upgrade/downgrade template map from Nside " + "{0} to {1}".format(self.template_nside, self.nside)) def _process_indexmap(self): """Upgrade/downgrade the spectral index map to have the same @@ -97,6 +110,8 @@ class Synchrotron: else: # upgrade/downgrade the resolution self.hpmap_index = hp.ud_grade(self.indexmap, nside_out=self.nside) + logger.info("Upgrade/downgrade spectral index map from Nside " + "{0} to {1}".format(self.indexmap_nside, self.nside)) def _add_smallscales(self): """Add fluctuations on small scales to the template map. @@ -128,11 +143,12 @@ class Synchrotron: 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) + gss = hp.synfast(cls=cl, nside=self.nside, new=True) # whiten the Gaussian random field gss = (gss - gss.mean()) / gss.std() self.hpmap_smallscales = alpha * gss * self.hpmap**beta self.hpmap += self.hpmap_smallscales + logger.info("Added small-scale fluctuations") def _transform_frequency(self, frequency): """Transform the template map to the requested frequency, @@ -149,6 +165,7 @@ class Synchrotron: header = fits.Header() header["COMP"] = ("Galactic synchrotron (unpolarized)", "Emission component") + header["CREATOR"] = (__name__, "File creator") # TODO: history = [] comments = [] @@ -157,6 +174,7 @@ class Synchrotron: for cmt in comments: header.add_comment(cmt) self.header = header + logger.info("Created FITS header") def output(self, hpmap, frequency): """Write the simulated synchrotron map to disk with proper @@ -167,8 +185,13 @@ class Synchrotron: np.dtype("float64"): "D", } # + 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() @@ -185,6 +208,7 @@ class Synchrotron: format=FITS_COLUMN_FORMATS.get(hpmap.dtype)) ], header=header) hdu.writeto(filepath, clobber=self.clobber, checksum=True) + logger.info("Write simulated map to file: {0}".format(filepath)) def simulate(self, frequencies): """Simulate the synchrotron map at the specified frequencies.""" @@ -196,6 +220,8 @@ class Synchrotron: # hpmaps = [] for f in np.array(frequencies, ndmin=1): + logger.info("Simulating synchrotron map at {0} ({1}) ...".format( + f, self.freq_unit)) hpmap_f = self._transform_frequency(f) hpmaps.append(hpmap_f) if self.save: @@ -45,6 +45,7 @@ setup( ], packages=find_packages(exclude=["docs", "tests"]), scripts=[ + "bin/fg21sim", "bin/healpix2hpx", "bin/hpx2healpix", ], |