Unverified Commit e93d99c3 authored by Akke Viitanen's avatar Akke Viitanen
Browse files

add schulze2015 docstrings

parent 154491c1

src/lsst_inaf_agile/schulze2015.py

+41 −22
Original line number Diff line number Diff line
@@ -20,6 +20,35 @@ def get_log_Psi( 
    log_M_bh_c=8,

    redshift_c=1.6,

):

    """

    Return the AGN bivariate (MBH-lambda) distribution function.



    The bivariate distribution function is the number density of AGN within

    some tiny dlogMBH, dloglambda interval. The distribution has the dimensions

    of 1/volume/dbin ** 2 (i.e. typically 1/Mpc3/dex2). One-dimensional

    distribution functions may be derived by integrating along a dimension of

    the bivariate distribution function.



    See: https://ui.adsabs.harvard.edu/abs/2015MNRAS.447.2085S/abstract



    Parameters

    ----------

    log_M_bh: float

        log10 of the SMBH mass in Msun.

    log_lambda_edd: float

        log10 of the Eddington ratio (dimensionless).

    redshift: float

        Redshift of the bivariate distribution function.

    *parameters: float

        The 13 parameters that fully define the Schulze+2015 model. Refer

        Schulze+2015 Sec. 4.2.



    Returns

    -------

    log_Psi: float

        The bivariate distribution function corresponding to the parameters.



    """

    # Calculate log_rho_bh

    log_M_bh_star = log_M_bh_star + c_bh * (redshift - redshift_c)

    log_x = log_M_bh - log_M_bh_star
@@ -43,34 +72,24 @@ def get_log_Psi( 




@np.vectorize

def get_log_Phi_bh(log_M_bh, redshift, log_lambda_edd_min=-2, log_lambda_edd_max=1):

def get_log_Phi_bh(log_M_bh, redshift, log_lambda_edd_min=-2, log_lambda_edd_max=1, *args, **kwargs):

    """

    Return 1-d BHMF by integrating along lambda_edd.

    """



    def fun(log_lambda_edd):

        return 10 ** get_log_Psi(log_M_bh, log_lambda_edd, redshift)

        return 10 ** get_log_Psi(log_M_bh, log_lambda_edd, redshift, *args, **kwargs)



    return np.log10(quad(fun, log_lambda_edd_min, log_lambda_edd_max)[0])





@np.vectorize

def get_log_Phi_lambda_edd(log_lambda_edd, redshift, log_M_bh_min=7, log_M_bh_max=11):

def get_log_Phi_lambda_edd(log_lambda_edd, redshift, log_M_bh_min=7, log_M_bh_max=11, *args, **kwargs):

    """

    Return 1-d ERDF by integrating along MBH.

    """



    def fun(log_M_bh):

        return 10 ** get_log_Psi(log_M_bh, log_lambda_edd, redshift)

        return 10 ** get_log_Psi(log_M_bh, log_lambda_edd, redshift, *args, **kwargs)



    return np.log10(quad(fun, log_M_bh_min, log_M_bh_max)[0])
 




if __name__ == "__main__":

    log_m_bh, log_lambda_edd = np.meshgrid(np.linspace(7, 11), np.linspace(-2, 1))

    redshifts = 1.3, 1.9



    import matplotlib as mpl

    import matplotlib.pyplot as plt



    fig, axes = plt.subplots(1, 2, figsize=(2 * 6.4, 4.8))

    cmap = mpl.cm.get_cmap("viridis")

    norm = mpl.colors.BoundaryNorm([-11, -10, -9, -8, -7, -6, -5, -4], cmap.N)

    for i, redshift in enumerate(redshifts):

        log_Psi = get_log_Psi(log_m_bh, log_lambda_edd, redshift)

        im = axes[i].imshow(log_Psi, cmap=cmap, norm=norm, extent=(7, 11, -2, 1))

        plt.colorbar(im, ax=axes[i])

        axes[i].set_title("z = %.2f" % redshift)

    plt.savefig("fig/schulze2015_fig10.pdf")

tests/lsst_inaf_agile/test_schulze2015.py

+29 −0
Original line number Diff line number Diff line 
from unittest import TestCase



import numpy as np

from lsst_inaf_agile import schulze2015
@@ -90,3 +92,30 @@ def test_get_log_Phi_lambda_edd(num_regression): 
        for i, args in enumerate(args_all)

    }

    num_regression.check(to_check)





class Schulze2015(TestCase):

    def test_plot_schulze2015(self):

        import matplotlib as mpl

        import matplotlib.pyplot as plt

        from lsst_inaf_agile import util



        log_m_bh, log_lambda_edd = np.meshgrid(np.linspace(7, 11), np.linspace(-2, 1))

        redshifts = 1.3, 1.9



        fig, axes = plt.subplots(1, 2, figsize=(2 * 6.4, 4.8), dpi=300)

        cmap = plt.get_cmap("Oranges")

        norm = mpl.colors.BoundaryNorm([-11, -10, -9, -8, -7, -6, -5, -4], cmap.N)

        for i, redshift in enumerate(redshifts):

            log_Psi = schulze2015.get_log_Psi(log_m_bh, log_lambda_edd, redshift)

            im = axes[i].imshow(log_Psi, cmap=cmap, norm=norm, extent=(7, 11, -2, 1))

            plt.colorbar(im, ax=axes[i])

            axes[i].set_title("z = %.2f" % redshift)



        axes[0].set_xlabel(r"$\log \left( M_\mathrm{BH} / M_\odot \right)$")

        axes[1].set_xlabel(r"$\log \left( M_\mathrm{BH} / M_\odot \right)$")

        axes[0].set_ylabel(r"$\log \left( \lambda_\mathrm{Edd} \right)$")

        fig.tight_layout()



        util.create_directory("data/tests/test_schulze2015/")

        plt.savefig("data/tests/test_schulze2015/fig10.pdf", bbox_inches="tight")