diff --git a/agnpy/absorption/absorption.py b/agnpy/absorption/absorption.py index f66958b..0c6e74e 100644 --- a/agnpy/absorption/absorption.py +++ b/agnpy/absorption/absorption.py @@ -23,10 +23,12 @@ x_re_shell_mu_s, ) from ..utils.conversion import nu_to_epsilon_prime, to_R_g_units -from ..targets import PointSourceBehindJet, SSDisk, SphericalShellBLR, RingDustTorus +from ..targets import PointSourceBehindJet, SSDisk, SphericalShellBLR, RingDustTorus, lines_dictionary from ..emission_regions import Blob from ..synchrotron import nu_synch_peak, Synchrotron +import cubepy as cp + __all__ = ["sigma", "Absorption", "ebl_files_dict", "EBL"] @@ -443,6 +445,111 @@ def evaluate_tau_blr_mu_s( ) return (prefactor * integral).to_value("") + @staticmethod + def compute_integrand_for_blr(cubepy_params_array, R_line, mu_s, epsilon_1, epsilon_line): + """ + Computes the integrand for the gamma-gamma absorption model with + a spherical shell Broad Line Region (BLR). + + Parameters + ---------- + cubepy_params_array : array-like + Array containing integration parameters [mu, phi, log_l]. + R_line : :class:`~astropy.units.Quantity` + Radius of the BLR spherical shell. + mu_s : float + Cosine of the angle between the blob motion and the jet axis. + epsilon_1 : float + Dimensionless energy parameter for the incoming photon. + epsilon_line : float + Dimensionless energy of the emitted line. + + Returns + ------- + float + The value of the integrand for the given parameters in cm^2. + """ + (mu, phi, log_l) = cubepy_params_array + l = np.exp(log_l) + mu_star = mu_star_shell(mu, R_line.to_value('cm'), l) + _cos_psi = cos_psi(mu_s, mu_star, phi) + x = x_re_shell(mu, R_line.to_value('cm'), l) + s = epsilon_1 * epsilon_line * (1 - _cos_psi) / 2 + integrand = (1 - _cos_psi) * l / x ** 2 * sigma(s).to_value("cm**2") + return integrand + @staticmethod + def evaluate_tau_blr_cubepy( + nu, + eps_abs, + z, + mu_s, + L_disk, + xi_line, + epsilon_line, + R_line, + r + ): + """Evaluates the gamma-gamma absorption produced by a spherical shell + BLR for a general set of model parameters using cubepy + + Parameters + ---------- + nu : :class:`~astropy.units.Quantity` + array of frequencies, in Hz, to compute the tau + **note** these are observed frequencies (observer frame) + z : float + redshift of the source + mu_s : float + cosine of the angle between the blob motion and the jet axis + L_disk : :class:`~astropy.units.Quantity` + Luminosity of the disk whose radiation is being reprocessed by the BLR + xi_line : float + fraction of the disk radiation reprocessed by the BLR + epsilon_line : string + dimensionless energy of the emitted line + R_line : :class:`~astropy.units.Quantity` + radius of the BLR spherical shell + r : :class:`~astropy.units.Quantity` + distance between the Broad Line Region and the blob + mu, phi : :class:`~numpy.ndarray` + arrays of cosine of zenith and azimuth angles to integrate over + + Returns + ------- + :class:float + array of the tau values corresponding to each frequency + """ + # conversions + epsilon_1 = nu_to_epsilon_prime(nu, z) + # set boundary for integration + low = np.array( + [ + [min(mu_to_integrate)], + [min(phi_to_integrate)], + [np.log(r.to_value('cm'))] + ] + ) + high = np.array( + [ + [max(mu_to_integrate)], + [max(phi_to_integrate)], + [np.log(1e5*r.to_value('cm'))] + ] + ) + prefactor = (L_disk * xi_line) / ( + (4 * np.pi) ** 2 * epsilon_line * m_e * c ** 3 + ) + + prefactor = prefactor.to("cm-1") + eps = eps_abs / prefactor.to_value("cm-1") + value, error = cp.integrate(Absorption.compute_integrand_for_blr, + low, + high, + abstol=eps, + args=(R_line, mu_s, epsilon_1, epsilon_line)) + integral = value[0]*u.cm + return (prefactor * integral[0]).to_value("") + def tau_blr(self, nu): """Evaluates the gamma-gamma absorption produced by a spherical shell BLR for a general set of model parameters @@ -479,6 +586,116 @@ def tau_blr_mu_s(self, nu): phi=self.phi, ) + def tau_blr_cubepy(self, nu, eps_abs=1e-6): + """Evaluates the gamma-gamma absorption produced by a spherical shell + BLR for a general set of model parameters with cubepy integration + method + """ + + tau_blr_list = [] + for freq in nu: + tau_blr_list.append( + self.evaluate_tau_blr_cubepy( + freq, + eps_abs, + self.z, + self.mu_s, + self.target.L_disk, + self.target.xi_line, + self.target.epsilon_line, + self.target.R_line, + self.r + ) + ) + return tau_blr_list + + def tau_blr_all_lines_cubepy(self, nu, eps_abs=1e-6): + """ + Calculate the absorption in all available BLR (Broad Line Region) shells. + This function scales the radius and luminosity to the Hbeta line + and applies the computation to all the available BLR lines (shells) using + their respective ratios. + + Parameters + ---------- + nu: numpy array + Frequency array with unit + eps_abs: float, optional + Absolute precision for the calculation, default is 1e-6 + """ + SUM_RATIO_LINES_BLR = 30.652 + XI_TARGET_LINE = self.target.xi_line + tau_z_all = np.zeros(nu.shape) + + blr_shells = dict() + absorption_shells = dict() + + # TODO write test + if self.target.line != "Hbeta": + raise KeyError("Please use Hbeta as reference") + + for line_key in lines_dictionary.keys(): + xi_line_current = lines_dictionary[line_key]['L_Hbeta_ratio'] + xi_line = (xi_line_current * XI_TARGET_LINE) / SUM_RATIO_LINES_BLR + R_line_current = self.target.R_line * lines_dictionary[line_key]['R_Hbeta_ratio'] + + blr_shells[line_key] = SphericalShellBLR(self.target.L_disk, xi_line, line_key, R_line_current) + absorption_shells[line_key] = Absorption(blr_shells[line_key], r=self.r, z=self.z) + + tau_z_current = absorption_shells[line_key].tau_blr_cubepy(nu, eps_abs=eps_abs) + tau_z_all += tau_z_current + + return tau_z_all + + def tau_blr_lines_cubepy(self, nu, lines_list, eps_abs=1e-6): + """ + Calculate the absorption in all available BLR (Broad Line Region) shells. + This function scales the radius and luminosity to the Hβ (H-beta) line + and computes the absorption for all available BLR lines (shells) using + their respective ratios. + + Parameters + ---------- + nu : numpy.ndarray + Frequency array for which the absorption will be calculated. + lines_list : list of str + List of BLR emission line names (e.g., 'Hbeta', 'MgII') for which the + absorption will be calculated. Each name corresponds to a predefined + line with associated properties (e.g., wavelength, radius ratio, and + luminosity ratio) that are stored as lines_dictionary. + eps_abs : float, optional + Absolute precision for the calculation (default is 1e-6). + + Returns + ------- + tau : numpy.ndarray + Optical depth values corresponding to the input frequencies. + """ + + sum_ratio_lines = sum(lines_dictionary[line]['L_Hbeta_ratio'] for line in lines_list) + XI_TARGET_LINE = self.target.xi_line + tau_z_lines = np.zeros(nu.shape) + + blr_shells = dict() + absorption_shells = dict() + + # TODO write test + if self.target.line != "Hbeta": + raise KeyError("Please use Hbeta as reference") + + for line_key in lines_list: + xi_line_current = lines_dictionary[line_key]['L_Hbeta_ratio'] + xi_line = (xi_line_current * XI_TARGET_LINE) / sum_ratio_lines + R_line_current = self.target.R_line * lines_dictionary[line_key]['R_Hbeta_ratio'] + + blr_shells[line_key] = SphericalShellBLR(self.target.L_disk, xi_line, line_key, R_line_current) + absorption_shells[line_key] = Absorption(blr_shells[line_key], r=self.r, z=self.z) + + tau_z_current = absorption_shells[line_key].tau_blr_cubepy(nu, eps_abs=eps_abs) + tau_z_lines += tau_z_current + + return tau_z_lines + @staticmethod def evaluate_tau_dt( nu, diff --git a/agnpy/absorption/tests/test_absorption.py b/agnpy/absorption/tests/test_absorption.py index 817191d..649a058 100644 --- a/agnpy/absorption/tests/test_absorption.py +++ b/agnpy/absorption/tests/test_absorption.py @@ -8,7 +8,7 @@ from agnpy.spectra import PowerLaw, BrokenPowerLaw from agnpy.emission_regions import Blob from agnpy.synchrotron import Synchrotron -from agnpy.targets import PointSourceBehindJet, SSDisk, SphericalShellBLR, RingDustTorus +from agnpy.targets import PointSourceBehindJet, SSDisk, SphericalShellBLR, RingDustTorus, lines_dictionary from agnpy.absorption import Absorption, EBL from agnpy.utils.math import axes_reshaper from agnpy.utils.validation_utils import ( @@ -228,6 +228,31 @@ def test_abs_dt_vs_point_source(self): # requires a 10% deviation from the two SED points assert check_deviation(nu, tau_dt, tau_ps_dt, 0.1) + def test_abs_blr_cubepy(self): + """Checks for any errors that may occur during the execution of the function that + calculates absorption in BLR with cubepy.""" + SUM_RATIO_LINES_BLR = 30.652 + E = np.logspace(-1, 1, 20) * u.TeV + nu = E.to("Hz", equivalencies=u.spectral()) + L_disk = 6e45 * u.Unit("erg s-1") + R_line = 2.45 * 1e17 * u.cm + r_blob_bh = 1e18 * u.cm + z = 1 + XI_LINE = 0.1 + blr_dict = dict() + ec_blr_dict = dict() + tau_z_all = np.zeros(nu.shape) + for line in list(lines_dictionary.keys()): + xi_line_this = lines_dictionary[line]['L_Hbeta_ratio'] + xi_line = (xi_line_this * XI_LINE) / SUM_RATIO_LINES_BLR + R_line_this = R_line * lines_dictionary[line]['R_Hbeta_ratio'] + blr_dict[line] = SphericalShellBLR(L_disk, xi_line, line, R_line_this) + ec_blr_dict[line] = Absorption(blr_dict[line], r=r_blob_bh, z=z) + tau_z_this = ec_blr_dict[line].tau_blr_cubepy(nu, eps_abs=1e-3) + tau_z_all += tau_z_this + assert True + + def sigma_pp(b): """pair production cross section""" diff --git a/environment.yml b/environment.yml index f2bfd43..43afe19 100644 --- a/environment.yml +++ b/environment.yml @@ -17,3 +17,4 @@ dependencies: - pip: - agnpy - gammapy + - cubepy<=1.0.2 diff --git a/setup.py b/setup.py index 683badf..82dc14b 100644 --- a/setup.py +++ b/setup.py @@ -20,6 +20,6 @@ "License :: OSI Approved :: BSD License", "Operating System :: OS Independent", ], - install_requires=["astropy>=5.0, <6.0", "numpy>=1.21", "scipy>=1.5,!=1.10", "matplotlib>=3.4"], + install_requires=["astropy>=5.0, <6.0", "numpy>=1.21", "scipy>=1.5,!=1.10", "matplotlib>=3.4", "cubepy<=1.0.2"], python_requires=">=3.8", )