aboutsummaryrefslogtreecommitdiffstats
path: root/fg21sim/extragalactic/pointsources/radioquiet.py
diff options
context:
space:
mode:
Diffstat (limited to 'fg21sim/extragalactic/pointsources/radioquiet.py')
-rw-r--r--fg21sim/extragalactic/pointsources/radioquiet.py189
1 files changed, 189 insertions, 0 deletions
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 <zxma_sjtu@qq.com>
+# 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