Unverified Commit 046e4437 authored by Akke Viitanen's avatar Akke Viitanen
Browse files

implements focc in the plots

parent 5d8957f4

.gitignore

+9 −0
Original line number Diff line number Diff line
@@ -34,6 +34,8 @@ opt/ivano_feature_extraction 
*~

.*.swo

.*.swn

.*.swm

.*.swp



# lincc files

example_*
@@ -47,3 +49,10 @@ _build 


# demo... might add back in later

demo



.venv/

doc/

notebooks/

overleaf/

paper/

table/

Makefile

+16 −0
Original line number Diff line number Diff line
@@ -162,6 +162,9 @@ fig/flux_sigma_eta_ratio_20250725.pdf: src/python/plot_flux_sigma_eta_ratio.py 
fig/flux_sigma_eta_ratio_20250911.pdf: src/python/plot_flux_sigma_eta_ratio.py

	python3 $< $@ 0 | tee table/table_flux_sigma_eta_ratio_20250911.tex



fig/flux_sigma_eta_ratio_20250930.pdf: src/scripts/plots/plot_flux_sigma_eta_ratio.py

	python3 $< $@ 0 | tee table/table_flux_sigma_eta_ratio_20250930.tex



fig/flux_sigma_eta_ratio_u.pdf: python/plot_flux_sigma_eta_ratio.py

	python3 $< $@ 0 u > table/flux_sigma_eta_u.table

fig/flux_sigma_eta_ratio_g.pdf: python/plot_flux_sigma_eta_ratio.py
@@ -254,3 +257,16 @@ tests: 
debug:

	PYTHONPATH=opt:src/python:$(PYTHONPATH) \

		python3 -m pdb src/python/main.py --config etc/config_test.ini





figures_20250930:

	# xlf

	python3 src/scripts/plots/plot_xray_luminosity_function.py data/catalog/dr1_new_new/ fig/xray_luminosity_function_20250930.pdf 0

	python3 src/scripts/plots/plot_xray_luminosity_function.py data/catalog/dr1_new_new/ fig/xray_luminosity_function__focc_20250930.pdf 1

	# qlf

	python3 src/scripts/plots/plot_quasar_luminosity_function.py fig/quasar_luminosity_function_20250930.pdf 0

	python3 src/scripts/plots/plot_quasar_luminosity_function.py fig/quasar_luminosity_function_focc_20250930.pdf 1

	# bhmf

	python3 src/scripts/plots/plot_black_hole_mass_function.py

	# gband number counts

	python3 src/scripts/plots/plot_number_counts_gband.py

pyproject.toml

+13 −0
Original line number Diff line number Diff line
@@ -17,6 +17,19 @@ classifiers = [ 
dynamic = ["version"]

requires-python = ">=3.9"

dependencies = [

    "astropy",

    "fitsio",

    "ipywidgets",

    "matplotlib",

    "numpy",

    "mypy",

    "pandas",

    "scipy",

    "sphinx",

    "sphinx_rtd_theme",

    "sphinx-autoapi",

    "sphinx_copybutton",

    "nbsphinx",

]



[project.urls]

src/lsst_inaf_agile/_version.py

+3 −3
Original line number Diff line number Diff line
@@ -27,7 +27,7 @@ version_tuple: VERSION_TUPLE 
commit_id: COMMIT_ID

__commit_id__: COMMIT_ID



__version__ = version = "0.1.dev27+g8ff1cd42d.d20250912"

__version_tuple__ = version_tuple = (0, 1, "dev27", "g8ff1cd42d.d20250912")

__version__ = version = "0.1.dev42+g5d8957f44.d20250930"

__version_tuple__ = version_tuple = (0, 1, "dev42", "g5d8957f44.d20250930")



__commit_id__ = commit_id = "g8ff1cd42d"
 
__commit_id__ = commit_id = "g5d8957f44"

src/lsst_inaf_agile/ananna2022.py

0 → 100644
+182 −0
Original line number Diff line number Diff line 
#!/usr/bin/env python3

# Author: Akke Viitanen

# Email: akke.viitanen@helsinki.fi

# Date: 2024-08-01 13:53:22



"""

Implements Ananna+2022

"""



import matplotlib.pyplot as plt

import numpy as np



labels = [

    r"Intrinsic ($\sigma=0.3$)",

    r"Intrinsic ($\sigma=0.3; \sigma_{\log L,{\rm scatt}} = 0.2$)",

    r"Intrinsic ($\sigma=0.5$)",

    r"Intrinsic ($\sigma=0.3; {\rm OA} = 35^{\circ}$)",

    r"$1/V_{\rm max}$",

]



###############################################################################

# Ananna+ 2022, Table 3

parameters = {

    "BHMF": {

        "All": [

            (labels[0], 10**7.88, 10**-3.52, -1.576, 0.593),

            (labels[1], 10**7.92, 10**-3.67, -1.530, 0.612),

            (labels[2], 10**7.67, 10**-3.37, -1.260, 0.630),

            (labels[3], 10**7.92, 10**-3.49, -1.576, 0.600),

            (labels[4], 10**8.12, 10**-4.33, -1.060, 0.574),

        ],

        "Type 1": [

            (labels[0], 10**7.97, 10**-4.19, -1.753, 0.561),

            (labels[1], 10**7.93, 10**-4.27, -1.730, 0.566),

            (labels[2], 10**7.91, 10**-4.27, -1.560, 0.590),

            (labels[4], 10**8.73, 10**-5.10, -1.350, 0.681),

        ],

        "Type 2": [

            (labels[0], 10**7.820, 10**-3.60, -1.16, 0.637),

            (labels[1], 10**7.790, 10**-3.64, -1.18, 0.617),

            (labels[2], 10**7.760, 10**-3.60, -0.99, 0.703),

            (labels[3], 10**7.730, 10**-3.44, -1.26, 0.635),

            (labels[4], 10**8.102, 10**-4.33, -1.04, 0.732),

        ],

    },

    "ERDF": {

        "All": [

            (labels[0], 10**-1.338, 10**-3.64, 0.38, 2.260),

            (labels[1], 10**-1.286, 10**-3.76, 0.40, 2.322),

            (labels[2], 10**-1.332, 10**-3.68, 0.484, 2.210),

            (labels[3], 10**-1.249, 10**-3.80, 0.28, 2.720),

            (labels[4], 10**-1.190, 10**-3.76, -0.02, 2.060),

        ],

        "Type 1": [

            (labels[0], 10**-1.152, 10**-4.08, 0.30, 2.51),

            (labels[1], 10**-1.138, 10**-4.09, 0.27, 2.57),

            (labels[2], 10**-1.103, 10**-4.23, 0.13, 2.97),

            (labels[4], 10**-1.060, 10**-4.02, -0.51, 2.57),

        ],

        "Type 2": [

            (labels[0], 10**-1.657, 10**-3.82, 0.376, 2.50),

            (labels[1], 10**-1.628, 10**-3.84, 0.320, 2.50),

            (labels[2], 10**-1.675, 10**-3.80, 0.330, 2.51),

            (labels[3], 10**-1.593, 10**-3.92, 0.300, 2.53),

            (labels[4], 10**-1.870, 10**-3.74, -0.500, 2.30),

        ],

    },

}





def get_phi_bh(m, m_star, phi_star, alpha, beta, h=1.0, sample=None):

    """

    Return the Schechter function form of BHMF

    """



    if sample:

        ## NOTE: errors for sigma=0.50 case

        # m_star = 10 ** (np.log10(m_star) + np.random.uniform(-0.20, 0.25))

        # alpha += np.random.uniform(-0.110, +0.190)

        # beta  += np.random.uniform(-0.086, +0.065)



        m_star = 10 ** (np.log10(m_star) + np.random.normal(scale=0.50 * (0.20 + 0.25)))

        alpha += np.random.normal(scale=0.50 * (0.110 + 0.190))

        beta += np.random.normal(scale=0.50 * (0.086 + 0.065))



    x = m / m_star

    ret = np.log(10) * phi_star * x ** (alpha + 1) * np.exp(-(x**beta))



    # NOTE: fix for h

    #   original unit is 1/(Mpc/h)^3

    #   so that e.g. h = 0.70 corresponds to 1/(Mpc/0.70)^3 = 1/Mpc3 * 0.70^3



    return ret * h**3





def get_phi_lambda(lambda_edd, lambda_edd_star, xi_star, delta1, epsilon_lambda, h=1.0):

    ratio = lambda_edd / lambda_edd_star

    return (

        np.ma.true_divide(xi_star, np.power(ratio, delta1) + np.power(ratio, delta1 + epsilon_lambda)) * h**3

    )





def get_phi_bh_fig10(m, is_type1=True, is_type2=True, log_ledd_gt=-3, h=1.0):

    x = np.log10(m)

    y = np.zeros_like(m)



    x1, y1 = np.loadtxt(f"data/ananna2022/fig10/bhmf/bhmf_ledd_gt_{log_ledd_gt:.1f}_type1.dat").T

    x2, y2 = np.loadtxt(f"data/ananna2022/fig10/bhmf/bhmf_ledd_gt_{log_ledd_gt:.1f}_type2.dat").T



    if is_type1:

        y += 10 ** np.interp(x, x1, y1, left=-np.inf, right=-np.inf)

    if is_type2:

        y += 10 ** np.interp(x, x2, y2, left=-np.inf, right=-np.inf)



    return y * h**3





def get_phi_lambda_fig10(lambda_edd, is_type1=True, is_type2=True, log_mbh_gt=6.5, h=1.0):

    x = np.log10(lambda_edd)

    y = np.zeros_like(lambda_edd)



    x1, y1 = np.loadtxt(f"data/ananna2022/fig10/erdf/erdf_mbh_gt_{log_mbh_gt:.1f}_type1.dat").T

    x2, y2 = np.loadtxt(f"data/ananna2022/fig10/erdf/erdf_mbh_gt_{log_mbh_gt:.1f}_type2.dat").T



    if is_type1:

        y += 10 ** np.interp(x, x1, y1, left=-np.inf, right=-np.inf)

    if is_type2:

        y += 10 ** np.interp(x, x2, y2, left=-np.inf, right=-np.inf)



    return y * h**3





if __name__ == "__main__":

    mbh = np.logspace(6.5, 10, 41)

    lambda_edd = np.logspace(-3.0, 0.5, 36)

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



    for i in range(5):

        for j, k in enumerate(["All", "Type 1", "Type 2"]):

            try:

                p = parameters["BHMF"][k][i]

                phi = get_phi_bh(mbh, *p[1:])

                axes[0, j].loglog(mbh, phi, label=p[0])

                axes[0, j].set_xlim(10**6.5, 10**9.5)

                axes[0, j].set_ylim(1e-10, 1e-2)

                axes[0, j].set_title(k)

                axes[0, j].set_xlabel(r"$M_{\rm BH}$ [Msun]")

                axes[0, j].set_ylabel(r"$\Phi_{\rm BH}$ [1/(Mpc3/h3)/dex")



                p = parameters["ERDF"][k][i]

                phi = get_phi_lambda(lambda_edd, *p[1:])

                axes[1, j].loglog(lambda_edd, phi, label=p[0])

                axes[1, j].set_xlim(10**-3.0, 10**+0.5)

                axes[1, j].set_ylim(10**-8.0, 10**-2.5)

                axes[1, j].set_xlabel(r"$\lambda{\rm Edd}$")

                axes[1, j].set_ylabel(r"$\Phi_{\lambda}$ [1/(Mpc3/h3)/dex")



            except IndexError:

                pass



    for ax in axes.flatten():

        ax.legend()



    plt.show()

    quit()



    plt.figure()

    for p in parameters["BHMF"]["All"]:

        plt.plot(mbh, get_phi_bh(mbh, *p[1:]), label=p[0])

    plt.xlabel(r"$M_{\rm BH}$ [Msun]")

    plt.ylabel(r"$\Phi_{M}$ [1/(Mpc/h)$^3$/dex]")

    plt.legend()

    plt.loglog()



    plt.figure()

    for p in parameters["ERDF"]["All"]:

        plt.plot(lambda_edd, get_phi_lambda(lambda_edd, *p[1:]), label=p[0])

    plt.xlabel(r"$\lambda$")

    plt.ylabel(r"$\Phi_{\lambda}$ [1/(Mpc/h)$^3$/dex]")

    plt.legend()

    plt.loglog()

    plt.show()