#!/usr/bin/env python3
# Author: Akke Viitanen
# Email: akke.viitanen@helsinki.fi
# Date: 2023-06-26 10:51:39
"""Implement Ueda+2014."""
import glob
import re
import matplotlib.pyplot as plt
import numpy as np
[docs]
def get_Psi(log_L_X, z, beta=0.24, Psi_min=0.20, Psi_max=0.84, *args, **kwargs):
"""Eq. 3."""
log_L_X = np.atleast_1d(log_L_X)
z = np.atleast_1d(z)
ret = []
Psi_4375 = get_Psi_4375(z, *args, **kwargs)
ret = Psi_4375 - beta * (log_L_X - 43.75)
ret[ret < Psi_min] = Psi_min
ret[ret > Psi_max] = Psi_max
return ret
[docs]
def get_Psi_4375(z, Psi_4375_0=0.43, a1=0.48):
"""Eq. 4."""
return np.where(
z < 2.0,
Psi_4375_0 * (1 + z) ** a1,
Psi_4375_0 * (1 + 2) ** a1,
)
[docs]
def get_f(log_L_X, z, log_N_H, epsilon=1.7, f_CTK=1.0, *args, **kwargs):
"""Eqs. 5 & 6."""
log_L_X = np.atleast_2d(log_L_X)
z = np.atleast_2d(z)
log_N_H = np.atleast_2d(log_N_H)
Psi = get_Psi(log_L_X, z)
select_Psi_1 = Psi < (1 + epsilon) / (3 + epsilon)
select_Psi_2 = Psi >= (1 + epsilon) / (3 + epsilon)
select_NH_00_21 = (log_N_H >= 0) * (log_N_H < 21)
select_NH_21_22 = (log_N_H >= 21) * (log_N_H < 22)
select_NH_22_23 = (log_N_H >= 22) * (log_N_H < 23)
select_NH_23_24 = (log_N_H >= 23) * (log_N_H < 24)
select_NH_24_99 = (log_N_H >= 24) * (log_N_H < 99)
ret = np.zeros_like(Psi)
# NOTE: the following contains all the conditions and their piecewise
# definitions from Ueda+2014
select_value = [
# Case Psi < (1 + epsilon) / (3 + epsilon)
(select_Psi_1 * select_NH_00_21, 1 - (2 + epsilon) / (1 + epsilon) * Psi),
(select_Psi_1 * select_NH_21_22, 1 / (1 + epsilon) * Psi),
(select_Psi_1 * select_NH_22_23, 1 / (1 + epsilon) * Psi),
(select_Psi_1 * select_NH_23_24, epsilon / (1 + epsilon) * Psi),
(select_Psi_1 * select_NH_24_99, f_CTK / 2 * Psi),
# Case Psi >= (1 + epsilon) / (3 + epsilon)
(select_Psi_2 * select_NH_00_21, 2 / 3.0 - (3 + 2 * epsilon) / (3 + 3 * epsilon) * Psi),
(select_Psi_2 * select_NH_21_22, 1 / 3.0 - epsilon / (3 + 3 * epsilon) * Psi),
(select_Psi_2 * select_NH_22_23, 1 / (1 + epsilon) * Psi),
(select_Psi_2 * select_NH_23_24, epsilon / (1 + epsilon) * Psi),
(select_Psi_2 * select_NH_24_99, f_CTK / 2 * Psi),
]
for select, value in select_value:
ret[select] = value[select]
return np.squeeze(ret)
# def get_single(log_L_X, z):
# Psi = get_Psi(log_L_X, z)
# if Psi < (1 + epsilon) / (3 + epsilon):
# if log_N_H < 21: return 1 - (2 + epsilon) / (1 + epsilon) * Psi
# if 21 <= log_N_H < 22: return 1 / (1 + epsilon) * Psi
# if 22 <= log_N_H < 23: return 1 / (1 + epsilon) * Psi
# if 23 <= log_N_H < 24: return epsilon / (1 + epsilon) * Psi
# if 24 <= log_N_H : return f_CTK / 2 * Psi
# if Psi >= (1 + epsilon) / (3 + epsilon):
# if log_N_H < 21: return 2 / 3. - (3 + 2 * epsilon) / (3 + 3 * epsilon) * Psi
# if 21 <= log_N_H < 22: return 1 / 3. - epsilon / (3 + 3 * epsilon) * Psi
# if 22 <= log_N_H < 23: return 1 / (1 + epsilon) * Psi
# if 23 <= log_N_H < 24: return epsilon / (1 + epsilon) * Psi
# if 24 <= log_N_H : return f_CTK / 2 * Psi
# return np.vectorize(get_single)(log_L_X, z)
[docs]
def get_fraction_agn_absorbed(log_L_X, z):
"""Return Fig. 5 from Ueda+2014 i.e. fraction of NH=22-24 sources to NH=20=24."""
f = [get_f(log_L_X, z, n) for n in [20, 21, 22, 23]]
return np.sum(f[-2:], axis=0) / np.sum(f, axis=0)
[docs]
def get_xray_luminosity_function_model(
log_lx,
z,
A=2.91 * 10**-6,
log_L_star=43.97,
gamma_1=0.96,
gamma_2=2.71,
p1=4.78,
p2=-1.5,
p3=-6.2,
log_L_p=44,
beta_1=0.84,
z_c1=1.86,
z_c2=3.00,
log_L_a1=44.61,
log_L_a2=44.00,
alpha_1=0.29,
alpha_2=-0.10,
):
"""
Estimate Ueda+2014 XLF model.
This function implements Ueda+2014 verbatim. Refer to Ueda+2014 and the
relevant equations for the details.
"""
# Estimate p1
p1 = p1 + beta_1 * (log_lx - log_L_p)
# Estimate the critical redshifts
z_c1 = np.where(log_lx <= log_L_a1, z_c1 * (10 ** (log_lx - log_L_a1)) ** alpha_1, z_c1)
z_c2 = np.where(log_lx <= log_L_a2, z_c2 * (10 ** (log_lx - log_L_a2)) ** alpha_2, z_c2)
# Estimate the evolution parameter
e = np.zeros_like(log_lx)
z1 = 1 + z
z1c1 = 1 + z_c1
z1c2 = 1 + z_c2
for s, v in [
((z >= 0.00) * (z < z_c1), z1**p1),
((z_c1 <= z) * (z < z_c2), z1c1**p1 * (z1 / z1c1) ** p2),
((z_c2 <= z) * (z < np.inf), z1c1**p1 * (z1c2 / z1c1) ** p2 * (z1 / z1c2) ** p3),
]:
e[s] = v[s]
# Estimate the local luminosity function
ratio = 10 ** (log_lx - log_L_star)
phi = A * (ratio**gamma_1 + ratio**gamma_2) ** -1
return phi * e
[docs]
def get_xray_luminosity_function_data(z, return_error=False):
"""
Get Ueda+2014 X-ray luminosity function data.
Returns
-------
x, y:
tuple of (logLX, phi) in typical units
"""
filename_ueda = sorted(glob.glob("opt/quasarlf/data/ueda2014/hard*.dat"))
zmin = np.array([float(re.findall("z(.*)-(.*).dat", f)[0][0]) for f in filename_ueda])
zmax = np.array([float(re.findall("z(.*)-(.*).dat", f)[0][1]) for f in filename_ueda])
z_ueda = (zmin * zmax) ** 0.5
idx = np.argmin(np.abs(z - z_ueda))
x, y, y_lo, y_hi = np.loadtxt(filename_ueda[idx]).T
if return_error:
return x, y, y_lo, y_hi
return x, y
[docs]
def main():
"""Produce Ueda+2014 debugging plots."""
z, lx = np.meshgrid(
np.linspace(0.0, 5.5, 56),
np.linspace(42, 46, 41),
)
def get_f_ctk(lx, z):
nhs = 20, 21, 22, 23, 24, 25
f_all = [get_f(lx, z, n) for n in nhs]
f_ctk = (f_all[-2] + f_all[-1]) / np.sum(f_all, axis=0)
return f_ctk
fig = plt.figure()
ax = fig.gca()
ax.set_title("frac_CTK = N_CTK / N_AGN")
ax.imshow(get_f_ctk(lx, z), extent=(0, 5.5, 42, 46))
plt.colorbar(ax=ax)
ax.set_xlabel(r"$z$")
ax.set_ylabel(r"$\log L_{\rm X}$ [erg/s]")
fig.savefig("fig/frac_ctk_ueda2014.pdf")
return 0
if __name__ == "__main__":
main()