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, 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