From 2f40f984a44c393b72a62b36cf054171dabf169d Mon Sep 17 00:00:00 2001 From: Jason Ma Date: Thu, 27 Oct 2016 12:13:07 +0800 Subject: extragalactic/pointsource (#3) Merge PR#3: Add new simulation component "extragalactic/pointsources", including SF, SB, RQ, FRI and FRII. * extragalactic/pointsource: add point source simulation module * extracgalactic/pointsource: add point source simulation module * Add configurating spec to extragalgactic point sources. * Modified some variables * base.py: modified * flux.py: modified * fr1.py: modified * fr2.py: modified * pointsources.py: modified * psparams.py: modified * radioquiet.py: modified * starforming.py: modified * starbursting.py: modified * Rewritten the comments. * base.py: modified * flux.py: modified * fr1.py: modified * fr2.py: modified * Modified * psparams.py: modified * radioquiet.py: modified * starforming.py: modified * starbursting.py: modified * Modified * Modified * Modified * Modified * Modified * MOdified * Modified * Modified * Modified * Modified * Modified * Modified * Modified * Modified * Modified * Changed pointsource to pointsources * Fixed some config keywords * Fixed some config keywords * Fixed some config keywords * Fixed some config keywords * Fixed some config keywords * Fixed some config keywords * Fixed some config keywords * base.py:rewrited * fr1.py:rewritten * pointsources.py: rewritten * radioquiet.py: rewritten * Rewritten * Rewritten * base.py: modified * fr1.py: modified * fr2.py: modified * radioquiet.py: modified * starbursting.py: modified * starforming.py: modified * Fix conflicts * fg21sim: fixed conflicts * base.py: modified frequencies loading in _get_configs() * Rewritten as forground.py * fg21sim: fixed conflicts * base.py: deteled loading for frequencies configurations. * fr1.py: modified * fr2.py: modified * pointsources.py: modified * radioquiet.py: modified * starbursting.py: modified * starforming.py: modified * Add new methods to calculate Tb. * Add new methods to calculate Tb. * Add new methods to calculate Tb. * Add new methods to calculate Tb. * Add new methods to calculate Tb. * Add new methods to calculate Tb. * Deleted useless comments. * Add pscomps to deal with multi-type PS problem. * Add a new key. * Fixed permission to 755. * Rejusted PS subsections. * Add methods to calcualte luminosity function and redshift distribution. * Rejusted generation of samples redshift and luminosity. * Fixed mistakes on FRII structure, added hotspots and offsets. * Reajusted generation of samples redshift and luminosity. * Readujsted generation of samples radii, redshifts and luminosity. * Readujsted generation of samples radii, redshifts and luminosity. * Fixed conflicts. * Fixed conficts. * Combined configurations of pointsources. * Removed the older extragalactic configuration file. * Fixed some mistakes. * Fixed mistakes of drawing PS. * Fixed mistakes of drawing PS. * Fixed code style by pep8 checking. * Fixed code style by pep8 checking. * Fixed code style by pep8 checking. * Fixed some coding style. * Reconfigured default redshift interval. * Reconfigured default redshift interval. * Reconfigured default redshift interval. * Fixed mistakes in method calc_single_Tb and changed resolution of grid. * Fixed mistakes in method calc_single_Tb and changed resolution of grid. * Deleted astropy.units style code to accelerate. * Deleted astropy.units style code to accelerate. * Deleted astropy.units style code to accelerate. * Deleted astropy.units style code to accelerate. * Deleted astropy.units style code to accelerate. * Deleted astropy.units style code to accelerate. * Deleted astropy.units style code to accelerate. * Fixed some mistakes. * Fixed some mistakes. * Changed dA from au.Mpc to float64. * Fixed some mistakes. * Fixed some mistakes. * Fixed some mistakes. * Reajusted grid resolution to generate discs. * Reajusted grid resolution to generate discs. * Reajusted loading strategy of parameter resolution. * Reajusted code style of configuration loading. * Reajusted code style of configuration loading. * Reajusted code style of configuration loading. * Reajusted code style of configuration loading. * Reajusted code style of configuration loading. --- fg21sim/configs/20-extragalactic.conf.spec | 42 +++ fg21sim/extragalactic/__init__.py | 5 + fg21sim/extragalactic/pointsources/__init__.py | 4 + fg21sim/extragalactic/pointsources/base.py | 203 +++++++++++ fg21sim/extragalactic/pointsources/fr1.py | 355 +++++++++++++++++++ fg21sim/extragalactic/pointsources/fr2.py | 380 +++++++++++++++++++++ fg21sim/extragalactic/pointsources/pointsources.py | 126 +++++++ fg21sim/extragalactic/pointsources/psparams.py | 77 +++++ fg21sim/extragalactic/pointsources/radioquiet.py | 189 ++++++++++ fg21sim/extragalactic/pointsources/starbursting.py | 238 +++++++++++++ fg21sim/extragalactic/pointsources/starforming.py | 257 ++++++++++++++ fg21sim/foregrounds.py | 5 + 12 files changed, 1881 insertions(+) create mode 100644 fg21sim/extragalactic/pointsources/__init__.py create mode 100644 fg21sim/extragalactic/pointsources/base.py create mode 100644 fg21sim/extragalactic/pointsources/fr1.py create mode 100644 fg21sim/extragalactic/pointsources/fr2.py create mode 100644 fg21sim/extragalactic/pointsources/pointsources.py create mode 100644 fg21sim/extragalactic/pointsources/psparams.py create mode 100644 fg21sim/extragalactic/pointsources/radioquiet.py create mode 100644 fg21sim/extragalactic/pointsources/starbursting.py create mode 100644 fg21sim/extragalactic/pointsources/starforming.py (limited to 'fg21sim') diff --git a/fg21sim/configs/20-extragalactic.conf.spec b/fg21sim/configs/20-extragalactic.conf.spec index bef7d66..c6be967 100644 --- a/fg21sim/configs/20-extragalactic.conf.spec +++ b/fg21sim/configs/20-extragalactic.conf.spec @@ -39,3 +39,45 @@ save = boolean(default=True) # Output directory to save the simulated results output_dir = string(default=None) + + # Extragalactic point sources + [[pointsources]] + # Whether save this point source catelogue to disk + save = boolean(default=True) + # Output directory to save the simulated catelogues + output_dir = string(default="PS_tables") + # PS components to be simulated + pscomponents=string_list(default=list()) + # Resolution [arcmin] + resolution=float(default=0.6) + # Number of each type of point source + # Star forming + [[[starforming]]] + # Number of samples + numps = integer(default=1000) + # Prefix + prefix = string(default="SF") + + [[[starbursting]]] + # Number of samples + numps = integer(default=1000) + # Prefix + prefix = string(default="SB") + + [[[radioquiet]]] + # Number of samples + numps = integer(default=1000) + # Prefix + prefix = string(default="RQ") + + [[[FRI]]] + # Number of samples + numps = integer(default=1000) + # Prefix + prefix = string(default="FRI") + + [[[FRII]]] + # Number of samples + numps = integer(default=1000) + # Prefix + prefix = string(default="FRII") diff --git a/fg21sim/extragalactic/__init__.py b/fg21sim/extragalactic/__init__.py index af976ca..f57928f 100644 --- a/fg21sim/extragalactic/__init__.py +++ b/fg21sim/extragalactic/__init__.py @@ -2,3 +2,8 @@ # MIT license from .clusters import GalaxyClusters + +# Copyright (c) 2016 Zhixian MA +# MIT license + +from .pointsources import PointSources diff --git a/fg21sim/extragalactic/pointsources/__init__.py b/fg21sim/extragalactic/pointsources/__init__.py new file mode 100644 index 0000000..4c64c98 --- /dev/null +++ b/fg21sim/extragalactic/pointsources/__init__.py @@ -0,0 +1,4 @@ +# Copyright (c) 2016 Zhixian MA +# MIT license + +from .pointsources import PointSources diff --git a/fg21sim/extragalactic/pointsources/base.py b/fg21sim/extragalactic/pointsources/base.py new file mode 100644 index 0000000..1544f5b --- /dev/null +++ b/fg21sim/extragalactic/pointsources/base.py @@ -0,0 +1,203 @@ +# Copyright (c) 2016 Zhixian MA +# MIT license + +""" +Simulation of point sources for 21cm signal detection + +Point sources types +------------------- +1. Star forming (SF) galaxies +2. Star bursting (SB) galaxies +3. Radio quiet AGN (RQ_AGN) +4. Faranoff-Riley I (FRI) +5. Faranoff-Riley II (FRII) + +References +---------- +[1] Wilman et al., + "A semi-empirical simulation of the extragalactic radio continuum + sky for next generation radio telescopes", + 2008, MNRAS, 388, 1335-1348. + http://adsabs.harvard.edu/abs/2008MNRAS.388.1335W +[2] Jelic et al., + "Foreground simulations for the LOFAR-Epoch of Reionization + Experiment", + 2008, MNRAS, 389, 1319-1335. + http://adsabs.harvard.edu/abs/2008MNRAS.389.1319W +[3] Spherical uniform distribution + https://www.jasondavies.com/maps/random-points/ +""" +import os +import numpy as np +import pandas as pd +import healpy as hp + +from .psparams import PixelParams + + +class BasePointSource: + """ + The basic class of point sources + + Parameters + ---------- + z: float; + Redshift, z ~ U(0,20) + dA: au.Mpc; + Angular diameter distance, which is calculated according to the + cosmology constants. In this work, it is calculated by module + basic_params + lumo: au.Jy; + Luminosity at the reference frequency. + lat: au.deg; + The colatitude angle in the spherical coordinate system + lon: au.deg; + The longtitude angle in the spherical coordinate system + area: au.sr; + Area of the point sources, sr = rad^2 + + """ + # Init + + def __init__(self, configs): + # configures + self.configs = configs + # PS_list information + self.columns = ['z', 'dA (Mpc)', 'luminosity (Jy)', 'Lat (deg)', + 'Lon (deg)', 'Area (sr)'] + self.nCols = len(self.columns) + self._set_configs() + + def _set_configs(self): + """ + Load the configs and set the corresponding class attributes. + """ + comp = "extragalactic/pointsources/" + # common + self.nside = self.configs.getn("common/nside") + # resolution + self.resolution = self.configs.getn(comp+"resolution") + # save flag + self.save = self.configs.getn(comp+"save") + # Output_dir + self.output_dir = self.configs.get_path(comp+"output_dir") + + def calc_number_density(self): + pass + + def calc_cdf(self): + """ + Calculate cumulative distribution functions for simulating of + samples with corresponding reshift and luminosity. + + Parameter + ----------- + rho_mat: np.ndarray rho(lumo,z) + The number density matrix (joint-distribution of z and flux) + of this type of PS. + + Returns + ------- + cdf_z, cdf_lumo: np.ndarray + Cumulative distribution functions of redshift and flux. + """ + # Normalization + rho_mat = self.rho_mat + rho_sum = np.sum(rho_mat) + rho_norm = rho_mat / rho_sum + # probability distribution of redshift + pdf_z = np.sum(rho_norm, axis=0) + pdf_lumo = np.sum(rho_norm, axis=1) + # Cumulative function + cdf_z = np.zeros(pdf_z.shape) + cdf_lumo = np.zeros(pdf_lumo.shape) + for i in range(len(pdf_z)): + cdf_z[i] = np.sum(pdf_z[:i]) + for i in range(len(pdf_lumo)): + cdf_lumo[i] = np.sum(pdf_lumo[:i]) + + return cdf_z, cdf_lumo + + def get_lumo_redshift(self): + """ + Randomly generate redshif and luminosity at ref frequency using + the CDF functions. + + Paramaters + ---------- + df_z, cdf_lumo: np.ndarray + Cumulative distribution functions of redshift and flux. + zbin,lumobin: np.ndarray + Bins of redshif and luminosity. + + Returns + ------- + z: float + Redshift. + lumo: au.W/Hz/sr + Luminosity. + """ + # Uniformlly generate random number in interval [0,1] + rnd_z = np.random.uniform(0, 1) + rnd_lumo = np.random.uniform(0, 1) + # Get redshift + dist_z = np.abs(self.cdf_z - rnd_z) + idx_z = np.where(dist_z == dist_z.min()) + z = self.zbin[idx_z[0]] + # Get luminosity + dist_lumo = np.abs(self.cdf_lumo - rnd_lumo) + idx_lumo = np.where(dist_lumo == dist_lumo.min()) + lumo = 10 ** self.lumobin[idx_lumo[0]] + + return float(z), float(lumo) + + def gen_single_ps(self): + """ + Generate single point source, and return its data as a list. + """ + # Redshift and luminosity + self.z, self.lumo = self.get_lumo_redshift() + # angular diameter distance + self.param = PixelParams(self.z) + self.dA = self.param.dA + # W/Hz/Sr to Jy + dA = self.dA * 3.0856775814671917E+22 # Mpc to meter + self.lumo = self.lumo / dA**2 / (10.0**-24) # [Jy] + # Position + x = np.random.uniform(0, 1) + self.lat = (np.arccos(2 * x - 1) / np.pi * 180 - 90) # [deg] + self.lon = np.random.uniform(0, np.pi * 2) / np.pi * 180 # [deg] + # Area + npix = hp.nside2npix(self.nside) + self.area = 4 * np.pi / npix # [sr] + + ps_list = [self.z, self.dA, self.lumo, self.lat, self.lon, self.area] + return ps_list + + def gen_catalog(self): + """ + Generate num_ps of point sources and save them into a csv file. + """ + # Init + ps_table = np.zeros((self.num_ps, self.nCols)) + for x in range(self.num_ps): + ps_table[x, :] = self.gen_single_ps() + + # Transform into Dataframe + self.ps_catalog = pd.DataFrame(ps_table, columns=self.columns, + index=list(range(self.num_ps))) + + def save_as_csv(self): + """Save the catalog""" + if not os.path.exists(self.output_dir): + os.mkdir(self.output_dir) + + pattern = "{prefix}.csv" + filename = pattern.format(prefix=self.prefix) + + # save to csv + if self.save: + file_name = os.path.join(self.output_dir, filename) + self.ps_catalog.to_csv(file_name) + + return file_name diff --git a/fg21sim/extragalactic/pointsources/fr1.py b/fg21sim/extragalactic/pointsources/fr1.py new file mode 100644 index 0000000..0e052cd --- /dev/null +++ b/fg21sim/extragalactic/pointsources/fr1.py @@ -0,0 +1,355 @@ +# Copyright (c) 2016 Zhixian MA +# MIT license + +import numpy as np +import healpy as hp + +from .psparams import PixelParams +from .base import BasePointSource +from ...utils import grid +from ...utils import convert + + +class FRI(BasePointSource): + """ + Generate Faranoff-Riley I (FRI) AGN + + Parameters + ---------- + lobe_maj: float + The major half axis of the lobe + lobe_min: float + The minor half axis of the lobe + lobe_ang: float + The rotation angle of the lobe with respect to line of sight + + Reference + ---------- + [1] Wang J et al., + "How to Identify and Separate Bright Galaxy Clusters from the + Low-frequency Radio Sky?", + 2010, ApJ, 723, 620-633. + http://adsabs.harvard.edu/abs/2010ApJ...723..620W + [2] Fast cirles drawing + https://github.com/liweitianux/fg21sim/fg21sim/utils/draw.py + https://github.com/liweitianux/fg21sim/fg21sim/utils/grid.py + """ + + def __init__(self, configs): + super().__init__(configs) + self.columns.extend( + ['lobe_maj (rad)', 'lobe_min (rad)', 'lobe_ang (deg)']) + self.nCols = len(self.columns) + self._set_configs() + # Paramters for core/lobe ratio + # Willman et al. 2008 Sec2.5.(iii)-(iv) + self.xmed = -2.6 + # Lorentz factor of the jet + self.gamma = 6 + # Number density matrix + self.rho_mat = self.calc_number_density() + # Cumulative distribution of z and lumo + self.cdf_z, self.cdf_lumo = self.calc_cdf() + + def _set_configs(self): + """Load the configs and set the corresponding class attributes""" + super()._set_configs() + pscomp = "extragalactic/pointsources/FRI/" + # point sources amount + self.num_ps = self.configs.getn(pscomp+"numps") + # prefix + self.prefix = self.configs.getn(pscomp+"prefix") + # redshift bin + z_type = self.configs.getn(pscomp+"z_type") + if z_type == 'custom': + start = self.configs.getn(pscomp+"z_start") + stop = self.configs.getn(pscomp+"z_stop") + step = self.configs.getn(pscomp+"z_step") + self.zbin = np.arange(start, stop + step, step) + else: + self.zbin = np.arange(0.1, 10, 0.05) + # luminosity bin + lumo_type = self.configs.getn(pscomp+"lumo_type") + if lumo_type == 'custom': + start = self.configs.getn(pscomp+"lumo_start") + stop = self.configs.getn(pscomp+"lumo_stop") + step = self.configs.getn(pscomp+"lumo_step") + self.lumobin = np.arange(start, stop + step, step) + else: + self.lumobin = np.arange(20, 28, 0.1) # [W/Hz/sr] + + def calc_number_density(self): + """ + Calculate number density rho(lumo,z) of FRI + + References + ---------- + [1] Wilman et al., + "A semi-empirical simulation of the extragalactic radio continuum + sky for next generation radio telescopes", + 2008, MNRAS, 388, 1335-1348. + http://adsabs.harvard.edu/abs/2008MNRAS.388.1335W + [2] Willott et al., + "The radio luminosity function from the low-frequency 3CRR, + 6CE and 7CRS complete samples", + 2001, MNRAS, 322, 536-552. + http://adsabs.harvard.edu/abs/2001MNRAS.322..536W + + Returns + ------- + rho_mat: np.ndarray + Number density matris (joint-distribution of luminosity and + reshift). + """ + # Init + rho_mat = np.zeros((len(self.lumobin), len(self.zbin))) + # Parameters + # Refer to [2] Table. 1 model C and Willman's section 2.4 + alpha = 0.539 # spectral index + lumo_star = 10.0**26.1 # critical luminosity + rho_l0 = 10.0**(-7.120) # normalization constant + z1 = 0.706 # cut-off redshift + z2 = 2.5 # cut-off redshift adviced by Willman + k1 = 4.30 # index of space density revolution + # Calculation + for i, z in enumerate(self.zbin): + if z <= z1: + rho_mat[:, i] = ((rho_l0 * (10**self.lumobin / lumo_star) ** + (-alpha) * + np.exp(-10**self.lumobin / lumo_star)) * + (1 + z)**k1) + elif z <= z2: + rho_mat[:, i] = ((rho_l0 * (10**self.lumobin / lumo_star) ** + (-alpha) * + np.exp(-10**self.lumobin / lumo_star)) * + (1 + z1)**k1) + else: + rho_mat[:, i] = ((rho_l0 * (10**self.lumobin / lumo_star) ** + (-alpha) * + np.exp(-10**self.lumobin / lumo_star)) * + (1 + z)**-z2) + + return rho_mat + + def gen_lobe(self): + """ + Calculate lobe parameters + + References + ---------- + [1] Wilman et al., + "A semi-empirical simulation of the extragalactic radio continuum + sky for next generation radio telescopes", + 2008, MNRAS, 388, 1335-1348. + http://adsabs.harvard.edu/abs/2008MNRAS.388.1335W + + Return + ------ + lobe: list + lobe = [lobe_maj, lobe_min, lobe_ang], which represent the major + and minor axes and the rotation angle. + """ + D0 = 1 # [Mpc] + self.lobe_maj = 0.5 * np.random.uniform(0, D0 * (1 + self.z)**(-1.4)) + self.lobe_min = self.lobe_maj * np.random.uniform(0.2, 1) + self.lobe_ang = np.random.uniform(0, np.pi) / np.pi * 180 + + # Transform to pixel + self.lobe_maj = self.param.get_angle(self.lobe_maj) + self.lobe_min = self.param.get_angle(self.lobe_min) + lobe = [self.lobe_maj, self.lobe_min, self.lobe_ang] + + return lobe + + def gen_single_ps(self): + """ + Generate single point source, and return its data as a list. + """ + # Redshift and luminosity + self.z, self.lumo = self.get_lumo_redshift() + # angular diameter distance + self.param = PixelParams(self.z) + self.dA = self.param.dA + # W/Hz/Sr to Jy + dA = self.dA * 3.0856775814671917E+22 # Mpc to meter + self.lumo = self.lumo / dA**2 / (10.0**-24) # [Jy] + # Position + x = np.random.uniform(0, 1) + self.lat = (np.arccos(2 * x - 1) / np.pi * 180 - 90) # [deg] + self.lon = np.random.uniform(0, np.pi * 2) / np.pi * 180 # [deg] + # lobe + lobe = self.gen_lobe() + # Area + self.area = np.pi * self.lobe_maj * self.lobe_min + + ps_list = [self.z, self.dA, self.lumo, self.lat, self.lon, self.area] + ps_list.extend(lobe) + + return ps_list + + def draw_single_ps(self, freq): + """ + Designed to draw the elliptical lobes of FRI and FRII + + Prameters + --------- + nside: int and dyadic + self.ps_catalog: pandas.core.frame.DataFrame + Data of the point sources + ps_type: int + Class type of the point soruces + freq: float + frequency + """ + # Init + resolution = self.resolution / 60 # [degree] + npix = hp.nside2npix(self.nside) + hpmap = np.zeros((npix,)) + num_ps = self.ps_catalog.shape[0] + # Gen flux list + Tb_list = self.calc_Tb(freq) + ps_lobe = Tb_list[:, 1] + # Iteratively draw ps + for i in range(num_ps): + # Parameters + c_lat = self.ps_catalog['Lat (deg)'][i] # core lat [au.deg] + c_lon = self.ps_catalog['Lon (deg)'][i] # core lon [au.deg] + lobe_maj = self.ps_catalog['lobe_maj (rad)'][ + i] * 180 / np.pi # [deg] + lobe_min = self.ps_catalog['lobe_min (rad)'][ + i] * 180 / np.pi # [deg] + lobe_ang = self.ps_catalog['lobe_ang (deg)'][ + i] / 180 * np.pi # [rad] + + # Lobe1 + lobe1_lat = (lobe_maj / 2) * np.cos(lobe_ang) + lobe1_lat = c_lat + lobe1_lat + lobe1_lon = (lobe_maj / 2) * np.sin(lobe_ang) + lobe1_lon = c_lon + lobe1_lon + # draw + # Fill with ellipse + lon, lat, gridmap = grid.make_grid_ellipse( + (lobe1_lon, lobe1_lat), (lobe_maj, lobe_min), + resolution, lobe_ang / np.pi * 180) + indices, values = grid.map_grid_to_healpix( + (lon, lat, gridmap), self.nside) + hpmap[indices] += ps_lobe[i] + + # Lobe2 + lobe2_lat = (lobe_maj / 2) * np.cos(lobe_ang + np.pi) + lobe2_lat = c_lat + lobe2_lat + lobe2_lon = (lobe_maj / 2) * np.sin(lobe_ang + np.pi) + lobe2_lon = c_lon + lobe2_lon + # draw + # Fill with ellipse + lon, lat, gridmap = grid.make_grid_ellipse( + (lobe2_lon, lobe2_lat), (lobe_maj, lobe_min), + resolution, lobe_ang / np.pi * 180) + indices, values = grid.map_grid_to_healpix( + (lon, lat, gridmap), self.nside) + hpmap[indices] += ps_lobe[i] + + # Core + pix_tmp = hp.ang2pix(self.nside, + (self.ps_catalog['Lat (deg)'] + 90) / + 180 * np.pi, self.ps_catalog['Lon (deg)'] / + 180 * np.pi) + ps_core = Tb_list[:, 0] + hpmap[pix_tmp] += ps_core + + return hpmap + + def draw_ps(self, freq): + """ + Read csv ps list file, and generate the healpix structure vector + with the respect frequency. + """ + # Init + num_freq = len(freq) + npix = hp.nside2npix(self.nside) + hpmaps = np.zeros((npix, num_freq)) + + # Gen ps_catalog + self.gen_catalog() + # get hpmaps + for i in range(num_freq): + hpmaps[:, i] = self.draw_single_ps(freq[i]) + + return hpmaps + + def calc_single_Tb(self, area, freq): + """ + Calculate brightness temperatur of a single ps + + Parameters + ------------ + area: `~astropy.units.Quantity` + Area of the PS, e.g., `1.0*au.sr` + freq: `~astropy.units.Quantity` + Frequency, e.g., `1.0*au.MHz` + + Return + ------- + Tb:`~astropy.units.Quantity` + Average brightness temperature, e.g., `1.0*au.K` + """ + # Init + freq_ref = 151 # [MHz] + freq = freq # [MHz] + # Luminosity at 151MHz + lumo_151 = self.lumo # [Jy] + # Calc flux + # core-to-extend ratio + ang = self.lobe_ang / 180 * np.pi + x = np.random.normal(self.xmed, 0.5) + beta = np.sqrt((self.gamma**2 - 1) / self.gamma) + B_theta = 0.5 * ((1 - beta * np.cos(ang))**-2 + + (1 + beta * np.cos(ang))**-2) + ratio_obs = 10**x * B_theta + # Core + lumo_core = ratio_obs / (1 + ratio_obs) * lumo_151 + a0 = (np.log10(lumo_core) - 0.07 * + np.log10(freq_ref * 10.0E-3) + + 0.29 * np.log10(freq_ref * 10.0E-3) * + np.log10(freq_ref * 10.0E-3)) + lgs = (a0 + 0.07 * np.log10(freq * 10.0E-3) - 0.29 * + np.log10(freq * 10.0E-3) * + np.log10(freq * 10.0E-3)) + flux_core = 10**lgs # [Jy] + # core area + npix = hp.nside2npix(self.nside) + core_area = 4 * np.pi / npix # [sr] + Tb_core = convert.Fnu_to_Tb_fast(flux_core, core_area, freq) # [K] + # lobe + lumo_lobe = lumo_151 * (1 - ratio_obs) / (1 + ratio_obs) # [Jy] + flux_lobe = (freq / freq_ref)**(-0.75) * lumo_lobe + Tb_lobe = convert.Fnu_to_Tb_fast(flux_lobe, area, freq) # [K] + + Tb = [Tb_core, Tb_lobe] + return Tb + + def calc_Tb(self, freq): + """ + Calculate the surface brightness temperature of the point sources. + + Parameters + ------------ + area: `~astropy.units.Quantity` + Area of the PS, e.g., `1.0*au.sr` + freq: `~astropy.units.Quantity` + Frequency, e.g., `1.0*au.MHz` + + Return + ------ + Tb_list: list + Point sources brightness temperature + """ + # Tb_list + num_ps = self.ps_catalog.shape[0] + Tb_list = np.zeros((num_ps, 2)) + # Iteratively calculate Tb + for i in range(num_ps): + ps_area = self.ps_catalog['Area (sr)'][i] # [sr] + Tb_list[i, :] = self.calc_single_Tb(ps_area, freq) + + return Tb_list diff --git a/fg21sim/extragalactic/pointsources/fr2.py b/fg21sim/extragalactic/pointsources/fr2.py new file mode 100644 index 0000000..b40f9ff --- /dev/null +++ b/fg21sim/extragalactic/pointsources/fr2.py @@ -0,0 +1,380 @@ +# Copyright (c) 2016 Zhixian MA +# MIT license + +import numpy as np +import healpy as hp + +from .base import BasePointSource +from .psparams import PixelParams +from ...utils import grid +from ...utils import convert + + +class FRII(BasePointSource): + """ + Generate Faranoff-Riley II (FRII) AGN + + Parameters + ---------- + lobe_maj: float + The major half axis of the lobe + lobe_min: float + The minor half axis of the lobe + lobe_ang: float + The rotation angle of the lobe correspoind to line of sight + + Reference + ---------- + [1] Wang J et al., + "How to Identify and Separate Bright Galaxy Clusters from the + Low-frequency Radio Sky?", + 2010, ApJ, 723, 620-633. + http://adsabs.harvard.edu/abs/2010ApJ...723..620W + [2] Fast cirles drawing + https://github.com/liweitianux/fg21sim/fg21sim/utils/draw.py + https://github.com/liweitianux/fg21sim/fg21sim/utils/grid.py + """ + + def __init__(self, configs): + super().__init__(configs) + self.columns.extend( + ['lobe_maj (rad)', 'lobe_min (rad)', 'lobe_ang (deg)']) + self.nCols = len(self.columns) + self._set_configs() + # Paramters for core/lobe ratio + # Willman et al. 2008 Sec2.5.(iii)-(iv) + self.xmed = -2.8 + # Lorentz factor of the jet + self.gamma = 8 + # Number density matrix + self.rho_mat = self.calc_number_density() + # Cumulative distribution of z and lumo + self.cdf_z, self.cdf_lumo = self.calc_cdf() + + def _set_configs(self): + """Load the configs and set the corresponding class attributes""" + super()._set_configs() + pscomp = "extragalactic/pointsources/FRII/" + # point sources amount + self.num_ps = self.configs.getn(pscomp+"numps") + # prefix + self.prefix = self.configs.getn(pscomp+"prefix") + # redshift bin + z_type = self.configs.getn(pscomp+"z_type") + if z_type == 'custom': + start = self.configs.getn(pscomp+"z_start") + stop = self.configs.getn(pscomp+"z_stop") + step = self.configs.getn(pscomp+"z_step") + self.zbin = np.arange(start, stop + step, step) + else: + self.zbin = np.arange(0.1, 10, 0.05) + # luminosity bin + lumo_type = self.configs.getn(pscomp+"lumo_type") + if lumo_type == 'custom': + start = self.configs.getn(pscomp+"lumo_start") + stop = self.configs.getn(pscomp+"lumo_stop") + step = self.configs.getn(pscomp+"lumo_step") + self.lumobin = np.arange(start, stop + step, step) + else: + self.lumobin = np.arange(25.5, 30.5, 0.1) # [W/Hz/sr] + + def calc_number_density(self): + """ + Calculate number density rho(lumo,z) of FRI + + References + ---------- + [1] Wilman et al., + "A semi-empirical simulation of the extragalactic radio continuum + sky for next generation radio telescopes", + 2008, MNRAS, 388, 1335-1348. + http://adsabs.harvard.edu/abs/2008MNRAS.388.1335W + [2] Willott et al., + "The radio luminosity function from the low-frequency 3CRR, + 6CE and 7CRS complete samples", + 2001, MNRAS, 322, 536-552. + http://adsabs.harvard.edu/abs/2001MNRAS.322..536W + + Returns + ------- + rho_mat: np.ndarray + Number density matris (joint-distribution of luminosity and + reshift). + """ + # Init + rho_mat = np.zeros((len(self.lumobin), len(self.zbin))) + # Parameters + # Refer to [2] Table. 1 model C and Willman's section 2.4 + alpha = 2.27 # spectral index + lumo_star = 10.0**26.95 # critical luminosity + rho_l0 = 10.0**(-6.196) # normalization constant + z0 = 1.91 # center redshift + z2 = 1.378 # variance + # Calculation + for i, z in enumerate(self.zbin): + # space density revolusion + fh = np.exp(-0.5 * (z - z0)**2 / z2**2) + rho_mat[:, i] = ((rho_l0 * (10**self.lumobin / lumo_star) ** + (-alpha) * + np.exp(-lumo_star / 10.0**self.lumobin)) * + fh) + + return rho_mat + + def gen_lobe(self): + """ + Calculate lobe parameters + + References + ---------- + [1] Wilman et al., + "A semi-empirical simulation of the extragalactic radio continuum + sky for next generation radio telescopes", + 2008, MNRAS, 388, 1335-1348. + http://adsabs.harvard.edu/abs/2008MNRAS.388.1335W + + Return + ------ + lobe: list + lobe = [lobe_maj, lobe_min, lobe_ang], which represent the major + and minor axes and the rotation angle. + """ + D0 = 1 # [Mpc] + self.lobe_maj = 0.5 * np.random.uniform(0, D0 * (1 + self.z)**(-1.4)) + self.lobe_min = self.lobe_maj * np.random.uniform(0.2, 1) + self.lobe_ang = np.random.uniform(0, np.pi) / np.pi * 180 + + # Transform to pixel + self.lobe_maj = self.param.get_angle(self.lobe_maj) + self.lobe_min = self.param.get_angle(self.lobe_min) + lobe = [self.lobe_maj, self.lobe_min, self.lobe_ang] + + return lobe + + def gen_single_ps(self): + """ + Generate single point source, and return its data as a list. + + """ + # Redshift and luminosity + self.z, self.lumo = self.get_lumo_redshift() + self.lumo_sr = self.lumo + # angular diameter distance + self.param = PixelParams(self.z) + self.dA = self.param.dA + # W/Hz/Sr to Jy + dA = self.dA * 3.0856775814671917E+22 # Mpc to meter + self.lumo = self.lumo / dA**2 / (10.0**-24) # [Jy] + # Position + x = np.random.uniform(0, 1) + self.lat = (np.arccos(2 * x - 1) / np.pi * 180 - 90) # [deg] + self.lon = np.random.uniform(0, np.pi * 2) / np.pi * 180 # [deg] + # lobe + lobe = self.gen_lobe() + # Area + self.area = np.pi * self.lobe_maj * self.lobe_min + + ps_list = [self.z, self.dA, self.lumo, self.lat, self.lon, self.area] + ps_list.extend(lobe) + + return ps_list + + def draw_single_ps(self, freq): + """ + Designed to draw the elliptical lobes of FRI and FRII + + Prameters + --------- + nside: int and dyadic + self.ps_catalog: pandas.core.frame.DataFrame + Data of the point sources + ps_type: int + Class type of the point soruces + freq: float + frequency + """ + # Init + resolution = self.resolution / 60 # [degree] + npix = hp.nside2npix(self.nside) + hpmap = np.zeros((npix,)) + num_ps = self.ps_catalog.shape[0] + # Gen flux list + Tb_list = self.calc_Tb(freq) + ps_lobe = Tb_list[:, 1] + # Iteratively draw ps + for i in range(num_ps): + # Parameters + c_lat = self.ps_catalog['Lat (deg)'][i] # core lat [deg] + c_lon = self.ps_catalog['Lon (deg)'][i] # core lon [au.deg] + lobe_maj = self.ps_catalog['lobe_maj (rad)'][ + i] * 180 / np.pi # [deg] + lobe_min = self.ps_catalog['lobe_min (rad)'][ + i] * 180 / np.pi # [deg] + lobe_ang = self.ps_catalog['lobe_ang (deg)'][ + i] / 180 * np.pi # [rad] + # Offset to the core, refer to Willman Sec2.5.vii + offset = lobe_maj * 2 * np.random.uniform(0.2, 0.8) + # Lobe1 + lobe1_lat = (lobe_maj / 2 + offset) * np.cos(lobe_ang) + lobe1_lat = c_lat + lobe1_lat + lobe1_lon = (lobe_maj / 2 + offset) * np.sin(lobe_ang) + lobe1_lon = c_lon + lobe1_lon + # draw + # Fill with ellipse + lon, lat, gridmap = grid.make_grid_ellipse( + (lobe1_lon, lobe1_lat), (lobe_maj, lobe_min), + resolution, lobe_ang / np.pi * 180) + indices, values = grid.map_grid_to_healpix( + (lon, lat, gridmap), self.nside) + hpmap[indices] += ps_lobe[i] + + # lobe1_hotspot + lobe1_hot_lat = (lobe_maj + offset) * np.cos(lobe_ang) + lobe1_hot_lat = (c_lat + 90 + lobe1_lat) / 180 * np.pi + lobe1_hot_lon = (lobe_maj + offset) * np.sin(lobe_ang) + lobe1_hot_lon = (c_lon + lobe1_lon) / 180 * np.pi + if lobe1_hot_lat < 0: + lobe1_hot_lat += np.pi + elif lobe1_hot_lat > np.pi: + lobe1_hot_lat -= np.pi + lobe1_hot_index = hp.ang2pix( + self.nside, lobe1_hot_lat, lobe1_hot_lon) + hpmap[lobe1_hot_index] += Tb_list[i, 2] + + # Lobe2 + lobe2_lat = (lobe_maj / 2) * np.cos(lobe_ang + np.pi) + lobe2_lat = c_lat + lobe2_lat + lobe2_lon = (lobe_maj / 2) * np.sin(lobe_ang + np.pi) + lobe2_lon = c_lon + lobe2_lon + # draw + # Fill with ellipse + lon, lat, gridmap = grid.make_grid_ellipse( + (lobe2_lon, lobe2_lat), (lobe_maj, lobe_min), + resolution, lobe_ang / np.pi * 180) + indices, values = grid.map_grid_to_healpix( + (lon, lat, gridmap), self.nside) + hpmap[indices] += ps_lobe[i] + + # lobe2_hotspot + lobe2_hot_lat = (lobe_maj + offset) * np.cos(lobe_ang + np.pi) + lobe2_hot_lat = (c_lat + 90 + lobe1_lat) / 180 * np.pi + lobe2_hot_lon = (lobe_maj + offset) * np.sin(lobe_ang + np.pi) + lobe2_hot_lon = (c_lon + lobe1_lon) / 180 * np.pi + if lobe2_hot_lat < 0: + lobe2_hot_lat += np.pi + elif lobe2_hot_lat > np.pi: + lobe2_hot_lat -= np.pi + lobe2_hot_index = hp.ang2pix( + self.nside, lobe2_hot_lat, lobe2_hot_lon) + hpmap[lobe2_hot_index] += Tb_list[i, 2] + + # Core + pix_tmp = hp.ang2pix(self.nside, + (self.ps_catalog['Lat (deg)'] + 90) / + 180 * np.pi, self.ps_catalog['Lon (deg)'] / + 180 * np.pi) + ps_core = Tb_list[:, 0] + hpmap[pix_tmp] += ps_core + + return hpmap + + def draw_ps(self, freq): + """ + Read csv ps list file, and generate the healpix structure vector + with the respect frequency. + """ + # Init + num_freq = len(freq) + npix = hp.nside2npix(self.nside) + hpmaps = np.zeros((npix, num_freq)) + + # Gen ps_catalog + self.gen_catalog() + # get hpmaps + for i in range(num_freq): + hpmaps[:, i] = self.draw_single_ps(freq[i]) + + return hpmaps + + def calc_single_Tb(self, area, freq): + """ + Calculate brightness temperatur of a single ps + + Parameters + ------------ + area: `~astropy.units.Quantity` + Area of the PS, e.g., `1.0*au.sr` + freq: `~astropy.units.Quantity` + Frequency, e.g., `1.0*au.MHz` + + Return + ------ + Tb:`~astropy.units.Quantity` + Average brightness temperature, e.g., `1.0*au.K` + """ + # Init + freq_ref = 151 # [MHz] + freq = freq # [MHz] + # Luminosity at 151MHz + lumo_151 = self.lumo # [Jy] + # Calc flux + # core-to-extend ratio + ang = self.lobe_ang / 180 * np.pi + x = np.random.normal(self.xmed, 0.5) + beta = np.sqrt((self.gamma**2 - 1) / self.gamma) + B_theta = 0.5 * ((1 - beta * np.cos(ang))**-2 + + (1 + beta * np.cos(ang))**-2) + ratio_obs = 10**x * B_theta + # Core + lumo_core = ratio_obs / (1 + ratio_obs) * lumo_151 + a0 = (np.log10(lumo_core) - 0.07 * + np.log10(freq_ref * 10.0E-3) + + 0.29 * np.log10(freq_ref * 10.0E-3) * + np.log10(freq_ref * 10.0E-3)) + lgs = (a0 + 0.07 * np.log10(freq * 10.0E-3) - 0.29 * + np.log10(freq * 10.0E-3) * + np.log10(freq * 10.0E-3)) + flux_core = 10**lgs # [Jy] + # core area + npix = hp.nside2npix(self.nside) + core_area = 4 * np.pi / npix # [sr] + Tb_core = convert.Fnu_to_Tb_fast(flux_core, core_area, freq) # [K] + # lobe + lumo_lobe = lumo_151 * (1 - ratio_obs) / (1 + ratio_obs) # [Jy] + flux_lobe = (freq / freq_ref)**(-0.75) * lumo_lobe + Tb_lobe = convert.Fnu_to_Tb_fast(flux_lobe, area, freq) # [K] + + # hotspots + # Willman Eq. (3) + f_hs = (0.4 * (np.log10(self.lumo_sr) - 25.5) + + np.random.uniform(-0.5, 0.5)) + Tb_hotspot = Tb_lobe * (1 + f_hs) + Tb = [Tb_core, Tb_lobe, Tb_hotspot] + + return Tb + + def calc_Tb(self, freq): + """ + Calculate the surface brightness temperature of the point sources. + + Parameters + ------------ + area: `~astropy.units.Quantity` + Area of the PS, e.g., `1.0*au.sr` + freq: `~astropy.units.Quantity` + Frequency, e.g., `1.0*au.MHz` + + Return + ------ + Tb_list: list + Point sources brightness temperature + """ + # Tb_list + num_ps = self.ps_catalog.shape[0] + Tb_list = np.zeros((num_ps, 3)) + # Iteratively calculate Tb + for i in range(num_ps): + ps_area = self.ps_catalog['Area (sr)'][i] # [sr] + Tb_list[i, :] = self.calc_single_Tb(ps_area, freq) + + return Tb_list diff --git a/fg21sim/extragalactic/pointsources/pointsources.py b/fg21sim/extragalactic/pointsources/pointsources.py new file mode 100644 index 0000000..85f6ec4 --- /dev/null +++ b/fg21sim/extragalactic/pointsources/pointsources.py @@ -0,0 +1,126 @@ +# Copyright (c) 2016 Zhixian MA +# MIT license + +""" +Extragalactic point sources (ps) simulation +""" + +import logging +import numpy as np +import healpy as hp +from collections import OrderedDict + +from .starforming import StarForming +from .starbursting import StarBursting +from .radioquiet import RadioQuiet +from .fr1 import FRI +from .fr2 import FRII + + +logger = logging.getLogger(__name__) + + +class PointSources: + """ + This class namely pointsource is designed to generate PS catalogs, + read csv format PS lists, calculate the flux and surface brightness + of the sources at different frequencies, and then ouput hpmaps + + Parameters + ---------- + configs: ConfigManager object + An 'ConfigManager' object contains default and user configurations. + For more details, see the example config specification. + + Functions + --------- + preprocessing + Generate the ps catalogs for each type. + simulate_frequency + Simualte point sources at provivded frequency + simulate + Simulate and project PSs to the healpix map. + postprocessing + Save catalogs + """ + PSCOMPONENTS_ALL = OrderedDict([ + ("starforming", StarForming), + ("starbursting", StarBursting), + ("radioquiet", RadioQuiet), + ("FRI", FRI), + ("FRII", FRII), + ]) + + def __init__(self, configs): + self.configs = configs + self._set_configs() + self.pscomps = OrderedDict() + for comp in self.pscomps_id: + logger.info("Initlalize PS component: {0}".format(comp)) + comp_type = self.PSCOMPONENTS_ALL[comp] + self.pscomps[comp] = comp_type(configs) + logger.info("Done initlalize %d PS components!" % + len(self.pscomps)) + + def _set_configs(self): + """Load configs and set the attributes""" + # Prefix of simulated point sources + self.pscomps_id = self.configs.getn("extragalactic/pscomponents") + if self.pscomps_id is None: + self.pscomps_id = ['starforming', 'starbursting', 'radioquiet', + 'FRI', 'FRII'] + print(self.pscomps_id) + # nside of the healpix cell + self.nside = self.configs.getn("common/nside") + # save flag + self.save = self.configs.getn("extragalactic/pointsources/save") + + def preprocess(self): + """Preprocess and generate the catalogs""" + logger.info("Generating PS catalogs...") + # Gen ps_catalog + for pscomp_obj in self.pscomps.values(): + pscomp_obj.gen_catalog() + logger.info("Generating PS catalogs done!") + + def simulate_frequency(self, freq): + """Simulate the point sources and output hpmaps""" + npix = hp.nside2npix(self.nside) + hpmap_f = np.zeros((npix,)) + # Projecting + logger.info("Generating PS hpmaps...") + for pscomp_obj in self.pscomps.values(): + hpmap_f += pscomp_obj.draw_single_ps(freq) + logger.info("Generating PS hpmaps done!") + + return hpmap_f + + def simulate(self, frequencies): + """Simulate the emission (HEALPix) maps of all Galactic SNRs for + every specified frequency. + + Parameters + ---------- + frequency : list[float] + List of frequencies (unit: `self.freq_unit`) where the + simulation performed. + + Returns + ------- + hpmaps : list[1D `~numpy.ndarray`] + List of HEALPix maps (in RING ordering) at each frequency. + """ + hpmaps = [] + for f in np.array(frequencies, ndmin=1): + hpmap_f = self.simulate_frequency(f) + hpmaps.append(hpmap_f) + return hpmaps + + def postprocess(self): + """Perform the post-simulation operations before the end.""" + # Save the catalog actually used in the simulation + if self.save: + logger.info("Saving simulated catalogs...") + for pscomp_obj in self.pscomps.values(): + pscomp_obj.save_as_csv() + logger.info("Saving simulated catalogs done!") diff --git a/fg21sim/extragalactic/pointsources/psparams.py b/fg21sim/extragalactic/pointsources/psparams.py new file mode 100644 index 0000000..7c3a93a --- /dev/null +++ b/fg21sim/extragalactic/pointsources/psparams.py @@ -0,0 +1,77 @@ +# Copyright (c) 2016 Zhixian MA +# MIT license + +""" +Basic parameters to be used in point sources simulation. + +References +---------- +[1] Angular diameter distance, + https://en.wikipedia.org/wiki/Angular_diameter_distance +[2] FlatLambdaCDM, + {http://docs.astropy.org/en/stable/api/astropy.cosmology. + FlatLambdaCDM.html#astropy.cosmology.FlatLambdaCDM} +""" + +from astropy.cosmology import FlatLambdaCDM + + +class PixelParams(): + """ + A class to transform cosmology distance to angles or pixels. + Parameters + ------------ + H0: float + The hubble constant at z = 0, whose unit is km/s/Mpc + Om0: float + The total matter density. + ang_res: float + Angular resolution, i.e. degree per pixel. (May be useless) + ang_total: list + Total angles of the simulated sky region,whose unit is degree (\deg) + z : float + Redshift + scale: float + The real object scale. + Example + ------------ + >>> PixelParams = PixelParams(img_size=(1024,1024),ang_total=(5,5)) + """ + + # Hubble constant at z = 0 + H0 = 71.0 + # Omega0, total matter density + Om0 = 0.27 + # Redshift + z = 0.0 + # Cosmology calculator + cosmo = 0.0 + # angular diameter distance, [Mpc] + dA = 0.0 + # angular resolution + ang_res = 0.0 + + def __init__(self, z=0.0): + self.z = z + self.cosmo = FlatLambdaCDM(H0=self.H0, Om0=self.Om0) + + # angular diameter distance, [Mpc] + self.dA = self.cosmo.angular_diameter_distance(self.z).value + + def get_angle(self, scale=1.0): + """ + Input real object scale, and output the respect observed + angle, and pixels. + """ + ang = scale / self.dA # [rac] + + return ang + + def get_scale(self, ang=1.0): + """ + Input real observed scale, and output the respect + real object scale, and pixels. + """ + scale = ang * self.dA # [Mpc] + + return scale diff --git a/fg21sim/extragalactic/pointsources/radioquiet.py b/fg21sim/extragalactic/pointsources/radioquiet.py new file mode 100644 index 0000000..2d2eefe --- /dev/null +++ b/fg21sim/extragalactic/pointsources/radioquiet.py @@ -0,0 +1,189 @@ +# Copyright (c) 2016 Zhixian MA +# MIT license + +import numpy as np +import healpy as hp + +from .base import BasePointSource +from ...utils import convert + + +class RadioQuiet(BasePointSource): + + def __init__(self, configs): + super().__init__(configs) + self._set_configs() + # Number density matrix + self.rho_mat = self.calc_number_density() + # Cumulative distribution of z and lumo + self.cdf_z, self.cdf_lumo = self.calc_cdf() + + def _set_configs(self): + """Load the configs and set the corresponding class attributes""" + super()._set_configs() + pscomp = "extragalactic/pointsources/radioquiet/" + # point sources amount + self.num_ps = self.configs.getn(pscomp+"numps") + # prefix + self.prefix = self.configs.getn(pscomp+"prefix") + # redshift bin + z_type = self.configs.getn(pscomp+"z_type") + if z_type == 'custom': + start = self.configs.getn(pscomp+"z_start") + stop = self.configs.getn(pscomp+"z_stop") + step = self.configs.getn(pscomp+"z_step") + self.zbin = np.arange(start, stop + step, step) + else: + self.zbin = np.arange(0.1, 10, 0.05) + # luminosity bin + lumo_type = self.configs.getn(pscomp+"lumo_type") + if lumo_type == 'custom': + start = self.configs.getn(pscomp+"lumo_start") + stop = self.configs.getn(pscomp+"lumo_stop") + step = self.configs.getn(pscomp+"lumo_step") + self.lumobin = np.arange(start, stop + step, step) + else: + self.lumobin = np.arange(18.7, 25.7, 0.1) # [W/Hz/sr] + + def calc_number_density(self): + """ + Calculate number density rho(lumo,z) of FRI + + References + ------------ + [1] Wilman et al., + "A semi-empirical simulation of the extragalactic radio continuum + sky for next generation radio telescopes", + 2008, MNRAS, 388, 1335-1348. + http://adsabs.harvard.edu/abs/2008MNRAS.388.1335W + + Returns + ------- + rho_mat: np.ndarray + Number density matris (joint-distribution of luminosity and + reshift). + """ + # Init + rho_mat = np.zeros((len(self.lumobin), len(self.zbin))) + # Parameters + # Refer to Willman's section 2.4 + alpha = 0.7 # spectral index + lumo_star = 10.0**21.3 # critical luminosity at 1400MHz + rho_l0 = 10.0**(-7) # normalization constant + z1 = 1.9 # cut-off redshift + k1 = -3.27 # index of space density revolution + # Calculation + for i, z in enumerate(self.zbin): + if z <= z1: + rho_mat[:, i] = ((rho_l0 * (10**self.lumobin / lumo_star) ** + (-alpha) * + np.exp(-10**self.lumobin / lumo_star)) * + (1 + z)**k1) + else: + rho_mat[:, i] = ((rho_l0 * (10**self.lumobin / lumo_star) ** + (-alpha) * + np.exp(-10**self.lumobin / lumo_star)) * + (1 + z1)**k1) + return rho_mat + + def draw_single_ps(self, freq): + """ + Designed to draw the radio quiet AGN + + Parameters + ---------- + ImgMat: np.ndarray + Two dimensional matrix, to describe the image + self.ps_catalog: pandas.core.frame.DataFrame + Data of the point sources + freq: float + frequency + """ + # Init + npix = hp.nside2npix(self.nside) + hpmap = np.zeros((npix,)) + # Gen Tb list + Tb_list = self.calc_Tb(freq) + # Iteratively draw the ps + num_ps = self.ps_catalog.shape[0] + for i in range(num_ps): + # Angle to pix + lat = (self.ps_catalog['Lat (deg)'] + 90) / 180 * np.pi + lon = self.ps_catalog['Lon (deg)'] / 180 * np.pi + pix = hp.ang2pix(self.nside, lat, lon) + # Gen hpmap + hpmap[pix] += Tb_list[i] + + return hpmap + + def draw_ps(self, freq): + """ + Read csv ps list file, and generate the healpix structure vector + with the respect frequency. + """ + # Init + num_freq = len(freq) + npix = hp.nside2npix(self.nside) + hpmaps = np.zeros((npix, num_freq)) + + # Gen ps_catalog + self.gen_catalog() + # get hpmaps + for i in range(num_freq): + hpmaps[:, i] = self.draw_single_ps(freq[i]) + + return hpmaps + + def calc_single_Tb(self, area, freq): + """ + Calculate brightness temperatur of a single ps + + Parameters + ------------ + area: `~astropy.units.Quantity` + Area of the PS, e.g., `1.0*au.sr2` + freq: `~astropy.units.Quantity` + Frequency, e.g., `1.0*au.MHz` + + Return + ------ + Tb:`~astropy.units.Quantity` + Average brightness temperature, e.g., `1.0*au.K` + """ + # Init + freq_ref = 1400 # [MHz] + freq = freq # [MHz] + # Luminosity at 1400MHz + lumo_1400 = self.lumo # [Jy] + # Calc flux + flux = (freq / freq_ref)**(-0.7) * lumo_1400 + # Calc brightness temperature + Tb = convert.Fnu_to_Tb_fast(flux, area, freq) + + return Tb + + def calc_Tb(self, freq): + """ + Calculate the surface brightness temperature of the point sources. + + Parameters + ------------ + area: `~astropy.units.Quantity` + Area of the PS, e.g., `1.0*au.sr` + freq: `~astropy.units.Quantity` + Frequency, e.g., `1.0*au.MHz` + + Return + ------ + Tb_list: list + Point sources brightness temperature + """ + # Tb_list + num_ps = self.ps_catalog.shape[0] + Tb_list = np.zeros((num_ps,)) + # Iteratively calculate Tb + for i in range(num_ps): + ps_area = self.ps_catalog['Area (sr)'][i] # [sr] + Tb_list[i] = self.calc_single_Tb(ps_area, freq) + + return Tb_list diff --git a/fg21sim/extragalactic/pointsources/starbursting.py b/fg21sim/extragalactic/pointsources/starbursting.py new file mode 100644 index 0000000..8b5c597 --- /dev/null +++ b/fg21sim/extragalactic/pointsources/starbursting.py @@ -0,0 +1,238 @@ +# Copyright (c) 2016 Zhixian MA +# MIT license + +import numpy as np +import healpy as hp + +from .psparams import PixelParams +from .base import BasePointSource +from ...utils import grid +from ...utils import convert + + +class StarBursting(BasePointSource): + + """ + Generate star forming point sources, inheritate from PointSource class. + """ + # Init + + def __init__(self, configs): + super().__init__(configs) + self.columns.append('radius (rad)') + self.nCols = len(self.columns) + self._set_configs() + # Number density matrix + self.rho_mat = self.calc_number_density() + # Cumulative distribution of z and lumo + self.cdf_z, self.cdf_lumo = self.calc_cdf() + + def _set_configs(self): + """Load the configs and set the corresponding class attributes""" + super()._set_configs() + pscomp = "extragalactic/pointsources/starbursting/" + # point sources amount + self.num_ps = self.configs.getn(pscomp+"numps") + # prefix + self.prefix = self.configs.getn(pscomp+"prefix") + # redshift bin + z_type = self.configs.getn(pscomp+"z_type") + if z_type == 'custom': + start = self.configs.getn(pscomp+"z_start") + stop = self.configs.getn(pscomp+"z_stop") + step = self.configs.getn(pscomp+"z_step") + self.zbin = np.arange(start, stop + step, step) + else: + self.zbin = np.arange(0.1, 10, 0.05) + # luminosity bin + lumo_type = self.configs.getn(pscomp+"lumo_type") + if lumo_type == 'custom': + start = self.configs.getn(pscomp+"lumo_start") + stop = self.configs.getn(pscomp+"lumo_stop") + step = self.configs.getn(pscomp+"lumo_step") + self.lumobin = np.arange(start, stop + step, step) + else: + self.lumobin = np.arange(21, 27, 0.1) # [W/Hz/sr] + + def calc_number_density(self): + """ + Calculate number density rho(lumo,z) of FRI + + References + ---------- + [1] Wilman et al., + "A semi-empirical simulation of the extragalactic radio continuum + sky for next generation radio telescopes", + 2008, MNRAS, 388, 1335-1348. + http://adsabs.harvard.edu/abs/2008MNRAS.388.1335W + + Returns + ------- + rho_mat: np.ndarray + Number density matris (joint-distribution of luminosity and + reshift). + """ + # Init + rho_mat = np.zeros((len(self.lumobin), len(self.zbin))) + # Parameters + # Refer to Willman's section 2.4 + alpha = 0.7 # spectral index + lumo_star = 10.0**22 # critical luminosity at 1400MHz + rho_l0 = 10.0**(-7) # normalization constant + z1 = 1.5 # cut-off redshift + k1 = 3.1 # index of space density revolution + # Calculation + for i, z in enumerate(self.zbin): + if z <= z1: + rho_mat[:, i] = ((rho_l0 * (10**self.lumobin / lumo_star) ** + (-alpha) * + np.exp(-10**self.lumobin / lumo_star)) * + (1 + z)**k1) + else: + rho_mat[:, i] = ((rho_l0 * (10**self.lumobin / lumo_star) ** + (-alpha) * + np.exp(-10**self.lumobin / lumo_star)) * + (1 + z1)**k1) + return rho_mat + + def get_radius(self): + if self.z <= 1.5: + self.radius = (1 + self.z)**2.5 * 1e-3 / 2 # [Mpc] + else: + self.radius = 10 * 1e-3 / 2 # [Mpc] + + return self.radius + + def gen_single_ps(self): + """ + Generate single point source, and return its data as a list. + """ + # Redshift and luminosity + self.z, self.lumo = self.get_lumo_redshift() + # angular diameter distance + self.param = PixelParams(self.z) + self.dA = self.param.dA + # W/Hz/Sr to Jy + dA = self.dA * 3.0856775814671917E+22 # Mpc to meter + self.lumo = self.lumo / dA**2 / (10.0**-24) # [Jy] + # Position + x = np.random.uniform(0, 1) + self.lat = (np.arccos(2 * x - 1) / np.pi * 180 - 90) # [deg] + self.lon = np.random.uniform(0, np.pi * 2) / np.pi * 180 # [deg] + # Radius + self.radius = self.param.get_angle(self.get_radius()) # [rad] + # Area + self.area = np.pi * self.radius**2 # [sr] + + ps_list = [self.z, self.dA, self.lumo, self.lat, self.lon, + self.area, self.radius] + + return ps_list + + def draw_single_ps(self, freq): + """ + Designed to draw the circular star forming and star bursting ps. + + Prameters + --------- + nside: int and dyadic + number of sub pixel in a cell of the healpix structure + self.ps_catalog: pandas.core.frame.DataFrame + Data of the point sources + freq: float + frequency + """ + # Init + npix = hp.nside2npix(self.nside) + hpmap = np.zeros((npix,)) + # Gen Tb list + Tb_list = self.calc_Tb(freq) + # Iteratively draw the ps + num_ps = self.ps_catalog.shape[0] + resolution = self.resolution / 60 # [degree] + for i in range(num_ps): + # grid + ps_radius = self.ps_catalog['radius (rad)'][i] # [rad] + ps_radius = ps_radius * 180 / np.pi # radius[rad] + c_lat = self.ps_catalog['Lat (deg)'][i] # core_lat [deg] + c_lon = self.ps_catalog['Lon (deg)'][i] # core_lon [deg] + # Fill with circle + lon, lat, gridmap = grid.make_grid_ellipse( + (c_lon, c_lat), (2 * ps_radius, 2 * ps_radius), resolution) + indices, values = grid.map_grid_to_healpix( + (lon, lat, gridmap), self.nside) + hpmap[indices] += Tb_list[i] + + return hpmap + + def draw_ps(self, freq): + """ + Read csv ps list file, and generate the healpix structure vector + with the respect frequency. + """ + # Init + num_freq = len(freq) + npix = hp.nside2npix(self.nside) + hpmaps = np.zeros((npix, num_freq)) + + # Gen ps_catalog + self.gen_catalog() + # get hpmaps + for i in range(num_freq): + hpmaps[:, i] = self.draw_single_ps(freq[i]) + + return hpmaps + + def calc_single_Tb(self, area, freq): + """ + Calculate brightness temperatur of a single ps + + Parameters + ------------ + area: `~astropy.units.Quantity` + Area of the PS, e.g., `1.0*au.sr2` + freq: `~astropy.units.Quantity` + Frequency, e.g., `1.0*au.MHz` + + Return + ------ + Tb:`~astropy.units.Quantity` + Average brightness temperature, e.g., `1.0*au.K` + """ + # Init + freq_ref = 1400 # [MHz] + freq = freq # [MHz] + # Luminosity at 1400MHz + lumo_1400 = self.lumo # [Jy] + # Calc flux + flux = (freq / freq_ref)**(-0.7) * lumo_1400 + # Calc brightness temperature + Tb = convert.Fnu_to_Tb_fast(flux, area, freq) + + return Tb + + def calc_Tb(self, freq): + """ + Calculate the surface brightness temperature of the point sources. + + Parameters + ------------ + area: `~astropy.units.Quantity` + Area of the PS, e.g., `1.0*au.sr` + freq: `~astropy.units.Quantity` + Frequency, e.g., `1.0*au.MHz` + + Return + ------ + Tb_list: list + Point sources brightness temperature + """ + # Tb_list + num_ps = self.ps_catalog.shape[0] + Tb_list = np.zeros((num_ps,)) + # Iteratively calculate Tb + for i in range(num_ps): + ps_area = self.ps_catalog['Area (sr)'][i] # [sr] + Tb_list[i] = self.calc_single_Tb(ps_area, freq) + + return Tb_list diff --git a/fg21sim/extragalactic/pointsources/starforming.py b/fg21sim/extragalactic/pointsources/starforming.py new file mode 100644 index 0000000..5d518d8 --- /dev/null +++ b/fg21sim/extragalactic/pointsources/starforming.py @@ -0,0 +1,257 @@ +# Copyright (c) 2016 Zhixian MA +# MIT license + +import numpy as np +import healpy as hp + +# from .psparams import PixelParams +from .base import BasePointSource +from ...utils import grid +from ...utils import convert +from .psparams import PixelParams + + +class StarForming(BasePointSource): + """ + Generate star forming point sources, inheritate from PointSource class. + + Reference + --------- + [1] Fast cirles drawing + https://github.com/liweitianux/fg21sim/fg21sim/utils/draw.py + https://github.com/liweitianux/fg21sim/fg21sim/utils/grid.py + """ + + def __init__(self, configs): + super().__init__(configs) + self.columns.append('radius (rad)') + self.nCols = len(self.columns) + self._set_configs() + # Number density matrix + self.rho_mat = self.calc_number_density() + # Cumulative distribution of z and lumo + self.cdf_z, self.cdf_lumo = self.calc_cdf() + + def _set_configs(self): + """ Load the configs and set the corresponding class attributes""" + super()._set_configs() + pscomp = "extragalactic/pointsources/starforming/" + # point sources amount + self.num_ps = self.configs.getn(pscomp+"numps") + # prefix + self.prefix = self.configs.getn(pscomp+"prefix") + # redshift bin + z_type = self.configs.getn(pscomp+"z_type") + if z_type == 'custom': + start = self.configs.getn(pscomp+"z_start") + stop = self.configs.getn(pscomp+"z_stop") + step = self.configs.getn(pscomp+"z_step") + self.zbin = np.arange(start, stop + step, step) + else: + self.zbin = np.arange(0.1, 10, 0.05) + # luminosity bin + lumo_type = self.configs.getn(pscomp+"lumo_type") + if lumo_type == 'custom': + start = self.configs.getn(pscomp+"lumo_start") + stop = self.configs.getn(pscomp+"lumo_stop") + step = self.configs.getn(pscomp+"lumo_step") + self.lumobin = np.arange(start, stop + step, step) + else: + self.lumobin = np.arange(17, 25.5, 0.1) # [W/Hz/sr] + + def calc_number_density(self): + """ + Calculate number density rho(lumo,z) of FRI + + References + ---------- + [1] Wilman et al., + "A semi-empirical simulation of the extragalactic radio continuum + sky for next generation radio telescopes", + 2008, MNRAS, 388, 1335-1348. + http://adsabs.harvard.edu/abs/2008MNRAS.388.1335W + + Returns + ------- + rho_mat: np.ndarray + Number density matris (joint-distribution of luminosity and + reshift). + """ + # Init + rho_mat = np.zeros((len(self.lumobin), len(self.zbin))) + # Parameters + # Refer to Willman's section 2.4 + alpha = 0.7 # spectral index + lumo_star = 10.0**22 # critical luminosity at 1400MHz + rho_l0 = 10.0**(-7) # normalization constant + z1 = 1.5 # cut-off redshift + k1 = 3.1 # index of space density revolution + # Calculation + for i, z in enumerate(self.zbin): + if z <= z1: + rho_mat[:, i] = (rho_l0 * (10**self.lumobin / lumo_star) ** + (-alpha) * np.exp(-10**self.lumobin / + lumo_star) * (1 + z)**k1) + else: + rho_mat[:, i] = (rho_l0 * (10**self.lumobin / lumo_star) ** + (-alpha) * np.exp(-10**self.lumobin / + lumo_star) * (1 + z1)**k1) + + return rho_mat + + def get_radius(self): + """ + Generate the disc diameter of normal starforming galaxies. + + Reference + --------- + [1] Wilman et al., Eq(7-9), + "A semi-empirical simulation of the extragalactic radio continuum + sky for next generation radio telescopes", + 2008, MNRAS, 388, 1335-1348. + http://adsabs.harvard.edu/abs/2008MNRAS.388.1335W + """ + # Willman Eq. (8) + delta = np.random.normal(0, 0.3) + log_M_HI = 0.44 * np.log10(self.lumo) + 0.48 + delta + # Willman Eq. (7) + log_D_HI = ((log_M_HI - (6.52 + np.random.uniform(-0.06, 0.06))) / + 1.96 + np.random.uniform(-0.04, 0.04)) + # Willman Eq. (9) + log_D = log_D_HI - 0.23 - np.log10(1 + self.z) + self.radius = 10**log_D / 2 * 1e-3 # [Mpc] + return self.radius + + def gen_single_ps(self): + """ + Generate single point source, and return its data as a list. + """ + # Redshift and luminosity + self.z, self.lumo = self.get_lumo_redshift() + # angular diameter distance + self.param = PixelParams(self.z) + self.dA = self.param.dA + # Radius + self.radius = self.param.get_angle(self.get_radius()) # [rad] + # W/Hz/Sr to Jy + dA = self.dA * 3.0856775814671917E+22 # Mpc to meter + self.lumo = self.lumo / dA**2 / (10.0**-24) # [Jy] + # Position + x = np.random.uniform(0, 1) + self.lat = (np.arccos(2 * x - 1) / np.pi * 180 - 90) # [deg] + self.lon = np.random.uniform(0, np.pi * 2) / np.pi * 180 # [deg] + # Area + self.area = np.pi * self.radius**2 # [sr] ? + + ps_list = [self.z, self.dA, self.lumo, self.lat, self.lon, + self.area, self.radius] + + return ps_list + + def draw_single_ps(self, freq): + """ + Designed to draw the circular star forming and star bursting ps. + + Prameters + --------- + nside: int and dyadic + number of sub pixel in a cell of the healpix structure + self.ps_catalog: pandas.core.frame.DataFrame + Data of the point sources + freq: float + frequency + """ + # Init + npix = hp.nside2npix(self.nside) + hpmap = np.zeros((npix,)) + # Gen Tb list + Tb_list = self.calc_Tb(freq) + # Iteratively draw the ps + num_ps = self.ps_catalog.shape[0] + resolution = self.resolution / 60 # [degree] + for i in range(num_ps): + # grid + ps_radius = self.ps_catalog['radius (rad)'][i] # [rad] + ps_radius = ps_radius * 180 / np.pi # radius [deg] + c_lat = self.ps_catalog['Lat (deg)'][i] # core_lat [deg] + c_lon = self.ps_catalog['Lon (deg)'][i] # core_lon [deg] + # Fill with circle + lon, lat, gridmap = grid.make_grid_ellipse( + (c_lon, c_lat), (2 * ps_radius, 2 * ps_radius), resolution) + indices, values = grid.map_grid_to_healpix( + (lon, lat, gridmap), self.nside) + hpmap[indices] += Tb_list[i] + + return hpmap + + def draw_ps(self, freq): + """ + Read csv ps list file, and generate the healpix structure vector + with the respect frequency. + """ + # Init + num_freq = len(freq) + npix = hp.nside2npix(self.nside) + hpmaps = np.zeros((npix, num_freq)) + + # Gen ps_catalog + self.gen_catalog() + # get hpmaps + for i in range(num_freq): + hpmaps[:, i] = self.draw_single_ps(freq[i]) + + return hpmaps + + def calc_single_Tb(self, area, freq): + """ + Calculate brightness temperatur of a single ps + + Parameters + ------------ + area: `~astropy.units.Quantity` + Area of the PS, e.g., `1.0*au.sr2` + freq: `~astropy.units.Quantity` + Frequency, e.g., `1.0*au.MHz` + + Return + ------ + Tb:`~astropy.units.Quantity` + Average brightness temperature, e.g., `1.0*au.K` + """ + # Init + freq_ref = 1400 # [MHz] + freq = freq # [MHz] + # Luminosity at 1400MHz + lumo_1400 = self.lumo # [Jy] + # Calc flux + flux = (freq / freq_ref)**(-0.7) * lumo_1400 + # Calc brightness temperature + Tb = convert.Fnu_to_Tb_fast(flux, area, freq) + + return Tb + + def calc_Tb(self, freq): + """ + Calculate the surface brightness temperature of the point sources. + + Parameters + ------------ + area: `~astropy.units.Quantity` + Area of the PS, e.g., `1.0*au.sr` + freq: `~astropy.units.Quantity` + Frequency, e.g., `1.0*au.MHz` + + Return + ------ + Tb_list: list + Point sources brightness temperature + """ + # Tb_list + num_ps = self.ps_catalog.shape[0] + Tb_list = np.zeros((num_ps,)) + # Iteratively calculate Tb + for i in range(num_ps): + ps_area = self.ps_catalog['Area (sr)'][i] # [sr] + Tb_list[i] = self.calc_single_Tb(ps_area, freq) + + return Tb_list diff --git a/fg21sim/foregrounds.py b/fg21sim/foregrounds.py index 9788d2c..e3e0f62 100644 --- a/fg21sim/foregrounds.py +++ b/fg21sim/foregrounds.py @@ -25,6 +25,7 @@ from .galactic import (Synchrotron as GalacticSynchrotron, FreeFree as GalacticFreeFree, SuperNovaRemnants as GalacticSNR) from .extragalactic import GalaxyClusters as EGGalaxyClusters +from .extragalactic import PointSources as ExtragalacticPointSources from .utils import write_fits_healpix @@ -67,6 +68,10 @@ class Foregrounds: ("galactic/freefree", GalacticFreeFree), ("galactic/snr", GalacticSNR), ("extragalactic/clusters", EGGalaxyClusters), + ("galactic/synchrotron", GalacticSynchrotron), + ("galactic/freefree", GalacticFreeFree), + ("galactic/snr", GalacticSNR), + ("extragalactic/pointsources", ExtragalacticPointSources), ]) def __init__(self, configs): -- cgit v1.2.2