kappa = 2 or 3 fails

[sigvaldm]

Electrons with kappa = 2 or 3 leads to divide-by-zero when using finite_radius_current

[sigvaldm]

On closer inspection, kappa<4 is not valid for finite_radius_current, because it is outside the domain of tabulated values. kappe=4, kappa=6, and kappa=inf are the tabulation points for kappa. However, the error also occur in OML_current, which, should be valid for kappa>D/2+1, where D is the number of dimensions, e.g. kappa>2.5 for spheres, and kappa>2 for cylinders. finite_radius_current uses OML_current internally, so probably it is enough to fix OML_current.

[sigvaldm]

Here is an example:

In [82]: OML_current(Sphere(r=1e-2), Electron(kappa=3+1e-5), eta=10)
Out[82]: -1.305821017623522e-05

In [83]: OML_current(Sphere(r=1e-2), Electron(kappa=3-1e-5), eta=10)
Out[83]: -1.3058236835311985e-05

In [84]: OML_current(Sphere(r=1e-2), Electron(kappa=3), eta=10)
---------------------------------------------------------------------------
ZeroDivisionError                         Traceback (most recent call last)
<ipython-input-84-9e851f1ac3f1> in <module>
----> 1 OML_current(Sphere(r=1e-2), Electron(kappa=3), eta=10)

~/Codes/langmuir/langmuir/analytical.py in OML_current(geometry, species, V, eta, normalization)
    195     else:
    196         C = np.sqrt(kappa-1.5)*gamma(kappa-1.)/gamma(kappa-0.5)
--> 197         D = (1.+24*alpha*((kappa-1.5)**2/((kappa-2.)*(kappa-3.))))/(1.+15*alpha*((kappa-1.5)/(kappa-2.5)))
    198         E = 4.*alpha*kappa*(kappa-1.)/( (kappa-2.)*(kappa-3.)+24*alpha*(kappa-1.5)**2 )
    199         F = ((kappa-1.)*(kappa-2.)*(kappa-3.)+8*alpha*(kappa-3.)*(kappa-1.5)**2) /( (kappa-2.)*(kappa-3.)*(kappa-1.5)+24*alpha*(kappa-1.5)**3 )

ZeroDivisionError: float division by zero

Although there is a divide-by-zero for kappa=3, the analytical expression used in OML_current clearly has a finite value when taking the limit as kappa tends towards 3. Looking at the analytical expressions in the paper, there are a lot of denominators with expressions like (kappa-3), (kappa-2.5) and so forth.

Perhaps one way to proceed would be to first make pytests that OML_current and finite_radius_current shouldn't throw divide-by-zero when exposed to such kappa-values, and then to correct OML_current by implementing limiting cases. Tests should exist for both sphere and cylinder, and both positive and negative eta.