Closed hagenw closed 7 years ago
The number of secondary sources (times 0.5) is used as maximum modal order. The imaginary part of the Hankel-function blows up. For values larger than double precision range, sph_jnyn seems to return 0.
Maybe it is then a good idea to check if this always happens for the same maximum modal order and add a warning or better error message to the code.
From what I've gathered:
sph_jnyn is deprecated in scipy 0.18.0 The new spherical_yn properly returns -inf for the blown-up Neumann function.
Unfortunately, this does not solve the problem:
In Matlab, large orders are fail-safe because "division by blown-up Hankel yields 0" works. This does not seem to be the case with Python's complex type, e.g.:
Matlab: a*(x + 1j*inf) = (a*x + 1j*inf)
Python: a*(x + 1j*inf) = (nan + 1j*inf)
Bottom line: A reasonable maximum order should be used.
Start with these settings
First, try this one, this should still work:
Now increase the number of secondary sources and we will get a division by zero: