add Kelly+2013 module

#!/usr/bin/env python3

# Author: Akke Viitanen

# Email: akke.viitanen@helsinki.fi

# Date: 2026-01-08 16:52:11

"""Implement Kelly+2013 BHMF."""

import os

import matplotlib.pyplot as plt

import numpy as np

from astropy.table import Table

REDSHIFTS = (

    0.40,

    0.60,

    0.80,

    1.00,

    1.20,

    1.40,

    1.60,

    1.80,

    2.15,

    2.65,

    3.20,

    3.75,

    4.25,

    4.75,

)

def get_bhmf(z: float, filename: str = "data/kelly2013/table1.dat") -> tuple[float, float, float, float]:

    """

    Return Kelly & Shen 2013 BHMF at the given redshift z.

    The BHMF is based on Table 1 of Kelly & Shen 2013. See 'See Also' below.

    The closest redshift bin to the input redshift is returned without any

    interpolation.

    Parameters

    ----------

    z: float

        Redshift of the black hole mass function.

    filename: str

        Path to the machine-readable Kelly & Shen 2013 table 1.

    Returns

    -------

    log_mbh, log_phi, log_phi_lo log_phi_hi: tuple[float, float, float, float]

        A 4-tuple containing the black hole mass function measurement.

    Raises

    ------

    ValueError

        If redshift is not within the Kelly & Shen 2013 range (z!=0-5).

    See Also

    -----

    https://ui.adsabs.harvard.edu/abs/2013ApJ...764...45K/abstract

    https://iopscience.iop.org/article/10.1088/0004-637X/764/1/45

    """

    if not 0.0 <= z <= 5.0:

        raise ValueError(f"Redshift {z} is not \\in [0.0, 5.0].")

    if not os.path.exists(filename):

        raise FileNotFoundError(f"Filename {filename} not found.")

    table = Table.read(filename, format="ascii")

    zkelly = table["zbar"][np.argmin((table["zbar"] - z) ** 2)]

    select = table["zbar"] == zkelly

    return (

        table["logMBH"][select],

        table["LologPhi"][select],

        table["logPhi"][select],

        table["UplogPhi"][select],

    )

def _get_figax():

    """Return Fig. 4 figure and axes."""

    fig, axes = plt.subplots(4, 4, sharex=True, sharey=True)

    fig.supxlabel(r"$\log \left( M_\mathrm{BH} / M_\odot \right)$")

    fig.supylabel(r"$\log \left[ \phi(M_\mathrm{BH} | z) / \mathrm{Mpc}^{-3}\,\mathrm{dex}^{-1} \right]$")

    axes[-1, -1].remove()

    axes[-1, -2].remove()

    # Set the redshift text

    for redshift, ax in zip(REDSHIFTS, axes.flatten(), strict=False):

        ax.text(0.9, 0.9, f"$z={redshift}$", ha="right", va="top", transform=ax.transAxes)

    return fig, axes

def _plot_bhmf(fig, axes, redshifts):

    """Plot Fig. 4 style BHMF given axes and redshifts."""

    for redshift, ax in zip(redshifts, axes.flatten(), strict=False):

        x, ylo, y, yhi = get_bhmf(redshift)

        ax.plot(x, y, label=f"$z={redshift}$")

        ax.fill_between(x, ylo, yhi, alpha=0.20)

    return fig, axes

#!/usr/bin/env python3

# Author: Akke Viitanen

# Email: akke.viitanen@helsinki.fi

# Date: 2026-01-08 16:52:11

"""Test Kelly+2013 BHMF."""

from unittest import TestCase

import numpy as np

from lsst_inaf_agile import kelly2013, util

class TestKelly2013(TestCase):

    def test_get_bhmf_one(self):

        t = np.asarray(kelly2013.get_bhmf(4.250)).T

        self.assertTrue(np.allclose(t[0], [8.000, -8.81, -7.93, -7.16]))

        self.assertTrue(np.allclose(t[-1], [11.00, -11.6, -10.9, -10.5]))

    def test_get_bhmf_close(self):

        for redshift in 0.39, 0.40, 0.41:

            t = np.asarray(kelly2013.get_bhmf(redshift)).T

            self.assertTrue(np.allclose(t[0, :], [8.000, -5.00, -4.64, -4.21]))

    def test_get_bhmf_invalid_redshift(self):

        with self.assertRaises(ValueError):

            kelly2013.get_bhmf(-1.0)

            kelly2013.get_bhmf(-0.1)

            kelly2013.get_bhmf(+5.5)

    def test_get_bhmf_nonexisting_filename(self):

        with self.assertRaises(FileNotFoundError):

            kelly2013.get_bhmf(0.0, "nonexistingfilename.tmp")

    def test_create_fig4(self):

        fig, axes = kelly2013._get_figax()

        kelly2013._plot_bhmf(fig, axes, kelly2013.REDSHIFTS)

        fig.tight_layout()

        for ax in axes.flatten()[:-2]:

            ax.set_xlim(8.0, 11.0)

            ax.set_ylim(-10.5, -3.0)

        savefig = "data/tests/kelly2013/fig4.pdf"

        util.create_directory(savefig)

        fig.savefig(savefig, bbox_inches="tight")