Graham McIntire
49cbe6789c
- Replace piecewise knife-edge loss with P.526-16 Eq. 31 single formula - Fix diffraction parameter ν to use standard formula instead of ad-hoc approximation - Implement Deygout 3-edge method for multiple obstacle diffraction - Add dynamic k-factor from HRRR refractivity gradient (falls back to 4/3) - Terrain worker now looks up nearest HRRR profile for atmospheric correction - Update algo.md with P.526-16 methods and k-factor table - Fix pre-existing map_live_test antenna height default (33 ft, not 8)
282 lines
8.9 KiB
Elixir
282 lines
8.9 KiB
Elixir
defmodule Microwaveprop.Terrain.TerrainAnalysisTest do
|
||
use ExUnit.Case, async: true
|
||
|
||
alias Microwaveprop.Terrain.TerrainAnalysis
|
||
|
||
describe "fresnel_radius/3" do
|
||
test "returns 0 when d1 or d2 is 0" do
|
||
assert TerrainAnalysis.fresnel_radius(0, 1000, 0.023) == 0
|
||
assert TerrainAnalysis.fresnel_radius(1000, 0, 0.023) == 0
|
||
end
|
||
|
||
test "computes correct Fresnel radius for known values" do
|
||
# 1296 MHz -> lambda = 0.2315 m
|
||
lambda = 0.3 / 1.296
|
||
d1 = 50_000.0
|
||
d2 = 50_000.0
|
||
|
||
r = TerrainAnalysis.fresnel_radius(d1, d2, lambda)
|
||
assert_in_delta r, 76.1, 1.0
|
||
end
|
||
end
|
||
|
||
describe "earth_bulge/3" do
|
||
test "returns 0 at endpoints" do
|
||
assert TerrainAnalysis.earth_bulge(0.0, 100.0) == 0.0
|
||
assert TerrainAnalysis.earth_bulge(1.0, 100.0) == 0.0
|
||
end
|
||
|
||
test "maximum bulge at midpoint of 100 km path" do
|
||
bulge = TerrainAnalysis.earth_bulge(0.5, 100.0)
|
||
assert_in_delta bulge, 147.0, 2.0
|
||
end
|
||
|
||
test "bulge increases with path length" do
|
||
bulge_50 = TerrainAnalysis.earth_bulge(0.5, 50.0)
|
||
bulge_100 = TerrainAnalysis.earth_bulge(0.5, 100.0)
|
||
assert bulge_100 > bulge_50
|
||
end
|
||
|
||
test "bulge decreases with larger k-factor (super-refraction)" do
|
||
bulge_std = TerrainAnalysis.earth_bulge(0.5, 100.0, 4 / 3)
|
||
bulge_super = TerrainAnalysis.earth_bulge(0.5, 100.0, 2.0)
|
||
assert bulge_super < bulge_std
|
||
end
|
||
end
|
||
|
||
# P.526-16 Eq. 31: J(ν) = 6.9 + 20*log10(sqrt((ν-0.1)² + 1) + ν - 0.1)
|
||
describe "knife_edge_loss/1 (P.526-16 Eq. 31)" do
|
||
test "returns 0 for ν ≤ -0.78 (well clear)" do
|
||
assert TerrainAnalysis.knife_edge_loss(-1.0) == 0
|
||
assert TerrainAnalysis.knife_edge_loss(-0.8) == 0
|
||
end
|
||
|
||
test "returns ~6 dB at grazing (ν = 0)" do
|
||
loss = TerrainAnalysis.knife_edge_loss(0.0)
|
||
assert_in_delta loss, 6.0, 0.1
|
||
end
|
||
|
||
test "returns ~13.9 dB at ν = 1" do
|
||
loss = TerrainAnalysis.knife_edge_loss(1.0)
|
||
assert_in_delta loss, 13.9, 0.2
|
||
end
|
||
|
||
test "returns ~19.0 dB at ν = 2" do
|
||
loss = TerrainAnalysis.knife_edge_loss(2.0)
|
||
assert_in_delta loss, 19.0, 0.2
|
||
end
|
||
|
||
test "returns ~22.5 dB at ν = 3" do
|
||
loss = TerrainAnalysis.knife_edge_loss(3.0)
|
||
assert_in_delta loss, 22.5, 0.3
|
||
end
|
||
|
||
test "monotonically increases with ν" do
|
||
values = [-0.5, 0.0, 0.5, 1.0, 1.5, 2.0, 3.0, 5.0]
|
||
|
||
losses = Enum.map(values, &TerrainAnalysis.knife_edge_loss/1)
|
||
|
||
losses
|
||
|> Enum.chunk_every(2, 1, :discard)
|
||
|> Enum.each(fn [a, b] -> assert b >= a end)
|
||
end
|
||
end
|
||
|
||
describe "diffraction_param/4" do
|
||
test "returns positive ν for obstacle above beam (blocked)" do
|
||
# h = 50m above beam, d1 = d2 = 25km, 10 GHz (λ = 0.03m)
|
||
nu = TerrainAnalysis.diffraction_param(50.0, 25_000.0, 25_000.0, 0.03)
|
||
assert nu > 0
|
||
end
|
||
|
||
test "returns negative ν for obstacle below beam (clear)" do
|
||
nu = TerrainAnalysis.diffraction_param(-50.0, 25_000.0, 25_000.0, 0.03)
|
||
assert nu < 0
|
||
end
|
||
|
||
test "returns 0 at grazing" do
|
||
assert TerrainAnalysis.diffraction_param(0.0, 25_000.0, 25_000.0, 0.03) == 0.0
|
||
end
|
||
|
||
test "ν increases with frequency (shorter wavelength = sharper shadow)" do
|
||
nu_10ghz = TerrainAnalysis.diffraction_param(50.0, 25_000.0, 25_000.0, 0.03)
|
||
nu_24ghz = TerrainAnalysis.diffraction_param(50.0, 25_000.0, 25_000.0, 0.0125)
|
||
assert nu_24ghz > nu_10ghz
|
||
end
|
||
end
|
||
|
||
describe "k_factor/1" do
|
||
test "returns 4/3 for standard atmosphere (-39 N/km)" do
|
||
k = TerrainAnalysis.k_factor(-39.0)
|
||
assert_in_delta k, 4 / 3, 0.01
|
||
end
|
||
|
||
test "returns ~1.0 for no refraction (0 N/km)" do
|
||
k = TerrainAnalysis.k_factor(0.0)
|
||
assert_in_delta k, 1.0, 0.01
|
||
end
|
||
|
||
test "returns large value approaching ducting (-157 N/km)" do
|
||
k = TerrainAnalysis.k_factor(-150.0)
|
||
assert k > 5
|
||
end
|
||
|
||
test "returns 4/3 for nil input" do
|
||
assert_in_delta TerrainAnalysis.k_factor(nil), 4 / 3, 0.001
|
||
end
|
||
end
|
||
|
||
describe "analyse/4" do
|
||
test "returns CLEAR for flat terrain with antenna heights" do
|
||
profile =
|
||
for i <- 0..10 do
|
||
f = i / 10
|
||
%{lat: 32.9 + f * 0.09, lon: -97.0, d: f, elev: 0.0, dist_km: f * 10.0}
|
||
end
|
||
|
||
result = TerrainAnalysis.analyse(profile, 10.0, 1.296, ant_ht_a: 30.0, ant_ht_b: 30.0)
|
||
|
||
assert result.verdict == "CLEAR"
|
||
assert result.obstructed_count == 0
|
||
assert result.diffraction_db == 0
|
||
end
|
||
|
||
test "returns BLOCKED for high obstacle" do
|
||
profile =
|
||
for i <- 0..10 do
|
||
f = i / 10
|
||
elev = if i == 5, do: 500.0, else: 100.0
|
||
%{lat: 32.9 + f, lon: -97.0, d: f, elev: elev, dist_km: f * 50.0}
|
||
end
|
||
|
||
result = TerrainAnalysis.analyse(profile, 50.0, 1.296)
|
||
|
||
assert result.verdict == "BLOCKED"
|
||
assert result.diffraction_db > 0
|
||
end
|
||
|
||
test "max_elevation_m reflects the terrain peak" do
|
||
profile =
|
||
for i <- 0..4 do
|
||
f = i / 4
|
||
elev = if i == 2, do: 500.0, else: 100.0
|
||
%{lat: 32.9 + f, lon: -97.0, d: f, elev: elev, dist_km: f * 50.0}
|
||
end
|
||
|
||
result = TerrainAnalysis.analyse(profile, 50.0, 1.296)
|
||
assert_in_delta result.max_elevation_m, 500.0, 0.1
|
||
end
|
||
|
||
test "elevated endpoints clear over flat terrain" do
|
||
profile =
|
||
for i <- 0..10 do
|
||
f = i / 10
|
||
%{lat: 32.9 + f * 0.18, lon: -97.0, d: f, elev: 0.0, dist_km: f * 20.0}
|
||
end
|
||
|
||
result = TerrainAnalysis.analyse(profile, 20.0, 1.296, ant_ht_a: 200.0, ant_ht_b: 200.0)
|
||
assert result.verdict == "CLEAR"
|
||
end
|
||
|
||
test "returns FRESNEL verdict for moderate ridge" do
|
||
profile =
|
||
for i <- 0..10 do
|
||
f = i / 10
|
||
elev = if i == 5, do: 165.0, else: 0.0
|
||
%{lat: 32.9 + f * 0.18, lon: -97.0, d: f, elev: elev, dist_km: f * 20.0}
|
||
end
|
||
|
||
result = TerrainAnalysis.analyse(profile, 20.0, 1.296, ant_ht_a: 200.0, ant_ht_b: 200.0)
|
||
assert result.verdict in ["FRESNEL_MINOR", "FRESNEL_PARTIAL"]
|
||
assert result.fresnel_hit_count >= 1
|
||
end
|
||
|
||
test "antenna heights raise beam above terrain" do
|
||
profile =
|
||
for i <- 0..10 do
|
||
f = i / 10
|
||
elev = if i == 5, do: 30.0, else: 0.0
|
||
%{lat: 32.9 + f * 0.09, lon: -97.0, d: f, elev: elev, dist_km: f * 10.0}
|
||
end
|
||
|
||
result_low = TerrainAnalysis.analyse(profile, 10.0, 1.296)
|
||
result_high = TerrainAnalysis.analyse(profile, 10.0, 1.296, ant_ht_a: 100.0, ant_ht_b: 100.0)
|
||
|
||
assert result_low.verdict == "BLOCKED"
|
||
assert result_high.verdict == "CLEAR"
|
||
end
|
||
|
||
test "min_clearance_m can be negative for obstructed paths" do
|
||
profile =
|
||
for i <- 0..4 do
|
||
f = i / 4
|
||
elev = if i == 2, do: 500.0, else: 100.0
|
||
%{lat: 32.9 + f, lon: -97.0, d: f, elev: elev, dist_km: f * 50.0}
|
||
end
|
||
|
||
result = TerrainAnalysis.analyse(profile, 50.0, 1.296)
|
||
assert result.min_clearance_m < 0
|
||
end
|
||
|
||
test "accepts k_factor option for atmospheric correction" do
|
||
profile =
|
||
for i <- 0..10 do
|
||
f = i / 10
|
||
elev = if i == 5, do: 100.0, else: 0.0
|
||
%{lat: 32.9 + f * 0.09, lon: -97.0, d: f, elev: elev, dist_km: f * 50.0}
|
||
end
|
||
|
||
# Super-refraction (k=2) reduces earth bulge → less total obstruction → less loss
|
||
result_std = TerrainAnalysis.analyse(profile, 50.0, 10.0, k_factor: 4 / 3)
|
||
result_super = TerrainAnalysis.analyse(profile, 50.0, 10.0, k_factor: 2.0)
|
||
|
||
assert result_super.diffraction_db <= result_std.diffraction_db
|
||
end
|
||
|
||
test "Deygout: two obstacles produce more loss than single worst" do
|
||
# Two 300m peaks at 1/3 and 2/3 of a 60km path, endpoints at 0m
|
||
profile =
|
||
for i <- 0..6 do
|
||
f = i / 6
|
||
|
||
elev =
|
||
case i do
|
||
2 -> 300.0
|
||
4 -> 300.0
|
||
_ -> 0.0
|
||
end
|
||
|
||
%{lat: 32.9 + f, lon: -97.0, d: f, elev: elev, dist_km: f * 60.0}
|
||
end
|
||
|
||
result = TerrainAnalysis.analyse(profile, 60.0, 10.0)
|
||
|
||
# Single 300m peak at midpoint of 60km path for comparison
|
||
single_profile =
|
||
for i <- 0..6 do
|
||
f = i / 6
|
||
elev = if i == 3, do: 300.0, else: 0.0
|
||
%{lat: 32.9 + f, lon: -97.0, d: f, elev: elev, dist_km: f * 60.0}
|
||
end
|
||
|
||
single_result = TerrainAnalysis.analyse(single_profile, 60.0, 10.0)
|
||
|
||
# Two obstacles should produce MORE diffraction loss than a single one
|
||
assert result.diffraction_db > single_result.diffraction_db
|
||
end
|
||
|
||
test "higher frequency produces more diffraction loss for same obstacle" do
|
||
profile =
|
||
for i <- 0..10 do
|
||
f = i / 10
|
||
elev = if i == 5, do: 200.0, else: 0.0
|
||
%{lat: 32.9 + f, lon: -97.0, d: f, elev: elev, dist_km: f * 30.0}
|
||
end
|
||
|
||
result_10ghz = TerrainAnalysis.analyse(profile, 30.0, 10.0)
|
||
result_24ghz = TerrainAnalysis.analyse(profile, 30.0, 24.0)
|
||
|
||
assert result_24ghz.diffraction_db > result_10ghz.diffraction_db
|
||
end
|
||
end
|
||
end
|