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/extragalactic/pointsources/radioquiet.py | 189 +++++++++++++++++++++++ 1 file changed, 189 insertions(+) create mode 100644 fg21sim/extragalactic/pointsources/radioquiet.py (limited to 'fg21sim/extragalactic/pointsources/radioquiet.py') 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 -- cgit v1.2.2