Numerical (round-off) errors/accuracy for a SegmentedWaveguide with identical neighboring Homogeneous segments

[fgasdia]

There sometimes appear to be errors in the amplitude and phase that begin at the boundary of segments in a SegmentedWaveguide when both sides of the transition are identical HomogeneousWaveguides. I suspect this is due to numerical errors produced by the mode conversion process when the eigenangles are approximately identical on both sides.

We could precheck for identical segments and merge them into a single longer HomogeneousWaveguide (or simply treat them this way in the mode sum), but this may not be possible in general. In the meantime, users should avoid using neighboring identical segments in SegmentedWaveguides.

[fgasdia]

I haven't been able to replicate this using sequential identical segments, but I have observed similar "bumps" in amplitude for some changing SegmentedWaveguides. It usually doesn't happen, but e5e88e9 has some scenarios where it occurs.

Work needed:

• identify what scenarios cause this to happen (it may be a secondary effect having to do with reflection coefficients)
• identify where the problem is actually occurring in the code. Is the problem localized to single distances or a lasting effect for the segment?
• Figure out a fix (could be a coding bug, physics problem, or implementation/numerical problem)

[fgasdia]

This may be a modelism that is a byproduct of the sharp boundaries between waveguide segments. I compared the electric field components (ignore the likely magnitude error) between FDTD and LongwaveModePropagator using the code below. There are 4 waveguide segments where the ionosphere and ground conductivity changes in each. There are small steps in the amplitude profile of Ex at the boundaries at 1000km and 1500km in both LMP and FDTD, and a step in Ey for LMP only.

using Dates
using LongwaveModePropagator
using LongwaveModePropagator: C0, QE, ME
const LMP = LongwaveModePropagator
using PropagationModelPrep
using LsqFit
using Plots, Printf, Statistics
using Plots.PlotMeasures

const TX = Transmitter(24e3)
const GS = GroundSampler(0:500:2000e3, Fields.E)
const SEGMENT_RANGES = [0, 500e3, 1000e3, 1500e3]
const GROUND_SIGMAS1 = [0.001, 0.001, 0.0005, 0.0001]
const GROUND_EPSRS1 = [4, 4, 10, 10]
const GROUND_SIGMAS2 = [0.001, 0.001, 0.001, 0.0001]
const GROUND_EPSRS2 = [4, 4, 4, 10]
const HPRIMES = [72.0, 74, 75, 76]
const BETAS = [0.28, 0.30, 0.32, 0.35]
const BFIELD = BField(50e-6, π/2, 0)

"Generate input file for FDTD"
function generate_fullfields()
    nsegments = length(SEGMENT_RANGES)

    input = ExponentialInput()
    input.name = "fullfields"
    input.description = "4 segments"
    input.datetime = Dates.now()

    input.segment_ranges = SEGMENT_RANGES
    input.hprimes = HPRIMES
    input.betas = BETAS
    input.b_mags = fill(BFIELD.B, nsegments)
    input.b_dips = fill(LMP.dip(BFIELD), nsegments)
    input.b_azs = fill(LMP.azimuth(BFIELD), nsegments)
    input.ground_sigmas = GROUND_SIGMAS1
    input.ground_epsrs = GROUND_EPSRS1
    input.frequency = TX.frequency.f
    input.output_ranges = collect(0:2e3:2000e3)

    return input
end

"Process FDTD output files"
function processfields(inputs, dfts::EMP2D.DFTFields)
    freq = first(dfts.DFTfreqs)
    freq_kHz = freq/1e3

    # Strip last index of each field because they're 0 (and then inf)
    dist = round.(dfts.dist[1:end-1], digits=3)  # fix floating point
    dist_km = dist./1e3

    drange_km = inputs.drange/1e3

    # 2D FDTD numerical dispersion correction - is this correct for Ep and Hr, Ht?
    p = [0.00136588, -0.026108, 0.154035]
    k = p[1]*freq_kHz^2 + p[2]*freq_kHz + p[3]  # deg/km²

    out = EMP2D.DFTFields(Vector{Float64})
    out.dist = copy(dist)
    out.DFTfreqs = copy(dfts.DFTfreqs)

    # Scale amplitude and phase of each field component
    for f in (:Er, :Et, :Ep, :Hr, :Ht, :Hp)
        amp = 20*log10.(getfield(dfts, f).amp[1:end-1]/1e-6)  # dB μV/m

        phase = rad2deg.(getfield(dfts, f).phase[1:end-1]) + (360*freq/C0*dist)
        phase += k*dist_km*drange_km^2

        setfield!(out, f, EMP2D.Phasor(amp, phase))
    end

    return out
end

"Build and run LMP"
function buildrunwaveguide(ground_epsrs, ground_sigmas)
    wvgs = Vector{HomogeneousWaveguide{Species}}(undef, length(SEGMENT_RANGES))
    for i in eachindex(SEGMENT_RANGES)
        ground = Ground(ground_epsrs[i], ground_sigmas[i])
        species = Species(QE, ME, z->waitprofile(z, HPRIMES[i], BETAS[i]), electroncollisionfrequency)

        wvgs[i] = HomogeneousWaveguide(BFIELD, species, ground, SEGMENT_RANGES[i])
    end
    wvg = SegmentedWaveguide(wvgs)

    propagate(wvg, TX, GS)
end

path = "fullfields"
s = generate_fullfields()
inputs = EMP2D.Inputs(2500e3, (s.frequency,))
inputs.savefields = Int32[0, 0, 0, 0, 0, 0]

# exepath = "fakepath"
# computejob = LocalParallel("fullfields", path, exepath, 4, 120)
# EMP2D.run(s, computejob; inputs)

dfts = EMP2D.readdfts(path)
out = processfields(inputs, dfts)

E, a, p = buildrunwaveguide(GROUND_EPSRS1, GROUND_SIGMAS1)
E2, a2, p2 = buildrunwaveguide(GROUND_EPSRS2, GROUND_SIGMAS2)

# Adjust bias of EMP2D to minimize L2 distance to LMP
lastidx = findfirst(x->x==2000e3, out.dist)
compmask = 2:lastidx  # exclude nan of LMP

model(x, p) = x .+ only(p)
fit = curve_fit(model, out.Er.amp[compmask], a[2:end,1], [0.0])
dEr = only(fit.param)
fit = curve_fit(model, out.Ep.amp[compmask], a[2:end,2], [0.0])
dEp = only(fit.param)
fit = curve_fit(model, out.Et.amp[compmask], a[2:end,3], [0.0])
dEt = only(fit.param)

plot(out.dist[compmask]/1000, out.Er.amp[compmask].+dEr, ylims=(-70, 150), label="Er"*@sprintf("%+.1f",dEr))
plot!(out.dist[compmask]/1000, out.Ep.amp[compmask].+dEp, label="Ep"*@sprintf("%+.1f",dEp))
plot!(out.dist[compmask]/1000, out.Et.amp[compmask].+dEt, label="Et"*@sprintf("%+.1f",dEt))
plot!(out.dist[compmask]/1000, a[2:end,1], label="LMP_Ez")
plot!(out.dist[compmask]/1000, a[2:end,2], label="LMP_Ey")
plot!(out.dist[compmask]/1000, a[2:end,3], label="LMP_Ex")
plot!(xlabel="range (km)", ylabel="amplitude (dB uV/m)", top_margin=6mm)
annotate!(0, 155, (string(GROUND_EPSRS1[1])*", "*string(GROUND_SIGMAS1[1]), 8, :left))
annotate!(500, 155, (string(GROUND_EPSRS1[2])*", "*string(GROUND_SIGMAS1[2]), 8, :left))
annotate!(1000, 155, (string(GROUND_EPSRS1[3])*", "*string(GROUND_SIGMAS1[3]), 8, :left))
annotate!(1500, 155, (string(GROUND_EPSRS1[4])*", "*string(GROUND_SIGMAS1[4]), 8, :left))
annotate!(0, 165, ("FDTD & LMP (ϵᵣ, σ)", 8, :left))
savefig("identicalground.png")

Then I ran LMP again but set the third waveguide segment (beginning at 1000km) with a ground identical to the ground in the second component.

plot(out.dist[compmask]/1000, out.Er.amp[compmask].+dEr, ylims=(-70, 150), label="Er"*@sprintf("%+.1f",dEr))
plot!(out.dist[compmask]/1000, out.Ep.amp[compmask].+dEp, label="Ep"*@sprintf("%+.1f",dEp))
plot!(out.dist[compmask]/1000, out.Et.amp[compmask].+dEt, label="Et"*@sprintf("%+.1f",dEt))
plot!(out.dist[compmask]/1000, a2[2:end,1], label="LMP_Ez")
plot!(out.dist[compmask]/1000, a2[2:end,2], label="LMP_Ey")
plot!(out.dist[compmask]/1000, a2[2:end,3], label="LMP_Ex")
plot!(xlabel="range (km)", ylabel="amplitude (dB uV/m)", top_margin=6mm)
annotate!(0, 155, (string(GROUND_EPSRS2[1])*", "*string(GROUND_SIGMAS2[1]), 8, :left))
annotate!(500, 155, (string(GROUND_EPSRS2[2])*", "*string(GROUND_SIGMAS2[2]), 8, :left))
annotate!(1000, 155, (string(GROUND_EPSRS2[3])*", "*string(GROUND_SIGMAS2[3]), 8, :left))
annotate!(1500, 155, (string(GROUND_EPSRS2[4])*", "*string(GROUND_SIGMAS2[4]), 8, :left))
annotate!(0, 165, ("LMP (ϵᵣ, σ)", 8, :left))
savefig("differentground.png")

The amplitude jumps at 1000km are no longer present.