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, 203 insertions, 0 deletions

diff --git a/fg21sim/extragalactic/pointsources/base.py b/fg21sim/extragalactic/pointsources/base.py
new file mode 100644
index 0000000..1544f5b
--- /dev/null
+++ b/fg21sim/extragalactic/pointsources/base.py
@@ -0,0 +1,203 @@
+# 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