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.

1 files changed, 238 insertions, 0 deletions

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 <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: `~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