Absorption with cubepy

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,112 @@ 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 +587,75 @@ 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_lines_cubepy(self, nu, lines_list=lines_dictionary.keys(), 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)
+
+        # 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 = SphericalShellBLR(self.target.L_disk, xi_line, line_key, R_line_current)
+            absorption_shells = Absorption(blr_shells, r=self.r, z=self.z)
+
+            tau_z_current = absorption_shells.tau_blr_cubepy(nu, eps_abs=eps_abs)
+            tau_z_lines += tau_z_current
+
+        return tau_z_lines
+
     @staticmethod
     def evaluate_tau_dt(
         nu,

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,49 @@ 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 test_tau_accuracy_standard_vs_cubepy(self):
+        """
+        This test aims to compare the calculation of tau using a standard method against
+        the implementation provided by the `cubepy` library.
+        The goal is to ensure that both methods yield consistent results within an acceptable tolerance.
+        """
+        L_disk = 2e45 * u.Unit("erg s-1")
+        xi_line = 0.024
+        R_line = 1e17 * u.cm
+        blr = SphericalShellBLR(L_disk, xi_line, "Lyalpha", R_line)
+        r = 2e17 * u.cm
+        z = 0.859
+        abs_blr = Absorption(blr, r, z)
+        E = np.logspace(2, 6) * u.GeV
+        nu = E.to("Hz", equivalencies=u.spectral())
+        tau_blr = abs_blr.tau(nu)
+        tau_blr_cubepy = abs_blr.tau_blr_cubepy(nu, eps_abs=1e-3)
+        assert np.allclose(tau_blr, tau_blr_cubepy, rtol=1e-1, atol=1e-1)
+
 
 def sigma_pp(b):
     """pair production cross section"""

environment.yml

@@ -13,4 +13,7 @@ dependencies:
   - matplotlib>=3.4,<=3.9
   - sherpa
   - pre-commit
-  - gammapy<1.2
+  - pip:
+    - gammapy<1.2
+    - cubepy<=1.0.2
+

setup.py

@@ -20,6 +20,7 @@
         "License :: OSI Approved :: BSD License",
         "Operating System :: OS Independent",
     ],
-    install_requires=["astropy>=5.0,<6.0", "numpy>=1.21", "scipy>=1.5,<1.10", "pyyaml", "matplotlib>=3.4,<=3.9", "sherpa", "pre-commit", "gammapy<1.2"],
-    python_requires=">=3.9,<3.12",
+
+    install_requires=["astropy>=5.0,<6.0", "numpy>=1.21", "scipy>=1.5,<1.10", "pyyaml", "matplotlib>=3.4,<=3.9", "sherpa", "pre-commit", "gammapy<1.2",  "cubepy<=1.0.2"],
+    python_requires=">=3.9,<3.12"
 )