From 435d853365bfbae568c197922f562551cb3c02c0 Mon Sep 17 00:00:00 2001 From: Julian Sitarek Date: Wed, 28 Apr 2021 13:10:00 +0200 Subject: [PATCH] added a special treatment for absorption on BLR if one of the grid points ends up very close to the shell also a test to check that it is working --- agnpy/absorption/absorption.py | 24 ++++++++++++++- agnpy/tests/test_absorption.py | 56 ++++++++++++++++++++++++++++++++++ agnpy/utils/math.py | 3 ++ 3 files changed, 82 insertions(+), 1 deletion(-) diff --git a/agnpy/absorption/absorption.py b/agnpy/absorption/absorption.py index 51ab54b2..f4175e0d 100644 --- a/agnpy/absorption/absorption.py +++ b/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 diff --git a/agnpy/tests/test_absorption.py b/agnpy/tests/test_absorption.py index 2035fe0c..cc4d416a 100644 --- a/agnpy/tests/test_absorption.py +++ b/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""" diff --git a/agnpy/utils/math.py b/agnpy/utils/math.py index ed49c586..ea489141 100644 --- a/agnpy/utils/math.py +++ b/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