RFC: neff for large step-index fibre

[chrisbrahms]

This is an attempt to fix #372

The problem with very large step-index fibres is that the w parameter becomes very large:

Luna.jl/src/StepIndexFibre.jl

Line 141 in 38ce40b

w = rk0*sqrt(neff^2 - nclad^2)

Then in two terms in the characteristic equations you need to calculate the ratio of modified Bessel functions:

Luna.jl/src/StepIndexFibre.jl

Line 147 in 38ce40b

kpart = ((ncore^2 - nclad^2)/(2*ncore^2))^2*(besselkp(n, w)/(w*besselk(n, w)))^2

For large w (>800), one or both of besselk and besselkp evaluates to zero within floating-point precision, so we cannot compute the ratio properly.

This PR fixes this by adopting an approximate expression for the ratio of besselkp and besselk when either are zero. The basic point is that K_n(x) approaches the same limit regardless of the value of n: https://dlmf.nist.gov/10.25#E3

So then the ratio of K_n and K_(n-1) approaches 1, which in turn leads to the approximate expression here.

With this I can run the following example, which is basically #372

using Luna
import PyPlot: plt
import Roots: find_zeros

a = 300e-6 # With the other parameters, the highest I can set this to is 5e-6
NA = 0.22
flength = 1
fr = 0.18
τfwhm = 200e-15
λ0 = 250e-9
energy = 1e-6

grid = Grid.EnvGrid(flength, λ0, (150e-9, 450e-9), 1e-12)

m = StepIndexFibre.StepIndexMode(a, NA)#, accellims=(150e-9, 450e-9, 100))

ω0 = PhysData.wlfreq(λ0)
nclad = m.cladn(ω0; z=0)
ncore = m.coren(ω0; z=0)
char = StepIndexFibre.make_char(a, 2π/λ0, ncore, nclad, m.n, m.kind)
z = find_zeros(char, nclad, ncore; no_pts=100, xatol=0)
##
plt.figure()
nrange = collect(range(nclad, ncore, 100000))
plt.plot(nrange, char.(nrange))
plt.scatter(z, zero(z); color="r", s=5)
plt.axhline(0; linestyle="--", color="k")
plt.ylim(-1, 1)
plt.xlim(nclad, ncore)
##
plt.figure()
nrange = collect(range(ncore-0.000001, ncore, 100000))
plt.plot(nrange, char.(nrange))
plt.axhline(0; linestyle="--", color="k", alpha=0.2)
plt.scatter(z, zero(z); color="r", s=5)
plt.ylim(-1, 1)
plt.xlim(extrema(nrange)...)

And I get these plots of the characteristic equation and its roots. Note that I had to set xatol=0 in find_zeros because the roots get very close together.