import numpy as np
from scipy.integrate import quad
[docs]
def get_log_Psi(
log_M_bh,
log_lambda_edd,
redshift,
log_Psi_star=-5.32,
log_M_bh_star=9.09,
alpha=-1.50,
beta=0.96,
log_lambda_edd_star0=-1.19,
alpha_lambda_edd=-0.29,
k_lambda_edd=0.099,
gamma=0.05,
c_bh=0.270,
c_lambda_edd=0.100,
c_alpha_lambda_edd=0.094,
log_M_bh_c=8,
redshift_c=1.6,
):
"""
Return the AGN bivariate (MBH-lambda) distribution function.
The bivariate distribution function is the number density of AGN within
some tiny dlogMBH, dloglambda interval. The distribution has the dimensions
of 1/volume/dbin ** 2 (i.e. typically 1/Mpc3/dex2). One-dimensional
distribution functions may be derived by integrating along a dimension of
the bivariate distribution function.
See: https://ui.adsabs.harvard.edu/abs/2015MNRAS.447.2085S/abstract
Parameters
----------
log_M_bh: float
log10 of the SMBH mass in Msun.
log_lambda_edd: float
log10 of the Eddington ratio (dimensionless).
redshift: float
Redshift of the bivariate distribution function.
*parameters: float
The 13 parameters that fully define the Schulze+2015 model. Refer
Schulze+2015 Sec. 4.2.
Returns
-------
log_Psi: float
The bivariate distribution function corresponding to the parameters.
"""
# Calculate log_rho_bh
log_M_bh_star = log_M_bh_star + c_bh * (redshift - redshift_c)
log_x = log_M_bh - log_M_bh_star
log_rho_bh = (alpha + 1) * log_x - 10 ** (beta * log_x) / np.log(10) + np.log10(np.log(10))
# Calculate log_rho_lambda_edd
log_lambda_edd_star = (
log_lambda_edd_star0 + k_lambda_edd * (log_M_bh - log_M_bh_c) + c_lambda_edd * (redshift - redshift_c)
)
# NOTE: typo in schulze? i.e. log_alpha_l should be alpha_l
# log_alpha_l = alpha_l + c_alpha_l * (z - z_c)
alpha_lambda_edd = alpha_lambda_edd + c_alpha_lambda_edd * (redshift - redshift_c)
log_x = log_lambda_edd - log_lambda_edd_star
log_rho_lambda_edd = (alpha_lambda_edd + 1) * log_x - 10**log_x / np.log(10) + np.log10(np.log(10))
# Calculate log_rho_z
log_rho_redshift = gamma * (redshift - redshift_c)
return log_Psi_star + log_rho_bh + log_rho_lambda_edd + log_rho_redshift
@np.vectorize
[docs]
def get_log_Phi_bh(log_M_bh, redshift, log_lambda_edd_min=-2, log_lambda_edd_max=1, *args, **kwargs):
"""
Return 1-d BHMF by integrating along lambda_edd.
"""
def fun(log_lambda_edd):
return 10 ** get_log_Psi(log_M_bh, log_lambda_edd, redshift, *args, **kwargs)
return np.log10(quad(fun, log_lambda_edd_min, log_lambda_edd_max)[0])
@np.vectorize
[docs]
def get_log_Phi_lambda_edd(log_lambda_edd, redshift, log_M_bh_min=7, log_M_bh_max=11, *args, **kwargs):
"""
Return 1-d ERDF by integrating along MBH.
"""
def fun(log_M_bh):
return 10 ** get_log_Psi(log_M_bh, log_lambda_edd, redshift, *args, **kwargs)
return np.log10(quad(fun, log_M_bh_min, log_M_bh_max)[0])