Merge pull request #43 from fgasdia/txpower

Monotonic splines and transmit power

Project.toml

@@ -1,7 +1,7 @@
 name = "LongwaveModePropagator"
 uuid = "38c6559a-e0a7-48c6-827d-3f02e2d19cec"
 authors = ["Forrest Gasdia <EP-Guy@users.noreply.github.com>"]
-version = "0.2.0"
+version = "0.3.0"
 
 [deps]
 Dates = "ade2ca70-3891-5945-98fb-dc099432e06a"

src/IO.jl

@@ -24,6 +24,9 @@ Type for Wait and Spies (1964) exponential ionosphere profiles defined by [`wait
 The [`electroncollisionfrequency`](@ref) is used for the electron-neutral collision frequency
 profile.
 
+- The electron density profile begins at 40 km altitude and extends to 110 km.
+- The transmitter power is 1 kW.
+
 # Fields
 
 - `name::String`
@@ -343,8 +346,7 @@ Return `HomogeneousWaveguide` from the `i`th entry in each field of `s`.
 """
 function buildwaveguide(s::ExponentialInput, i)
     bfield = BField(s.b_mags[i], s.b_dips[i], s.b_azs[i])
-    species = Species(QE, ME, z -> waitprofile(z, s.hprimes[i], s.betas[i];
-                                               cutoff_low=40e3),
+    species = Species(QE, ME, z -> waitprofile(z, s.hprimes[i], s.betas[i]; cutoff_low=40e3),
                       electroncollisionfrequency)
     ground = Ground(s.ground_epsrs[i], s.ground_sigmas[i])
     return HomogeneousWaveguide(bfield, species, ground, s.segment_ranges[i])
@@ -353,13 +355,19 @@ end
 """
     buildwaveguide(s::TableInput, i)
 
-Return `HomogeneousWaveguide` from the `i`th entry in each field of `s` with a linear
-interpolation over `density` and `collision_frequency`.
+Return `HomogeneousWaveguide` from the `i`th entry in each field of `s` with a
+FritschButland monotonic interpolation over `density` and `collision_frequency`.
+
+Outside of `s.altitude` the nearest `s.density` or `s.collision_frequency` is used.
 """
 function buildwaveguide(s::TableInput, i)
     bfield = BField(s.b_mags[i], s.b_dips[i], s.b_azs[i])
-    density_itp = LinearInterpolation(s.altitude, s.density[i]; extrapolation_bc=Line())
-    collision_itp = LinearInterpolation(s.altitude, s.collision_frequency[i]; extrapolation_bc=Line())
+
+    ditp = interpolate(s.altitude, s.density[i], FritschButlandMonotonicInterpolation())
+    citp = interpolate(s.altitude, s.collision_frequency[i], FritschButlandMonotonicInterpolation())
+
+    density_itp = extrapolate(ditp, Flat())
+    collision_itp = extrapolate(citp, Flat())
     species = Species(QE, ME, density_itp, collision_itp)
     ground = Ground(s.ground_epsrs[i], s.ground_sigmas[i])
     return HomogeneousWaveguide(bfield, species, ground, s.segment_ranges[i])
@@ -371,6 +379,9 @@ end
     buildrun(s::BatchInput; mesh=nothing, unwrap=true, params=LMPParams())
 
 Build LMP structs from an `Input` and run `LMP`.
+
+For `TableInput`s, a FritschButland monotonic interpolation is performed over `density` and
+`collision_frequency`.
 """
 buildrun
 
@@ -384,17 +395,17 @@ function buildrun(s::ExponentialInput; mesh=nothing, unwrap=true, params=LMPPara
         ground = Ground(only(s.ground_epsrs), only(s.ground_sigmas))
         waveguide = HomogeneousWaveguide(bfield, species, ground)
 
-        tx = Transmitter(VerticalDipole(), Frequency(s.frequency), 100e3)
+        tx = Transmitter(s.frequency)
         rx = GroundSampler(s.output_ranges, Fields.Ez)
     else
         # SegmentedWaveguide
         waveguide = SegmentedWaveguide([buildwaveguide(s, i) for i in
                                         eachindex(s.segment_ranges)])
-        tx = Transmitter(VerticalDipole(), Frequency(s.frequency), 100e3)
+        tx = Transmitter(s.frequency)
         rx = GroundSampler(s.output_ranges, Fields.Ez)
     end
 
-    E, amp, phase = propagate(waveguide, tx, rx; mesh=mesh, unwrap=unwrap, params=params)
+    _, amp, phase = propagate(waveguide, tx, rx; mesh=mesh, unwrap=unwrap, params=params)
 
     output = BasicOutput()
     output.name = s.name
@@ -416,23 +427,26 @@ function buildrun(s::TableInput; mesh=nothing, unwrap=true, params=LMPParams())
     if length(s.segment_ranges) == 1
         # HomogeneousWaveguide
         bfield = BField(only(s.b_mags), only(s.b_dips), only(s.b_azs))
-        density_itp = LinearInterpolation(s.altitude, only(s.density); extrapolation_bc=Line())
-        collision_itp = LinearInterpolation(s.altitude, only(s.collision_frequency); extrapolation_bc=Line())
+
+        ditp = interpolate(s.altitude, only(s.density), FritschButlandMonotonicInterpolation())
+        citp = interpolate(s.altitude, only(s.collision_frequency), FritschButlandMonotonicInterpolation())
+
+        density_itp = extrapolate(ditp, Flat())
+        collision_itp = extrapolate(citp, Flat())
         species = Species(QE, ME, density_itp, collision_itp)
         ground = Ground(only(s.ground_epsrs), only(s.ground_sigmas))
         waveguide = HomogeneousWaveguide(bfield, species, ground)
 
-        tx = Transmitter(VerticalDipole(), Frequency(s.frequency), 100e3)
+        tx = Transmitter(s.frequency)
         rx = GroundSampler(s.output_ranges, Fields.Ez)
     else
         # SegmentedWaveguide
-        waveguide = SegmentedWaveguide([buildwaveguide(s, i) for i in
-                                        eachindex(s.segment_ranges)])
-        tx = Transmitter(VerticalDipole(), Frequency(s.frequency), 100e3)
+        waveguide = SegmentedWaveguide([buildwaveguide(s, i) for i in eachindex(s.segment_ranges)])
+        tx = Transmitter(s.frequency)
         rx = GroundSampler(s.output_ranges, Fields.Ez)
     end
 
-    E, amp, phase = propagate(waveguide, tx, rx; mesh=mesh, unwrap=unwrap, params=params)
+    _, amp, phase = propagate(waveguide, tx, rx; mesh=mesh, unwrap=unwrap, params=params)
 
     output = BasicOutput()
     output.name = s.name

src/LongwaveModePropagator.jl

@@ -87,7 +87,7 @@ Parameters for the `LongwaveModePropagator` module with defaults:
 - `approxsusceptibility::Bool = false`: use a cubic interpolating spline representation of
     [`susceptibility`](@ref) during the integration of [`dRdz`](@ref).
 - `susceptibilitysplinestep::Float64 = 10.0`: altitude step in meters used to build the
-    cubic spline representation of [`susceptibility`](@ref) if `approxsusceptibility == true`.
+    spline representation of [`susceptibility`](@ref) if `approxsusceptibility == true`.
 - `grpfparams::GRPFParams = GRPFParams(100000, 1e-5, true)`: parameters for the `GRPF`
     complex root-finding algorithm.
 - `integrationparams::IntegrationParams{T} =

test/IO.jl

@@ -283,6 +283,35 @@ end
 
 function basic_lmp()
     output = LMP.propagate("basic.json")
+
+    # Compare to a `propagate` call with waveguide arguments
+    tx = Transmitter(24e3)
+    rx = GroundSampler(0:20e3:2000e3, Fields.Ez)
+
+    # From `generate_basic()`
+    segment_ranges = [0, 500e3, 1000e3, 1500e3]
+    hprimes = [72.0, 74, 75, 76]
+    betas = [0.28, 0.30, 0.32, 0.35]
+    b_mags = fill(50e-6, length(segment_ranges))
+    b_dips = fill(π/2, length(segment_ranges))
+    b_azs = fill(0.0, length(segment_ranges))
+    ground_sigmas = [0.001, 0.001, 0.0005, 0.0001]
+    ground_epsrs = [4, 4, 10, 10]
+
+    wvg = SegmentedWaveguide(
+        [HomogeneousWaveguide(
+            BField(b_mags[i], b_dips[i], b_azs[i]),
+            Species(QE, ME, z->waitprofile(z, hprimes[i], betas[i]; cutoff_low=40e3), electroncollisionfrequency),
+            Ground(ground_epsrs[i], ground_sigmas[i]),
+            segment_ranges[i]
+        ) for i = 1:4]
+    )
+
+    _, a, p = propagate(wvg, tx, rx)
+
+    @test isapprox(a, output.amplitude, atol=0.1)
+    @test isapprox(p, output.phase, atol=deg2rad(1))
+
     return output
 end
 

test/LongwaveModePropagator.jl

@@ -52,7 +52,6 @@ function test_propagate_segmented(scenario)
     @test eltype(amp) == Float64
     @test eltype(phase) == Float64
 
-    mesh = LMP.defaultmesh(tx.frequency)
     E1, amp1, phase1 = propagate(waveguide, tx, rx)
     @test E1 ≈ E    rtol=1e-3
     @test amp1 ≈ amp    rtol=1e-3
@@ -76,11 +75,11 @@ function test_mcranges_segmented(scenario)
     waveguide = SegmentedWaveguide([HomogeneousWaveguide(bfield[i], species[i], ground[i],
                                                          distances[i]) for i in 1:2])
 
-    Eref, ampref, phaseref = propagate(waveguide, tx, rx; unwrap=false)
+    _, ampref, phaseref = propagate(waveguide, tx, rx; unwrap=false)
 
-    E1, a1, p1 = propagate(waveguide, tx, GroundSampler(600e3, Fields.Ez))
-    E2, a2, p2 = propagate(waveguide, tx, GroundSampler(1400e3, Fields.Ez))
-    E3, a3, p3 = propagate(waveguide, tx, GroundSampler(1800e3, Fields.Ez))
+    _, a1, p1 = propagate(waveguide, tx, GroundSampler(600e3, Fields.Ez))
+    _, a2, p2 = propagate(waveguide, tx, GroundSampler(1400e3, Fields.Ez))
+    _, a3, p3 = propagate(waveguide, tx, GroundSampler(1800e3, Fields.Ez))
 
     m1 = findfirst(isequal(600e3), rx.distance)
     m2 = findfirst(isequal(1400e3), rx.distance)

test/lwpc_comparisons.jl

@@ -48,6 +48,7 @@ end
         resonant_horizontal_scenario=>"resonant_horizontal.log",
         segmented_scenario=>"segmented.log")
 
+    @info "  Running:"
     for scn in (verticalB_scenario, resonant_scenario, nonresonant_scenario,
         resonant_elevatedrx_scenario, resonant_horizontal_scenario,)
 

test/magnetoionic.jl

@@ -1,10 +1,9 @@
 function test_susceptibility(scenario)
-    @unpack ea, tx, bfield, species, ground = scenario
+    @unpack tx, bfield, species, ground = scenario
 
     M1 = LMP.susceptibility(70e3, tx.frequency, bfield, species)
     M2 = LMP.susceptibility(70e3, tx.frequency, bfield, species; params=LMPParams())
-    M3 = LMP.susceptibility(70e3, tx.frequency, bfield, species;
-                            params=LMPParams(earthradius=6350e3))
+    M3 = LMP.susceptibility(70e3, tx.frequency, bfield, species; params=LMPParams(earthradius=6350e3))
     @test M1 == M2
     @test !(M2 ≈ M3)
 
@@ -13,13 +12,11 @@ function test_susceptibility(scenario)
 
     M4 = LMP.susceptibility(70e3, tx.frequency, waveguide)
     M5 = LMP.susceptibility(70e3, tx.frequency, waveguide; params=LMPParams())
-    M6 = LMP.susceptibility(70e3, tx.frequency, waveguide;
-                            params=LMPParams(earthradius=6350e3))
+    M6 = LMP.susceptibility(70e3, tx.frequency, waveguide; params=LMPParams(earthradius=6350e3))
 
     M7 = LMP.susceptibility(70e3, modeequation)
     M8 = LMP.susceptibility(70e3, modeequation; params=LMPParams())
-    M9 = LMP.susceptibility(70e3, modeequation;
-                            params=LMPParams(earthradius=6350e3))
+    M9 = LMP.susceptibility(70e3, modeequation; params=LMPParams(earthradius=6350e3))
 
     @test M4 == M5 == M1
     @test M6 == M3
@@ -31,7 +28,7 @@ function test_susceptibility(scenario)
 end
 
 function test_spline(scenario)
-    @unpack ea, tx, bfield, species, ground = scenario
+    @unpack tx, bfield, species = scenario
 
     itp = LMP.susceptibilityspline(tx.frequency, bfield, species)
 

test/modeconversion.jl

@@ -1,5 +1,5 @@
 function test_modeconversion_segmented(scenario)
-    @unpack distances, ea, tx, bfield, species, ground = scenario()
+    @unpack distances, tx, bfield, species, ground = scenario()
     params = LMPParams()
 
     waveguide = SegmentedWaveguide([HomogeneousWaveguide(bfield[i], species[i],

test/modefinder.jl

@@ -102,7 +102,7 @@ function test_dRdθdz(scenario)
 end
 
 function test_integratedreflection_vertical(scenario)
-    @unpack ea, tx, ground, bfield, species = scenario
+    @unpack tx, ground, bfield, species = scenario
 
     waveguide = HomogeneousWaveguide(bfield, species, ground)
     me = PhysicalModeEquation(tx.frequency, waveguide)

test/modesum.jl

@@ -32,7 +32,6 @@ function test_Efield(scenario)
     modes = TEST_MODES[scenario]
 
     X = LMP.distance(rx, tx)
-    E = zeros(ComplexF64, length(X))
 
     E1 = LMP.Efield(modes, waveguide, tx, rx)  # out-of-place
 

test/runtests.jl

@@ -1,5 +1,4 @@
-using Test, Dates, Random
-using LinearAlgebra, Statistics
+using Test, Dates, Random, LinearAlgebra, Statistics
 using StaticArrays
 using Parameters
 using OrdinaryDiffEq, DiffEqCallbacks

test/wavefields.jl

@@ -9,7 +9,7 @@ function randomwavefields()
 end
 
 function test_Wavefields()
-    modes, wavefields = randomwavefields()
+    _, wavefields = randomwavefields()
     @unpack wavefieldheights = LMPParams()
 
     @test LMP.numheights(wavefields) == length(wavefieldheights)
@@ -155,7 +155,7 @@ function test_boundaryscalars(scenario)
 end
 
 function test_fieldstrengths(scenario)
-    @unpack ea, bfield, tx, ground, species = scenario
+    @unpack bfield, tx, ground, species = scenario
     modes, wavefields = randomwavefields()
 
     modes = LMP.eigenangles(wavefields)