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.

[JamesMCannon]

I can recreate this issue using a "real" path from NML to Gillam MB Canada.

     using LongwaveModePropagator
    using LongwaveModePropagator: QE, ME
    using Plots
    using Dates, DateFormats
    using LMPTools

    function WGSegmentation(tx,rx;x_res=20e3,dec_year=2024,gs_res=50e3,max_WG_length=400e3)
    #=Segments a waveguide by first checking for where the ground conductivity changes along a path from tx to rx.
    After ground conductivity, any waveguide segment longer than max_WG_length is split at max_WG_length from the previous segment.
    Background magnetic fields are calculated for each segement and a ground sampler is setup. 
    =#
    first_gnds, first_dists = groundsegments(tx,rx,resolution=x_res,require_permittivity_change=false)
    max_dist=range(tx,rx)

    append!(first_dists,max_dist)
    append!(first_gnds,first_gnds[[end]])

    max_idx = length(first_dists)
    dists=[]
    gnds=[]
    
    for i in 1:1:max_idx-1
        d, g = recursive_dists_segmentation(first_dists[i],first_dists[i+1],max_WG_length,first_gnds[i],first_gnds[i+1])
        if i > 1
            if dists[end]==d[1]
                    append!(dists,d[2:end])
                    append!(gnds,g[2:end])
            else
                    append!(dists,d)
                    append!(gnds,g)
            end
        else
            append!(dists,d)
            append!(gnds,g)
        end
    end

    pop!(dists)
    pop!(gnds)
    mag = igrf(tx,rx,dec_year,dists)
    num_segments = length(dists)

    sample_dists = collect(0:gs_res:max_dist)
    push!(sample_dists,max_dist)

    gs = GroundSampler(sample_dists,Fields.Ez)

    return num_segments, dists, gnds, mag, gs
end

function recursive_dists_segmentation(current,next,max_diff,curr_gnd,next_gnd)
    #Recusively segments the distances. Coded while tired, should be easy to re-write with a while loop
    diff_x = next-current
    if diff_x>max_diff
        dists,associated_gnds=recursive_dists_segmentation(current+max_diff,next,max_diff,curr_gnd,next_gnd)
        list = vcat(current,dists)
        gnds = vcat(curr_gnd,associated_gnds)
    else
        list = vcat(current,next)
        gnds = vcat(curr_gnd,next_gnd)
    end
    return list, gnds
end

    dec_year = yeardecimal(DateTime(2024,07,02))

    const LMP = LongwaveModePropagator
    
    rx = Receiver("GIL", 56.377043, -94.64403, 0.0, VerticalDipole())
    tx = Transmitter{VerticalDipole}("NML", 46.366014, -98.335544, VerticalDipole(), Frequency(25.2e3), 233000)

    x_res = 1e3
    gs_res = 1e3
    max_WG_dist=250e3 #20 wavelengths of 24.8 kHz

    println("Starting ",tx.name," to ",rx.name)
    num_segments, dists, gnds, mag, gs = WGSegmentation(tx,rx,dec_year=dec_year,x_res=x_res,gs_res=gs_res,max_WG_length=max_WG_dist)

   #gnds[6] = Ground(15,0.003)

    sim_params = LMPParams(grpfparams=GRPFParams(500, 500000, 3, 100000, 1e-5, true))

    # Set Wait profile parameters 
    h=75.7
    beta=0.37

    pha_plt=plot(dpi=400)
    amp_plt=plot(dpi=400)

    for h in 75:.1:76.5

        electrons = Species(QE, ME, z->waitprofile(z, h, beta,threshold=1e15), electroncollisionfrequency)

        waveguide = SegmentedWaveguide([HomogeneousWaveguide(mag[i], electrons, gnds[i],dists[i]) for i in 1:num_segments])

        # return the complex electric field, amplitude (dB uV/m), and phase (radians)
        Ez, az, pz = propagate(waveguide, tx, gs, params =sim_params);

        #convert amp to V/m
        az_uV = 10 .^(az./20)
        az_V = az_uV.*(10^-6)
        #convert pha to degrees
        pha_deg = pz .* (180/pi)
            
    #----------------------------------------------------------------------------
    
        global amp_plt=plot!(amp_plt,gs.distance/1000, az_V, line_z=h,label="",c=:plasma,xlabel="distance (km)", ylabel="amplitude (|E| V/m)",dpi=400,ylims = (1e-4, 1e1),xlims = (0,dists[end]/1000),yscale=:log10,minorgrid=true)
        global pha_plt=plot!(pha_plt,gs.distance/1000, pz,line_z=h,c=:plasma, xlabel="distance (km)", ylabel="phase (rad)",dpi=400,xlims = (0,dists[end]/1000),minorgrid=true,label="")
    end

    display(pha_plt)
    display(amp_plt)
   

Uncommenting the line gnds[6] = Ground(15,0.003) to force the second to last segment to match the last segment and running the same code results in a much smaller jump in amplitude/phase at the waveguide segment border.

[JamesMCannon]

Across several paths this behavior can be seen for some ionospheres but not others. For this set, it appears constrained in H' to between 73 and 78 and either beta between 0.3 and 0.4 or greater than 0.5. A map is included as a reference for these paths.