C++ get_attenuation occasionally fails and returns nans

[sjoerd-bouma]

It seems that for some geometries, the C++ analytic raytracer returns np.nan for the attenuation (at least when using SP1 as the attenuation model), which then results in traces that are all nans. Using the Python version works. I'm a bit surprised as I hadn't run into this before, but also can't seem to find an earlier version of NuRadioMC on which this does not happen.

MWE below:

from NuRadioMC.SignalProp import analyticraytracing 
from NuRadioMC.utilities import medium
import numpy as np

use_cpp = {'C++':True, 'Python':False}
x_channel = [17.3, -10.0, -142.0]
x_nu = [210.75915126907, 267.081824987715, -86.5823449330406]

for key, value in use_cpp.items():
    print(f"Using {key} raytracer")
    analyticraytracing.cpp_available = value
    rt = analyticraytracing.ray_tracing(medium.get_ice_model('southpole_2015'), attenuation_model='SP1',)
    rt.set_start_and_end_point(x_nu, x_channel)
    rt.find_solutions()

    for iS in range(rt.get_number_of_solutions()):
        attenuation = rt.get_attenuation(iS, np.fft.rfftfreq(2048, .2), 5)
        if np.any(np.isnan(attenuation)):
            print('attenuation contains nans!', attenuation)
        else:
            print('no nans! Hooray!')

Not sure if this is similar to the issue in #521?

[colemanalan]

I tracked this down to a failure in the GSL status when it tries to run gsl_integration_qags, here.

However, if you look at the actual value of result there, the values seem to be fine. So maybe the integrator worried that it has not converged, but they actually have converged.

Potentially, we can just return exp(-1 * result) regardless of the GSL status.

[colemanalan]

Actually, we may have some numerical instabilities here. I was looking at how the result for the first 3 freq bins are affected by changing z. If you use

-85.5823449330406

• 0.0377103
• 0.0875237
• 0.10958

For -86.5823449330406 (the one that gives nans)

• 0.0320964
• 0.0744749
• 0.0932363

And for -87.5823449330406

• 0.0378141
• 0.0877644
• 0.109881

So clearly the bad values aren't bounded properly. Probably a fix of this will involve finding the weirdness of the integrand that cuases this to fail or to slightly perturb the start/end points (which hopefully fixes the issue) and return some kind of interpolation