In using SciPy's scipy.special.ellipeinc and ellipkinc, there seem to be some islands of numerical instability. For example,
>>> from scipy.special import ellipkinc
>>> ellipkinc(0.9272952180016123, 0.68359375000000011)
nan
>>> ellipkinc(0.9272952180016123, 0.6835937500000002)
2.0518660200390668
>>> ellipkinc(0.9272952180016123, 0.68359375)
1.0259330100195332
>>> ellipkinc(0.9272952180016123, 0.68359374)
1.0259330081166262
This happens where k^2.sin^2(phi) is close to 0.3, but there is nothing unusual about the elliptic integrals themselves here, so presumably it's a numerical thing. I don't know enough about the inner workings of this algorithm to say what's wrong, so what's my best option? I thought of:
round(0.68359375000000011,8)
say, but this is surely going to slow my code?
(This is more an extended comment about the problem than an answer to what to do about it.)
This looks like a bug in ellipkinc. I get just one floating point value where the function returns nan, and four adjacent floating point values where the function returns twice the "correct" value:
In [91]: phi = 0.9272952180016123
In [92]: mbad = 0.68359375000000011
In [93]: m = np.nextafter(mbad, 0)
In [94]: mvals = []
In [95]: for j in range(10):
....: mvals.append(m)
....: m = np.nextafter(m, 1)
....:
In [96]: mvals
Out[96]:
[0.68359375,
0.68359375000000011,
0.68359375000000022,
0.68359375000000033,
0.68359375000000044,
0.68359375000000056,
0.68359375000000067,
0.68359375000000078,
0.68359375000000089,
0.683593750000001]
In [97]: ellipkinc(phi, mvals)
Out[97]:
array([ 1.02593301, nan, 2.05186602, 2.05186602, 2.05186602,
2.05186602, 1.02593301, 1.02593301, 1.02593301, 1.02593301])