Nan Output 7-10 kHz

[DavidGoldak]

Hi Dr. Gasdia,

Check out the attached screenshot, this is your "getting started" example, but using a different transmitter frequency. As you can see, a frequency of 8 kHz results in Nan output from "propagate" but 7 kHz works (as does frequencies above 10 kHz).

I've found that this almost always occurs whether I use wait profile, ferguson profile, FIRI model, distance to source doesn't seem to matter, number of segments in a segmented waveguide doesn't seem to matter, using a single homogenous waveguide doesn't help either. Sometimes night-time seems to be better than day?

But I very often get Nan's between ~ 8 - 10 kHz inclusively, so 7 kHz and below is usually okay, and 11 kHz and above is usually okay.

Any ideas what might be going on, or a possible solution to fix?

Thanks,

David

[JamesMCannon]

I can confirm I also have this issue using just the Running your first model, although for me it spans [7.074, 11.999] kHz. I found downgrading from LongwaveModePropagator v0.4.1 to v0.4.0 eliminates this issue. This is likely related to the changes in 0.4.1 to use a RootsAndPoles.rectanglemesh at frequencies less than 12 kHz.

[fgasdia]

Thanks for the extra info @JamesMCannon! I don't know how I didn't run into this when I selected the 12 kHz...

I'm working on a major update to RootsAndPoles which will require me to change/simplify the defaultmesh function anyways, but for now a simple workaround is for the user to provide their own mesh.

In the example code below, I build a rectangular mesh that goes down to 1 deg in the real axis and -60 deg in the imaginary axis.

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

using Plots
using Printf

# Generate mesh using a higher initial resolution than the default for < 12 kHz
# (parameters are what's used for above 12 kHz)
rmin, rmax = deg2rad(1), deg2rad(89.9)
imin, imax = deg2rad(-60), deg2rad(0.0)
Δr_coarse = deg2rad(0.5)
Δr_fine = deg2rad(0.1)
rtransition, itransition = deg2rad(75), deg2rad(-1.5)
zbl_coarse = complex(rmin, imin)
ztr_coarse = complex(rmax, imax)

mesh = LMP.rectangulardomain(zbl_coarse, ztr_coarse, Δr_coarse)
filter!(z->(real(z) < rtransition || imag(z) < itransition), mesh)
zbl_fine = complex(rtransition, itransition)
ztr_fine = complex(rmax, imax)
append!(mesh, LMP.rectangulardomain(zbl_fine, ztr_fine, Δr_fine))

# Set up the waveguide
rx = GroundSampler(0:5e3:2000e3, Fields.Ez)
bfield = BField(50e-6, π/2, 0)
electrons = Species(QE, ME, z->waitprofile(z, 75, 0.35), electroncollisionfrequency)
ground = Ground(10, 1e-4)
waveguide = HomogeneousWaveguide(bfield, electrons, ground)

# Try several transmitter frequencies
freqs = [3e3, 6e3, 9e3, 12e3, 15e3, 18e3]
as = Matrix{Float64}(undef, length(rx.distance), length(freqs))
ps = similar(as)
for (i, f) in pairs(freqs)
    tx = Transmitter(f)
    E, a, p = propagate(waveguide, tx, rx; mesh)
    as[:,i] .= a
    ps[:,i] .= p
end

plot(rx.distance/1000, as, label=permutedims(string.(freqs/1000)), legend_title="freq (kHz)",
     xlabel="Distance (km)", ylabel="Amplitude (dB μV/m)")

The resulting plot is:

What we see is explained a little bit in the docs pages: Meshgrid Part 1 and Meshgrid Part 2.

The default mesh uses a higher resolution at large real and small (negative) imaginary angles because at high frequencies the roots and poles of the mode equation are closely spaced there. We could use a high resolution throughout the full domain, but that is inefficient (many integrations are being done unnecessarily). The default at frequencies above 12 kHz also cuts off the large real/large (negative) imaginary part of the domain because there are no eigenangles there. At lower frequencies, the modes are not closely separated so there is no need for a higher resolution at large real/small imaginary angles. However, there are modes at larger (negative) imaginary angles in general, so the mesh needs to cover a larger domain.

I suspect the NaNs are coming from overflow of the Modified Hankel function. You may continue to see this warning when you run the example above:

┌ Warning: imaginary component less than -0.5410 rad (-31°) may cause wavefields calculated with modified Hankel functions to overflow.
└ @ LongwaveModePropagator ~/dev/LongwaveModePropagator.jl/src/LongwaveModePropagator.jl:201

but the results are still good. If you make imin even more negative, you'll eventually get the NaN's again.

I also ran the above example at all frequencies using a copy-paste of the mesh default for >= 12 kHz.

rmin = deg2rad(30)
imin = deg2rad(-10)

mesh = trianglemesh(zbl_coarse, ztr_coarse, Δr_coarse)
filter!(z->(real(z) < rtransition || imag(z) < itransition), mesh)
zbl_fine = complex(rtransition, itransition)
ztr_fine = complex(rmax, imax)
append!(mesh, trianglemesh(zbl_fine, ztr_fine, Δr_fine))

All the frequencies except for 3 kHz match the results above. At 3 kHz, there are only a few widely-spaced modes. With this smaller domain we aren't seeing many modes and get an amplitude curve which only captures the lowest order features. Overall, this suggests the threshold in defaultmesh should be less than 12 kHz (probably 6 to 7 kHz).

[JamesMCannon]

With an update coming to meshing in general, this likely a moot point but perhaps the issue relates to this block in modefinder.jl

elseif meshshape == "rectanglemesh"
            mesh = rectangulardomain(zbl_coarse, ztr_coarse, Δr_coarse)

            filter!(z->(real(z) < rtransition || imag(z) < itransition), mesh)
            zbl_fine = complex(rtransition, itransition)
            ztr_fine = complex(rmax, imax)
            append!(mesh, rectangulardomain(zbl_coarse, ztr_coarse, Δr_fine))
end

Is it a typo that with a rectangular mesh, you append the zbl_coarse, and ztr_coarse rather than zbl_fine and ztr_fine?