Calculation of BLR absorption fails for r = R_line

[cosimoNigro]

The calculation of the absorption (and I presume also of the External Compton) on the SphericalShellBLR fails when r=R_line.

This is due to the definition of x (the distance from the reprocessing material to the blob), from Finke 2016:
x^2 = r^2 + R_line^2 - 2 * r * R_re * mu_re
where

• r = height of the blob above the BH;
• R_re = radial coordinate of the reprocessing material;
• mu_re = zenith angle of the reprocessing material.

Clearly when r=R_re, x=0 and the absorption (or EC) calculations gives nan, as the optical depth (or the SED) is dependent on x^(-2).

Here a small snippet to trigger the issue:

import numpy as np
import astropy.units as u
import astropy.constants as const
from agnpy.targets import SphericalShellBLR
from agnpy.absorption import Absorption

# BLR definition
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)

# absorption
abs_blr = Absorption(blr, r=R_line)

# tau
nu = np.logspace(22, 28) * u.Hz
tau_blr = abs_blr.tau(nu)

print(tau_blr)

returning:

[nan nan nan nan nan nan nan nan nan nan nan nan nan nan nan nan nan nan
 nan nan nan nan nan nan nan nan nan nan nan nan nan nan nan nan nan nan
 nan nan nan nan nan nan nan nan nan nan nan nan nan nan]

@jsitarek, shall we put a check in the x calculation?

[jsitarek]

Hi, I think this is the same thing which I mentioned in the second paragraph of #86 and this can also get triggered in other cases.
I can make a fix, but there are a few different ways to do it:

1. set the value of integrant to 0 (or some bordering value) for the points where you end up with position ~R_line
2. shift the actual integration points that are ~R_line a bit e.g. 0.1 %)

I think the second solution is probably easier. Any preference?

[cosimoNigro]

I agree, the second solution is easier, and probably also computationally faster.

[jsitarek]

added a PR #91 that should solve this.
I also made a fix for the case of mu_s !=1 - in that case we would not end up with nans usually because you would need to hit not only position of R, but also of the angles, but either way if one of integration points lies on the BLR sphere it can make the integration messy so I shifted those points as well.

There might be also an issue with the IC part, I guess it would also have the same problem, but a bit simpler, because you do not integrate over the photon path, so you could have a problem only in r = R_line