aboutsummaryrefslogtreecommitdiffstats
path: root/fg21sim/extragalactic/pointsources
diff options
context:
space:
mode:
Diffstat (limited to 'fg21sim/extragalactic/pointsources')
-rw-r--r--fg21sim/extragalactic/pointsources/__init__.py4
-rw-r--r--fg21sim/extragalactic/pointsources/base.py203
-rw-r--r--fg21sim/extragalactic/pointsources/fr1.py357
-rw-r--r--fg21sim/extragalactic/pointsources/fr2.py382
-rw-r--r--fg21sim/extragalactic/pointsources/pointsources.py128
-rw-r--r--fg21sim/extragalactic/pointsources/psparams.py79
-rw-r--r--fg21sim/extragalactic/pointsources/radioquiet.py190
-rw-r--r--fg21sim/extragalactic/pointsources/starbursting.py239
-rw-r--r--fg21sim/extragalactic/pointsources/starforming.py258
9 files changed, 0 insertions, 1840 deletions
diff --git a/fg21sim/extragalactic/pointsources/__init__.py b/fg21sim/extragalactic/pointsources/__init__.py
deleted file mode 100644
index ea60028..0000000
--- a/fg21sim/extragalactic/pointsources/__init__.py
+++ /dev/null
@@ -1,4 +0,0 @@
-# Copyright (c) 2016 Zhixian MA <zxma_sjtu@qq.com>
-# MIT license
-
-from .pointsources import PointSources # noqa: F401
diff --git a/fg21sim/extragalactic/pointsources/base.py b/fg21sim/extragalactic/pointsources/base.py
deleted file mode 100644
index 1544f5b..0000000
--- a/fg21sim/extragalactic/pointsources/base.py
+++ /dev/null
@@ -1,203 +0,0 @@
-# Copyright (c) 2016 Zhixian MA <zxma_sjtu@qq.com>
-# 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
deleted file mode 100644
index 4d337eb..0000000
--- a/fg21sim/extragalactic/pointsources/fr1.py
+++ /dev/null
@@ -1,357 +0,0 @@
-# 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: float
- Area of the PS
- Unit: [arcsec^2]
- 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)
- sr_to_arcsec2 = (np.rad2deg(1) * 3600) ** 2 # [sr] -> [arcsec^2]
- core_area = 4 * np.pi / npix * sr_to_arcsec2 # [arcsec^2]
- Tb_core = convert.Fnu_to_Tb(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(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
- ------------
- 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))
- sr_to_arcsec2 = (np.rad2deg(1) * 3600) ** 2 # [sr] -> [arcsec^2]
- # Iteratively calculate Tb
- for i in range(num_ps):
- ps_area = self.ps_catalog['Area (sr)'][i] # [sr]
- area = ps_area * sr_to_arcsec2
- Tb_list[i, :] = self.calc_single_Tb(area, freq)
-
- return Tb_list
diff --git a/fg21sim/extragalactic/pointsources/fr2.py b/fg21sim/extragalactic/pointsources/fr2.py
deleted file mode 100644
index ec565d1..0000000
--- a/fg21sim/extragalactic/pointsources/fr2.py
+++ /dev/null
@@ -1,382 +0,0 @@
-# Copyright (c) 2016 Zhixian MA <zxma_sjtu@qq.com>
-# 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: float
- Area of the PS
- Unit: [arcsec^2]
- 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)
- sr_to_arcsec2 = (np.rad2deg(1) * 3600) ** 2 # [sr] -> [arcsec^2]
- core_area = 4 * np.pi / npix * sr_to_arcsec2 # [arcsec^2]
- Tb_core = convert.Fnu_to_Tb(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(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
- ------------
- 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))
- sr_to_arcsec2 = (np.rad2deg(1) * 3600) ** 2 # [sr] -> [arcsec^2]
- # Iteratively calculate Tb
- for i in range(num_ps):
- ps_area = self.ps_catalog['Area (sr)'][i] # [sr]
- area = ps_area * sr_to_arcsec2
- Tb_list[i, :] = self.calc_single_Tb(area, freq)
-
- return Tb_list
diff --git a/fg21sim/extragalactic/pointsources/pointsources.py b/fg21sim/extragalactic/pointsources/pointsources.py
deleted file mode 100644
index a921e03..0000000
--- a/fg21sim/extragalactic/pointsources/pointsources.py
+++ /dev/null
@@ -1,128 +0,0 @@
-# Copyright (c) 2016 Zhixian MA <zxma_sjtu@qq.com>
-# 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!")
- # XXX/TODO:
- # * Output the HEALPix map to file
- # * Also return the path to the output file
- return (hpmap_f, None)
-
- 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
deleted file mode 100644
index 823af72..0000000
--- a/fg21sim/extragalactic/pointsources/psparams.py
+++ /dev/null
@@ -1,79 +0,0 @@
-# Copyright (c) 2016 Zhixian MA <zxma_sjtu@qq.com>
-# 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
- 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
deleted file mode 100644
index 0fbb6b2..0000000
--- a/fg21sim/extragalactic/pointsources/radioquiet.py
+++ /dev/null
@@ -1,190 +0,0 @@
-# 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: float
- Area of the PS
- Unit: [arcsec^2]
- 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(flux, area, freq)
-
- return Tb
-
- def calc_Tb(self, freq):
- """
- Calculate the surface brightness temperature of the point sources.
-
- Parameters
- ------------
- 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,))
- sr_to_arcsec2 = (np.rad2deg(1) * 3600) ** 2 # [sr] -> [arcsec^2]
- # Iteratively calculate Tb
- for i in range(num_ps):
- ps_area = self.ps_catalog['Area (sr)'][i] # [sr]
- area = ps_area * sr_to_arcsec2
- Tb_list[i] = self.calc_single_Tb(area, freq)
-
- return Tb_list
diff --git a/fg21sim/extragalactic/pointsources/starbursting.py b/fg21sim/extragalactic/pointsources/starbursting.py
deleted file mode 100644
index 897d53f..0000000
--- a/fg21sim/extragalactic/pointsources/starbursting.py
+++ /dev/null
@@ -1,239 +0,0 @@
-# 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 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: float
- Area of the PS
- Unit: [arcsec^2]
- 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(flux, area, freq)
-
- return Tb
-
- def calc_Tb(self, freq):
- """
- Calculate the surface brightness temperature of the point sources.
-
- Parameters
- ------------
- 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,))
- sr_to_arcsec2 = (np.rad2deg(1) * 3600) ** 2 # [sr] -> [arcsec^2]
- # Iteratively calculate Tb
- for i in range(num_ps):
- ps_area = self.ps_catalog['Area (sr)'][i] # [sr]
- area = ps_area * sr_to_arcsec2
- Tb_list[i] = self.calc_single_Tb(area, freq)
-
- return Tb_list
diff --git a/fg21sim/extragalactic/pointsources/starforming.py b/fg21sim/extragalactic/pointsources/starforming.py
deleted file mode 100644
index d4664bb..0000000
--- a/fg21sim/extragalactic/pointsources/starforming.py
+++ /dev/null
@@ -1,258 +0,0 @@
-# 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
-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: float
- Area of the PS
- Unit: [arcsec^2]
- 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(flux, area, freq)
-
- return Tb
-
- def calc_Tb(self, freq):
- """
- Calculate the surface brightness temperature of the point sources.
-
- Parameters
- ------------
- 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,))
- sr_to_arcsec2 = (np.rad2deg(1) * 3600) ** 2 # [sr] -> [arcsec^2]
- # Iteratively calculate Tb
- for i in range(num_ps):
- ps_area = self.ps_catalog['Area (sr)'][i] # [sr]
- area = ps_area * sr_to_arcsec2
- Tb_list[i] = self.calc_single_Tb(area, freq)
-
- return Tb_list
Copyright © 2017-2018 Aaron LI. All Rights Reserved.