Modes not being identifed at low VLF and ELF

[fgasdia]

For some reason grpf from RootsAndPoles.jl, the mode finder used by LongwaveModePropagator.jl, is not identifying roots and poles at low VLF/ELF frequencies. See the example with a 1 kHz transmitter below.

using Plots

using LongwaveModePropagator
using LongwaveModePropagator: QE, ME
const LMP = LongwaveModePropagator

using RootsAndPoles

tx = Transmitter(1e3)
rx = GroundSampler(0:1e3:2000e3, Fields.Ez)

bfield = BField(5.5e-5, π/2, 0)
electrons = Species(QE, ME, z->waitprofile(z, 83, 0.5), electroncollisionfrequency)
ground = GROUND[3]

waveguide = HomogeneousWaveguide(bfield, electrons, ground)
modeequation = PhysicalModeEquation(tx.frequency, waveguide)

# GRPF `mesh` covers almost the full lower right quadrant of the complex plane
mesh = LMP.defaultmesh(tx.frequency;
    rmin=deg2rad(1), imin=deg2rad(-89), imax=deg2rad(0), rmax=deg2rad(89.9))

roots, poles = LMP.grpf(θ->LMP.solvemodalequation(θ, modeequation;
    susceptibilityfcn=z->LMP.susceptibility(z, modeequation)), mesh, LMPParams().grpfparams)

# Dense mesh
Δr = 0.2
x = 0:Δr:90
y = -90:Δr:0
densemesh = x .+ im*y';

function modeequationphase(me, mesh)
    phase = Vector{Float64}(undef, length(mesh))
    Threads.@threads for i in eachindex(mesh)
        f = LMP.solvemodalequation(deg2rad(mesh[i]), me)
        phase[i] = rad2deg(angle(f))
    end
    return phase
end

phase = modeequationphase(modeequation, densemesh);

img = heatmap(x, y, reshape(phase, length(x), length(y))';
        color=:twilight, clims=(-180, 180),
        xlims=(0, 90), ylims=(-90, 0),
        xlabel="real(θ)", ylabel="imag(θ)")

plot!(img, real(rad2deg.(roots)), imag(rad2deg.(roots)); color="red",
      seriestype=:scatter, markersize=5, label="roots");
plot!(img, real(rad2deg.(poles)), imag(rad2deg.(poles)); color="red",
      seriestype=:scatter, markershape=:utriangle, markersize=5, labels="poles")

[fgasdia]

This was due to defaultmesh using a triangle shaped mesh that doesn't include the lower right half of the quadrant. This is for efficiency at higher frequencies, but the whole quadrant needs to be searched at lower frequencies.

I additionally checked to see if topheight of LMPParams() needed to be raised, but even at nighttime, raising topheight makes a very small difference for a 100 Hz ELF transmitter.

using Plots

using LongwaveModePropagator
using LongwaveModePropagator: QE, ME
const LMP = LongwaveModePropagator

using FaradayInternationalReferenceIonosphere
const FIRI = FaradayInternationalReferenceIonosphere

using Interpolations

# Constants
ranges = 0:1e3:2000e3
rx = GroundSampler(ranges, Fields.Ez)
bfield = BField(5.5e-5, π/2, 0)

profile = firi(120, 45)
newprofile = FIRI.extrapolate(profile, 0:1e3:200e3)
itp =  Interpolations.interpolate(0:1e3:200e3, newprofile, FritschButlandMonotonicInterpolation())

electrons = Species(QE, ME, itp, electroncollisionfrequency)
ground = GROUND[5]

waveguide = HomogeneousWaveguide(bfield, electrons, ground)


# Varying topheight

tx = Transmitter(100)

function varytopheight(topheights, waveguide, tx, rx)
    amps = Matrix{Float64}(undef, length(ranges), length(topheights))
    phases = similar(amps)
    for i in eachindex(topheights)
        E, a, p = propagate(waveguide, tx, rx; params=LMPParams(topheight=topheights[i]))
        
        amps[:,i] = a
        phases[:,i] = p
    end
    return amps, phases
end

theights = (110e3, 130e3, 150e3, 180e3, 200e3)
amps, phases = varytopheight(theights, waveguide, tx, rx)

plot(rx.distance/1000, amps, labels=[theights...]'/1e3,
    xlabel="range (km)", ylabel="amplitude (dB)", title="top height (km)",
    linewidth=1.5)
plot(rx.distance/1000, rad2deg.(phases), labels=[theights...]',
    xlabel="range (km)", ylabel="phase (deg)",
    linewidth=1.5, legend=false)