#!/usr/bin/env python3
# Author: Akke Viitanen
# Email: akke.viitanen@helsinki.fi
# Date: 2024-08-01 13:53:22
"""Implement Ananna+ 2022."""
import matplotlib.pyplot as plt
import numpy as np
from numpy.typing import ArrayLike
[docs]
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
[docs]
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),
],
},
}
[docs]
def get_phi_bh(
m: ArrayLike,
m_star: float,
phi_star: float,
alpha: float,
beta: float,
h: float = 1.0,
sample: bool = False,
) -> ArrayLike:
"""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 = np.asarray(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
[docs]
def get_phi_lambda(
lambda_edd: ArrayLike,
lambda_edd_star: float,
xi_star: float,
delta1: float,
epsilon_lambda: float,
h: float = 1.0,
) -> ArrayLike:
"""Return phi_lambda following the functional form in Ananna."""
ratio = np.asarray(lambda_edd) / lambda_edd_star
return (
np.ma.true_divide(xi_star, np.power(ratio, delta1) + np.power(ratio, delta1 + epsilon_lambda)) * h**3
)
[docs]
def get_phi_bh_fig10(
m: ArrayLike, is_type1: bool = True, is_type2: bool = True, log_ledd_gt: float = -3, h: float = 1.0
):
"""Get phi_bh from Ananna Fig10."""
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
[docs]
def get_phi_lambda_fig10(
lambda_edd: ArrayLike,
is_type1: bool = True,
is_type2: bool = True,
log_mbh_gt: float = 6.5,
h: float = 1.0,
) -> ArrayLike:
"""Get phi_lambda from Ananna Fig10."""
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__":
[docs]
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()