fgasdia / longwavemodepropagator.jl Goto Github PK

#55 @JuliaRegistrator register

The default tolerance for the dR/dz solver (LMPParams().integrationparams.tolerance) is currently 1e-8. However, the tolerance used by the GRPF algorithm (currently default at 1e-5) dominates the affect of the differential equation solver tolerance when it comes to identifying eigenangles. In fact, even an integration tolerance of 1e-3 finds the same modes as a tolerance of 1e-10, but is more than a factor of 2 faster. Sometimes 1e-3 results in instabilities, so a default tolerance of 1e-4 or 1e-5 should be used instead.

Additionally, findmodes currently performs filtering of the candidate modes by 1) uniqueness of the modes to some tolerance and 2) closeness of the mode condition to 0. However, it appears that there is rarely (never?) a need to remove any modes based on 2 (GRPF takes care of this correctly). In fact, it appears that valid modes are unnecessarily removed for nearly any waveguide scenario with the current settings. This filtering should be removed.

In a test of 30 random waveguides there were not any duplicate modes either, but I know in earlier versions of LMP this occurred. To some extent, it is down to the user to ensure the tolerances on GRPF and integrationparams and appropriate (if the tolerance of the mode finder is too small relative to the integration tolerance, identical modes may be identified).

This is breaking because propagation results will change.

It would be convenient to read the .inp and .ndx files.

Magnetic dip angles of +/- 1 deg (horizontal field) are leading to numerical issues.

Using the GroundSampler to set up equivalent samplers for Ex, Ey, and Ez, for a variety of ionospheres and grounds, Ex dominates several orders of magnitude above Ez. Ey behaves predictably with much lower amplitude and Ez behaves appropriately when compared to Lehtninen's Full Wave Method. Code used to create the figure is included below.

    using LongwaveModePropagator
    using LongwaveModePropagator: QE, ME
    using Plots
    using LaTeXStrings
    using Colors, ColorSchemes

    const LMP = LongwaveModePropagator


    txpower = 1e6
    # "standard" vertical dipole transmitter at 24 kHz
    tx = Transmitter(VerticalDipole(),24e3,txpower)

    # sample vertical electric field every 5 km out to 2000 km from tx
    rx_z = GroundSampler(0:1e3:1500e3, Fields.Ez)
    rx_y = GroundSampler(0:1e3:1500e3, Fields.Ey)
    rx_x = GroundSampler(0:1e3:1500e3, Fields.Ex)

    # vertical magnetic field
    bfield = BField(50e-6, π/2, 0)

    # "typical" earth ground
    ground = Ground(10, 1e-4)
    h=75
    electrons = Species(QE, ME, z->waitprofile(z, h, 0.4), electroncollisionfrequency)

    waveguide = HomogeneousWaveguide(bfield, electrons, ground)

    # return the complex electric field, amplitude (dB uV/m), and phase (radians)
    Ez, az, pz = propagate(waveguide, tx, rx_z)
    Ey, ay, py = propagate(waveguide, tx, rx_y)
    Ex, ax, px = propagate(waveguide, tx, rx_x)
    
    # transform units to V/m
    az_V = 10 .^((az.-120) ./20)
    ay_V = 10 .^((ay.-120) ./20)
    ax_V = 10 .^((ax.-120) ./20)
   
#----------------------------------------------------------------------------
    #Plot

    plot(rx_x.distance/1000, [ax_V ay_V az_V],label=["x" "y" "z"], xlabel="distance (km)", ylabel="amplitude (E V/m)",title="E, h' = 75, \$ \\beta \$ = 0.4, Ground",dpi=400,ylims = (1e-8, 1e4),xlims = (0,1500),yscale=:log10,minorgrid=true)

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.

Consider adding IO file support for multiple field components.

Currently the Modified Hankel Functions of Order One-Third are used to compute the height gain functions (and excitation factors) for the mode sum. Consider integrating the wavefields directly. This may be easier to extend fields into the ionosphere and above.

As part of this, study the accuracy of the Modified Hankel Functions relative to the full wave solution.

The increase in runtime will be negligible relative to mode finding.

The machinery currently exists to automatically calculate the distance between a Transmitter and Receiver, but the geodesic calculation would either require Proj4.jl (a wrapper for the PROJ C library) or GeographicLib.jl, which is pure Julia. I'd prefer to use GeographicLib, but it must be registered to be a dependency of this package (which is registered). Waiting on anowacki/GeographicLib.jl#2.

The current workaround is to use a GroundSampler or Sampler with a single pre-computed distance from the Transmitter.

@JuliaRegistrator register

This could potentially be implemented by interpolating the location of eigenangles as a function of frequency. It would be especially useful if we could specify the frequency content of the transmitter. See the papers on "pulse" programs from NOSC.

@JuliaRegistrator register

There is a specialized form for StaticMatrix types.

The default finemesh range is new

BasicInput is used for Wait density profile - maybe call it ExponentialInput or something similar.

Look at improving the original MODESRCH algorithm which has robustness problems but is much faster than GRPF.

By default, a ProgressMeter prints when propagating a SegmentedWaveguide.

Problems:

• There is a higher level ProgressMeter printed when running a collection of SegmentedWaveguides, but it is not shown until the entire collection is complete.
• In fact, ProgressMeters cannot currently be "stacked" with multiple levels of progress, and even if they could be, the user code would need to interact with LongwaveModePropagator's ProgressMeter objects.
• When estimating with LongwaveModePropagator, there will be many, many progress bars printed. This crowds the console and the user cannot use their own ProgressMeter to monitor overall progress.

The ProgressMeter should probably be replaced with Logging utilities. See e.g. MicroLogging.jl, LoggingExtras.jl, TerminalLoggers.jl, or ProgressLogging.jl.

I'm slightly confused by the way the dev version of docs is rendering. Literate.jl may be adding text that is printed (like nothing #hide shortly after the text line "Then we simply define the ODEProblem and solve" in Solvers for ionospheric reflection coefficient. There are also figures not being displayed, even though they are in the stable version.

When setting up a GroundSampler, if a single (or set of few) distance(s) is provided (at range for a rx for instance), the phase returned by "propagate()" may be ambiguous by multiples of 2pi radians.

Code used to generate this figure:

using LongwaveModePropagator
using LongwaveModePropagator: QE, ME
using Plots
using LaTeXStrings


# "standard" vertical dipole transmitter at 24 kHz 1MW
txpower = 1e6
tx = Transmitter(VerticalDipole(),24e3,txpower)

# sample vertical electric field every 1 km, 250km, and 500km, out to 1500 km from tx
rx_z = GroundSampler(0:1e3:1500e3, Fields.Ez)
rx_250 = GroundSampler([250e3,500e3,750e3, 1000e3,1250e3,1500e3],Fields.Ez)
rx_500 = GroundSampler([500e3,1000e3,1500e3],Fields.Ez)

# vertical magnetic field
bfield = BField(50e-6, π/2, 0)

# "typical" seawater ground
ground = Ground(81, 4)

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

waveguide = HomogeneousWaveguide(bfield, electrons, ground)

# return the complex electric field, amplitude (dB uV/m), and phase (radians)
Ez, az, pz = propagate(waveguide, tx, rx_z);
E, a, p_250 = propagate(waveguide, tx, rx_250);
E, a, p_500 = propagate(waveguide, tx, rx_500);

#----------------------------------------------------------------------------

plt=plot(rx_z.distance/1000,pz.*(180/pi), ylims=(-90,450),label="1km Sampler",dpi=400)
scatter!(rx_250.distance/1000,p_250.*(180/pi),label="250km Sampler")
scatter!(rx_500.distance/1000,p_500.*(180/pi),label="500km Sampler")
plot!(xlabel="Distance (km)",ylabel="Phase(\$ \\degree \$)")

savefig(plt,"PhaseTestPlot.png")

See https://github.com/github/renaming

expected January 2021

It would be nice to allow creation with syntax like this:

wvg = SegmentedWaveguide() do w
        for i in eachindex(hps)
            species = Species(QE, ME, z->waitprofile(z, hps[i], βs[i], cutoff_low=40e3), electroncollisionfrequency)
            w = HomogeneousWaveguide(BFIELD, species, GND, dists[i])
        end
end

EigenAngle currently precomputes and stores trigonometric functions of ea. Computation of these functions may be faster than memory access of the struct fields. See e.g. https://viralinstruction.com/posts/optimise/#reduce_memory_consumption.

This needs to be specified on read and write.

Although accepting NaN is supported by Python, the Matlab built-in JSON library does not support NaN. Instead, Matlab uses null for NaN, but we can't specify null because Julia interprets that as nothing.

Also see PropagationModelPrep.

@JuliaRegistrator register

Release notes:

• Radiation resistance correction for transmitting antennas elevated above the ground
  • This can be toggled on by setting radiationresistancecorrection = true in LMPParams(). By default, radiationresistancecorrection = false until the implementation can be quantitatively verified against another model. The radiation resistance correction assumes a perfectly conductive ground, rather than the realistic Ground defined by the user at the transmitter.
• The returned electric field E is no longer set to a value matching LWPC if the x distance from the transmitter is less than 1 meter. Instead, if the distance is 0 meters from the transmitter, than the returned field will be Inf.
• compat version bumps, including minimum Julia version of 1.10

Dear Forrest,

I exchanged a couple of emails with Bob Marshall with questions related to LMP. Bob contacted you but did not get any answer. Bob could not answer my questions and tried to contact you, but failed. He said your CU email([email protected]) might not work.
I hope you can get this message through GitHub and can answer the questions and help us. I attach the mails.

Thank you,

Janos
lmp.zip

In 29ed585 the length of LMPParams().wavefieldheights was increased to match the small integration step size needed to replicate the wavefields from Pitteway 1965 (documented here). Was this really necessary for the integration methods being used, or is it only something needed for plotting?

Add a Fields.E which returns Ex, Ey, and Ez field components at once.

When using LinearInterpolation for numberdensity and collisionfrequency fields of Species, the total runtime for a "Wait profile" is approximately twice the runtime for the same discretely sampled Wait profile with a CubicSplineInterpolation. This is likely because the DifferentialEquations solver is trying to resolve the "knees" between each piecewise linear segment.

Cubic splines may overshoot the local true density, so PCHIP interpolation is preferred, but there is no up-to-date Julia package for PCHIP right now. I have a private GitHub repo with three PCHIP algorithms but they can't be released under MIT license, so I'd rather not register them.

@JuliaRegistrator register

Release notes:

• LongwaveModePropagator.defaultmesh now has different default behavior for transmitter frequencies below 12 kHz.
  • The previous default used the internal trianglemesh function and searched for eigenangles in the region of the lower right complex quadrant near the positive real and negative imaginary axes. At low frequencies the eigenangles can have large real and large negative imaginary components, which were not contained in the trianglemesh. Now RootsAndPoles.rectanglemesh is used for these lower frequencies.
  • Fixed #68

@JuliaRegistrator register

• Bump Interpolations.jl version and remove deprecated functions.

@JuliaRegistrator register

This issue is used to trigger TagBot; feel free to unsubscribe.

If you haven't already, you should update your TagBot.yml to include issue comment triggers.
Please see this post on Discourse for instructions and more details.

If you'd like for me to do this for you, comment TagBot fix on this issue.
I'll open a PR within a few hours, please be patient!

The appendix of Jacobson 2009 (10.1029/2008JA013642) has several tests that are relevant to portions of this code. Consider implementing them as formal tests.

I suggest pointing out https://github.com/fgasdia/LMPTools.jl in the README in case anyone's looking for ground conductivity, magnetic field models, etc. However, LMPTools is not registered at this time.

@JuliaRegistrator register

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")

One way to improve performance may be to precompute/build an interpolation object for the susceptibility tensor as a function of height. This becomes more efficient relative to computing susceptibility live during integration as the number of species increases. The susceptibility tensor is a function of height only (not eigenangle), so as we solve the mode equation for many different eigenangles we are computing the same susceptibilities over and over (but each call is <200 ns for one species).

There may also be instances where interpolation is less efficient (low frequency, so not many integrations are required & simple profile, so integration is fast). There is some overhead to computing the interpolator.