#!/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
[docs]
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,
)
[docs]
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],
)
[docs]
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
[docs]
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