diff options
Diffstat (limited to 'fg21sim/extragalactic/pointsources/fr1.py')
| -rw-r--r-- | fg21sim/extragalactic/pointsources/fr1.py | 355 |
1 files changed, 355 insertions, 0 deletions
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 <zxma_sjtu@qq.com> +# 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 |