diff options
| -rwxr-xr-x | bin/get-healpix-patch | 10 | ||||
| -rw-r--r-- | fg21sim/configs/00-general.conf.spec | 2 | ||||
| -rw-r--r-- | fg21sim/configs/10-galactic.conf.spec | 7 | ||||
| -rw-r--r-- | fg21sim/configs/20-extragalactic.conf.spec | 11 | ||||
| -rw-r--r-- | fg21sim/extragalactic/clusters/emission.py | 4 | ||||
| -rw-r--r-- | fg21sim/galactic/snr.py | 12 | ||||
| -rw-r--r-- | fg21sim/sky.py | 27 | ||||
| -rw-r--r-- | fg21sim/utils/units.py | 7 | ||||
| -rw-r--r-- | fg21sim/utils/wcs.py | 14 |
9 files changed, 62 insertions, 32 deletions
diff --git a/bin/get-healpix-patch b/bin/get-healpix-patch index 7f01f78..0dc27f6 100755 --- a/bin/get-healpix-patch +++ b/bin/get-healpix-patch @@ -22,6 +22,7 @@ from fg21sim.configs import configs from fg21sim.sky import SkyPatch from fg21sim.utils import setup_logging from fg21sim.utils.fits import read_fits_healpix +from fg21sim.utils.units import UnitConversions as AUC def main(): @@ -39,7 +40,7 @@ def main(): help="size of the sky patch; " + "format: xsize,ysize; unit: pixel") parser.add_argument("--pixelsize", dest="pixelsize", type=float, - help="pixel size of the sky patch; unit: arcmin") + help="pixel size of the sky patch; unit: [arcsec]") parser.add_argument("-S", "--smooth", dest="smooth", action="store_true", help="Smooth the output patch using a Gaussian " + "filter whose sigma is 2x the pixel size " + @@ -59,7 +60,7 @@ def main(): configs.getn("sky/patch/ycenter")) # [ deg ] size = (configs.getn("sky/patch/xsize"), configs.getn("sky/patch/ysize")) - pixelsize = configs.getn("sky/patch/pixelsize") # [ arcmin ] + pixelsize = configs.getn("sky/patch/pixelsize") # [ arcsec ] elif not all([args.center, args.size, args.pixelsize]): raise ValueError("--center, --size, and --pixelsize are " + "required when --config is missing!") @@ -75,7 +76,7 @@ def main(): logger.info("patch center: (%.3f, %.3f) [deg]" % center) logger.info("patch size: (%d, %d) pixels" % size) - logger.info("patch pixel size: %.1f [arcmin]" % pixelsize) + logger.info("patch pixel size: %.1f [arcsec]" % pixelsize) sky = SkyPatch(size=size, pixelsize=pixelsize, center=center) logger.info("Read HEALPix map from file: %s" % args.infile) @@ -86,7 +87,8 @@ def main(): input_data=(hpdata, hpheader["COORDSYS"]), output_projection=sky.wcs, shape_out=size) if args.smooth: - sigma = 2 * hp.nside2resol(nside, arcmin=True) / pixelsize + sigma = (2.0 * hp.nside2resol(nside, arcmin=True) * + AUC.arcmin2arcsec / pixelsize) image = scipy.ndimage.gaussian_filter(image, sigma=sigma) logger.info("Smoothed sky patch using Gaussian filter of " + "sigma = %.2f [pixel]" % sigma) diff --git a/fg21sim/configs/00-general.conf.spec b/fg21sim/configs/00-general.conf.spec index b5bfd3f..906be3e 100644 --- a/fg21sim/configs/00-general.conf.spec +++ b/fg21sim/configs/00-general.conf.spec @@ -49,7 +49,7 @@ type = option("patch", "healpix", default="patch") xsize = integer(default=None, min=1) ysize = integer(default=None, min=1) - # Pixel size [ arcmin ] + # Pixel size [ arcsec ] pixelsize = float(default=None, min=0.0) # Configurations for input HEALPix sky diff --git a/fg21sim/configs/10-galactic.conf.spec b/fg21sim/configs/10-galactic.conf.spec index c01aed7..1ef15eb 100644 --- a/fg21sim/configs/10-galactic.conf.spec +++ b/fg21sim/configs/10-galactic.conf.spec @@ -76,9 +76,10 @@ # Output the effective/inuse SNRs catalog data (CSV file) catalog_outfile = string(default=None) - # Resolution (unit: arcmin) for simulating each SNR, which are finally - # mapped to the HEALPix map of Nside specified in "[common]" section. - resolution = float(default=0.5, min=0.0) + # Resolution for simulating each SNR template, which are finally + # mapped to the all-sky HEALPix map if used. + # Unit: [arcsec] + resolution = float(default=30.0, min=5.0) # Filename prefix for this component prefix = string(default="gsnr") diff --git a/fg21sim/configs/20-extragalactic.conf.spec b/fg21sim/configs/20-extragalactic.conf.spec index 7678a6a..9216663 100644 --- a/fg21sim/configs/20-extragalactic.conf.spec +++ b/fg21sim/configs/20-extragalactic.conf.spec @@ -31,9 +31,10 @@ # The fraction that a cluster hosts a radio halo halo_fraction = float(default=None, min=0.0, max=1.0) - # Resolution (unit: arcmin) for simulating each cluster, which are finally - # mapped to the HEALPix map of Nside specified in "[common]" section. - resolution = float(default=0.5, min=0.0) + # Resolution for simulating each cluster templates, which are finally + # mapped to the all-sky HEALPix map if used. + # Unit: [arcsec] + resolution = float(default=30.0, min=5.0) # Filename prefix for this component prefix = string(default="egcluster") @@ -59,6 +60,10 @@ # event rather than accretion. (unit: Msun) merger_mass_min = float(default=1e12, min=1e10) + # Mass ratio of the main and sub clusters, below which is regarded as + # a major merger event. + ratio_major = float(default=3.0, min=1.0, max=10.0) + # Radius of the giant radio halo in clusters (unit: kpc) # XXX: currently only support a constant radius of halos radius = float(default=500.0, min=100.0) diff --git a/fg21sim/extragalactic/clusters/emission.py b/fg21sim/extragalactic/clusters/emission.py index 345e3db..93e04ff 100644 --- a/fg21sim/extragalactic/clusters/emission.py +++ b/fg21sim/extragalactic/clusters/emission.py @@ -184,7 +184,7 @@ class SynchrotronEmission: ---------- pixelsize : float The pixel size of the output simulated sky image - Unit: [arcmin] + Unit: [arcsec] Returns ------- @@ -195,7 +195,7 @@ class SynchrotronEmission: DA = self.cosmo.DL(self.z) * AUC.Mpc2cm # [cm] radius = self.radius * AUC.kpc2cm # [cm] omega = (np.pi * radius**2 / DA**2) * AUC.rad2deg**2 # [deg^2] - pixelarea = (pixelsize/60.0) ** 2 # [deg^2] + pixelarea = (pixelsize * AUC.arcsec2deg) ** 2 # [deg^2] if omega < pixelarea: omega = pixelarea F_nu = self.flux(nu) diff --git a/fg21sim/galactic/snr.py b/fg21sim/galactic/snr.py index 6897ce4..c44fc33 100644 --- a/fg21sim/galactic/snr.py +++ b/fg21sim/galactic/snr.py @@ -19,6 +19,7 @@ from ..sky import get_sky from ..utils.wcs import make_wcs from ..utils.convert import Fnu_to_Tb_fast from ..utils.grid import make_ellipse +from ..utils.units import UnitConversions as AUC logger = logging.getLogger(__name__) @@ -80,7 +81,7 @@ class SuperNovaRemnants: 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") # [ arcmin ] + self.resolution = self.configs.getn(comp+"/resolution") # [ arcsec ] self.prefix = self.configs.getn(comp+"/prefix") self.save = self.configs.getn(comp+"/save") self.output_dir = self.configs.get_path(comp+"/output_dir") @@ -100,7 +101,8 @@ class SuperNovaRemnants: * glon, glat : SNR coordinate, Galactic coordinate, [deg] * size_major, size_minor : SNR angular sizes; major and minor axes of the ellipse fitted to the SNR, or the diameter of the fitted - circle if the SNR is nearly circular; [arcmin] + circle if the SNR is nearly circular; + originally in [arcmin], converted to [arcsec] * flux : Flux density at 1 GHz, [Jy] """ self.catalog = pd.read_csv(self.catalog_path) @@ -111,6 +113,10 @@ class SuperNovaRemnants: nrow, ncol)) # The flux densities are given at 1 GHz self.catalog_flux_freq = (1.0*au.GHz).to(self.freq_unit).value + # Convert ``size_major`` and ``size_minor`` from unit [arcmin] + # to [arcsec] + self.catalog["size_major"] *= AUC.arcmin2arcsec + self.catalog["size_minor"] *= AUC.arcmin2arcsec def _save_catalog_inuse(self): """ @@ -239,7 +245,7 @@ class SuperNovaRemnants: freq_ref = self.catalog_flux_freq # [ MHz ] Fnu = flux * (frequency / freq_ref) ** (-specindex) # [ Jy ] omega = size[0] * size[1] # [ deg^2 ] - pixelarea = (self.sky.pixelsize/60.0) ** 2 # [ deg^2 ] + pixelarea = (self.sky.pixelsize * AUC.arcsec2deg) ** 2 # [ deg^2 ] if omega < pixelarea: # The object is smaller than a pixel, so round up to a pixel area omega = pixelarea diff --git a/fg21sim/sky.py b/fg21sim/sky.py index 637c990..1460ec2 100644 --- a/fg21sim/sky.py +++ b/fg21sim/sky.py @@ -21,6 +21,7 @@ import healpy as hp from .utils.wcs import make_wcs from .utils.fits import read_fits_healpix, write_fits_healpix +from .utils.units import UnitConversions as AUC logger = logging.getLogger(__name__) @@ -45,29 +46,25 @@ class SkyPatch: has shape (ysize, xsize). pixelsize : float The pixel size of the sky patch, will be used to determine - the sky coordinates. [ arcmin ] + the sky coordinates. + Unit: [arcsec] center : (ra, dec) tuple - The (R.A., Dec.) coordinate of the sky patch center. [ deg ] + The (R.A., Dec.) coordinate of the sky patch center. + Unit: [deg] infile : str The path to the input sky patch Attributes ---------- - type : str, "patch" or "healpix" + type_ : str, "patch" or "healpix" The type of this sky map data : 1D `numpy.ndarray` The flattened 1D array of map data NOTE: The 2D image is flattened to 1D, making it easier to be manipulated in a similar way as the HEALPix map. - size : int tuple, (width, height) - The dimensions of the FITS image shape : int tuple, (nrow*ncol, ) The shape of the flattened image array NOTE: nrow=height, ncol=width - pixelsize : float - The pixel size of the sky map [ arcmin ] - center : float tuple, (ra, dec) - The (R.A., Dec.) coordinate of the sky patch center. [ deg ] """ type_ = "patch" # Input sky patch and its frequency [ MHz ] @@ -93,9 +90,10 @@ class SkyPatch: XXX/FIXME --------- + Assumed a flat sky! Consider the spherical coordination and WCS sky projection!! """ - ps = self.pixelsize / 60.0 # [ deg ] + ps = self.pixelsize * AUC.arcsec2deg # [deg] size = self.xsize * self.ysize return size * ps**2 @@ -175,7 +173,7 @@ class SkyPatch: header = self.header.copy(strip=True) wcs_header = self.wcs.to_header() header.extend(wcs_header, update=True) - header["PIXSIZE"] = (self.pixelsize, "Pixel size [arcmin]") + header["PIXSIZE"] = (self.pixelsize, "Pixel size [arcsec]") header["RA0"] = (self.center[0], "R.A. of patch center [deg]") header["DEC0"] = (self.center[1], "Dec. of patch center [deg]") hdu = fits.PrimaryHDU(data=image, header=header) @@ -284,7 +282,8 @@ class SkyHealpix: shape : int tuple, (npix,) The shape (i.e., length) of the HEALPix array pixelsize : float - The pixel size of the HEALPix map [ arcmin ] + The pixel size of the HEALPix map + Unit: [arcsec] """ type_ = "healpix" # Input sky patch and its frequency [ MHz ] @@ -317,7 +316,9 @@ class SkyHealpix: @property def pixelsize(self): - return hp.nside2resol(self.nside, arcmin=True) + ps = hp.nside2resol(self.nside, arcmin=True) + ps *= AUC.arcmin2arcsec + return ps def read(self, infile, frequency=None): """ diff --git a/fg21sim/utils/units.py b/fg21sim/utils/units.py index aea0f92..bb8f670 100644 --- a/fg21sim/utils/units.py +++ b/fg21sim/utils/units.py @@ -50,6 +50,13 @@ class UnitConversions: keV2erg = au.keV.to(au.erg) # Angle rad2deg = au.rad.to(au.deg) + deg2rad = au.deg.to(au.rad) + arcmin2deg = au.arcmin.to(au.deg) + deg2arcmin = au.deg.to(au.arcmin) + arcsec2deg = au.arcsec.to(au.deg) + deg2arcsec = au.deg.to(au.arcsec) + arcmin2arcsec = au.arcmin.to(au.arcsec) + arcsec2arcmin = au.arcsec.to(au.arcmin) class Constants: diff --git a/fg21sim/utils/wcs.py b/fg21sim/utils/wcs.py index b4a9ed1..3cb76b9 100644 --- a/fg21sim/utils/wcs.py +++ b/fg21sim/utils/wcs.py @@ -8,6 +8,8 @@ Create WCS for sky projection. import numpy as np from astropy.wcs import WCS +from .units import UnitConversions as AUC + def make_wcs(center, size, pixelsize, frame="ICRS", projection="TAN"): """ @@ -16,11 +18,13 @@ def make_wcs(center, size, pixelsize, frame="ICRS", projection="TAN"): Parameters ---------- center : (xcenter, ycenter) float tuple - The equatorial/galactic coordinate of the sky/image center [ deg ]. + The equatorial/galactic coordinate of the sky/image center. + Unit: [deg] size : (xsize, ysize) int tuple The size (width, height) of the sky/image. pixelsize : float - The pixel size of the sky/image [ arcmin ] + The pixel size of the sky/image. + Unit: [arcsec] frame : str, "ICRS" or "Galactic" The coordinate frame, only one of ``ICRS`` or ``Galactic``. projection : str, "TAN" or "CAR" @@ -31,10 +35,14 @@ def make_wcs(center, size, pixelsize, frame="ICRS", projection="TAN"): ------- w : `~astropy.wcs.WCS` Created WCS header/object + + TODO/XXX + --------- + To support other projections, e.g., ``SIN`` (common in radio astronomy). """ xcenter, ycenter = center # [ deg ] xsize, ysize = size - delt = pixelsize / 60.0 # [ deg ] + delt = pixelsize * AUC.arcsec2deg # [ deg ] if projection.upper() not in ["TAN", "CAR"]: raise ValueError("unsupported projection: " % projection) if frame.upper() == "ICRS": |