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