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.

10 files changed, 1834 insertions, 0 deletions

diff --git a/fg21sim/extragalactic/__init__.py b/fg21sim/extragalactic/__init__.py
index af976ca..f57928f 100644
--- a/fg21sim/extragalactic/__init__.py
+++ b/fg21sim/extragalactic/__init__.py
@@ -2,3 +2,8 @@
 # MIT license
 
 from .clusters import GalaxyClusters
+
+# Copyright (c) 2016 Zhixian MA <zxma_sjtu@qq.com>
+# MIT license
+
+from .pointsources import PointSources
diff --git a/fg21sim/extragalactic/pointsources/__init__.py b/fg21sim/extragalactic/pointsources/__init__.py
new file mode 100644
index 0000000..4c64c98
--- /dev/null
+++ b/fg21sim/extragalactic/pointsources/__init__.py
@@ -0,0 +1,4 @@
+# Copyright (c) 2016 Zhixian MA <zxma_sjtu@qq.com>

+# MIT license

+

+from .pointsources import PointSources

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
diff --git a/fg21sim/extragalactic/pointsources/fr1.py b/fg21sim/extragalactic/pointsources/fr1.py
new file mode 100644
index 0000000..0e052cd
--- /dev/null
+++ b/fg21sim/extragalactic/pointsources/fr1.py
@@ -0,0 +1,355 @@
+# 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: `~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:`~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)
+        core_area = 4 * np.pi / npix  # [sr]
+        Tb_core = convert.Fnu_to_Tb_fast(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_fast(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
+        ------------
+        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, 2))
+        # 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
diff --git a/fg21sim/extragalactic/pointsources/fr2.py b/fg21sim/extragalactic/pointsources/fr2.py
new file mode 100644
index 0000000..b40f9ff
--- /dev/null
+++ b/fg21sim/extragalactic/pointsources/fr2.py
@@ -0,0 +1,380 @@
+# 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: `~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:`~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)
+        core_area = 4 * np.pi / npix  # [sr]
+        Tb_core = convert.Fnu_to_Tb_fast(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_fast(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
+        ------------
+        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, 3))
+        # 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
diff --git a/fg21sim/extragalactic/pointsources/pointsources.py b/fg21sim/extragalactic/pointsources/pointsources.py
new file mode 100644
index 0000000..85f6ec4
--- /dev/null
+++ b/fg21sim/extragalactic/pointsources/pointsources.py
@@ -0,0 +1,126 @@
+# 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!")
+
+        return hpmap_f
+
+    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
new file mode 100644
index 0000000..7c3a93a
--- /dev/null
+++ b/fg21sim/extragalactic/pointsources/psparams.py
@@ -0,0 +1,77 @@
+# 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 (\deg)
+    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
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
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
diff --git a/fg21sim/extragalactic/pointsources/starforming.py b/fg21sim/extragalactic/pointsources/starforming.py
new file mode 100644
index 0000000..5d518d8
--- /dev/null
+++ b/fg21sim/extragalactic/pointsources/starforming.py
@@ -0,0 +1,257 @@
+# 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: `~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