Merge pull request #91 from jsitarek/abs_blr_RLine_bug

added a special treatment for absorption on BLR if one of the grid points ends up very close to the shell

agnpy/absorption/absorption.py

@@ -5,7 +5,13 @@
 from astropy.io import fits
 from astropy.constants import c, G, h, m_e, M_sun, sigma_T
 from scipy.interpolate import interp2d
-from ..utils.math import axes_reshaper, log, mu_to_integrate, phi_to_integrate
+from ..utils.math import (
+    axes_reshaper,
+    log,
+    mu_to_integrate,
+    phi_to_integrate,
+    min_rel_distance,
+)
 from ..utils.geometry import (
     cos_psi,
     x_re_shell,
@@ -321,6 +327,12 @@ def evaluate_tau_blr(
         epsilon_1 = nu_to_epsilon_prime(nu, z)
         # multidimensional integration
         l = np.logspace(0, 5, l_size) * r
+
+        # check if any point is too close to R_line, the function works only for mu=1, so
+        # we can check directly if R_line is within 'l' array
+        idx = np.isclose(l, R_line, rtol=min_rel_distance)
+        l[idx] += min_rel_distance * R_line
+
         _mu, _phi, _l, _epsilon_1 = axes_reshaper(mu, phi, l, epsilon_1)
         x = x_re_shell(_mu, R_line, _l)
         _mu_star = mu_star_shell(_mu, R_line, _l)
@@ -387,6 +399,16 @@ def evaluate_tau_blr_mu_s(
         # multidimensional integration
         # here uu is the distance that the photon traversed
         uu = np.logspace(-5, 5, u_size) * r
+
+        # check if for any uu value the position of the photon is too close to the BLR
+        x_cross = np.sqrt(r ** 2 + uu ** 2 + 2 * uu * r * mu_s)
+        idx = np.isclose(x_cross, R_line, rtol=min_rel_distance)
+        if idx.any():
+            uu[idx] += min_rel_distance * R_line
+            # it might happen that some of the points get more shifted then the next one,
+            # possibly making integration messy, so we sort the points
+            uu = np.sort(uu)
+
         _mu_re, _phi_re, _u, _epsilon_1 = axes_reshaper(mu, phi, uu, epsilon_1)
 
         # distance between soft photon and gamma ray

agnpy/tests/test_absorption.py

@@ -374,6 +374,62 @@ def test_abs_blr_mu_s_vs_on_axis(self, r_to_R):
         # which are probably due to numerical uncertainties in the integrals
         assert check_deviation(nu, tau_blr, tau_blr_on_axis, 0.25, x_range=xrange)
 
+    @pytest.mark.parametrize("r_R_line", [1, 1.0e5 ** (-10.0 / 49)])
+    def test_tau_blr_Rline(self, r_R_line):
+        """
+        Checks if absorption on BLR works also fine
+        if one of the integration points falls on R_line
+        """
+        L_disk = 2 * 1e46 * u.Unit("erg s-1")
+        xi_line = 0.024
+        R_line = 1.1 * 1e17 * u.cm
+        blr = SphericalShellBLR(L_disk, xi_line, "Lyalpha", R_line)
+
+        abs_blr = Absorption(blr, r=r_R_line * R_line)
+
+        # reference (without problem)
+        abs_blr0 = Absorption(blr, r=r_R_line * R_line * 1.001)
+
+        nu = np.logspace(22, 28) * u.Hz
+        tau_blr = abs_blr.tau(nu)
+        tau_blr0 = abs_blr0.tau(nu)
+
+        # the current integrals are not very precise, and this includes
+        # tricky part to integrate, so allowing for rather large error margin
+        assert np.allclose(tau_blr, tau_blr0, atol=0.01, rtol=0.04)
+
+    def find_r_for_x_cross(mu_s, ipoint, npoints):
+        """
+        Finding which r should be used to fall with one of the points on top of BLR sphere.
+        """
+        # u = alpha * r
+        alpha = 1.0e-5 * 1.0e10 ** (ipoint / (npoints - 1.0))
+        return 1 / np.sqrt(alpha ** 2 + 2 * alpha * mu_s + 1)
+
+    @pytest.mark.parametrize("r_R_line", [1, find_r_for_x_cross(0.99, 60, 100)])
+    def test_tau_blr_mus_Rline(self, r_R_line):
+        """
+        Checks if absorption on BLR works also fine
+        if one of the integration points falls on R_line
+        """
+        L_disk = 2 * 1e46 * u.Unit("erg s-1")
+        xi_line = 0.024
+        R_line = 1.1 * 1e17 * u.cm
+        blr = SphericalShellBLR(L_disk, xi_line, "Lyalpha", R_line)
+
+        abs_blr = Absorption(blr, r=r_R_line * R_line, mu_s=0.99)
+
+        # reference (without problem)
+        abs_blr0 = Absorption(blr, r=r_R_line * R_line * 1.001, mu_s=0.99)
+
+        nu = np.logspace(22, 28) * u.Hz
+        tau_blr = abs_blr.tau(nu)
+        tau_blr0 = abs_blr0.tau(nu)
+
+        # the current integrals are not very precise, and this includes
+        # tricky part to integrate, so allowing for rather large error margin
+        assert np.allclose(tau_blr, tau_blr0, atol=0.01, rtol=1.04)
+
 
 class TestEBL:
     """class grouping all tests related to the EBL class"""

agnpy/utils/math.py

@@ -9,6 +9,9 @@
 mu_to_integrate = np.linspace(-1, 1, 100)
 phi_to_integrate = np.linspace(0, 2 * np.pi, 50)
 
+# minimum relative distance to the absorber (to avoid infinite integrals)
+min_rel_distance = 1.e-4
+
 # type of float to be used for math operation
 numpy_type = np.float64
 # smallest positive float