graham/prop
1
0
Fork 0
Code Issues Pull requests Projects Releases Packages Wiki Activity Actions

test(property): add property tests for terrain + atmospheric math

Browse source 
Cover the pure-math invariants of the ITU-R P.526-16 diffraction
helpers, the haversine metric on Geo, and the Magnus / dry-adiabat
helpers in the skew-T renderer. 29 properties, one invariant each.
This commit is contained in:
Graham McIntire 2026-04-21 13:56:39 -05:00
parent 71e8f53142
commit 4c1013cdbd
No known key found for this signature in database
GPG key ID: F4ABF488E6029E59
3 changed files with 467 additions and 0 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

97
test/microwaveprop/geo_property_test.exs Normal file
View file

@ -0,0 +1,97 @@
defmodule Microwaveprop.GeoPropertyTest do
@moduledoc """
StreamData property tests for `Microwaveprop.Geo`.
The haversine formula defines a metric on the surface of the Earth,
so it must satisfy the three metric axioms: non-negativity with
identity-of-indiscernibles, symmetry, and the triangle inequality.
These properties encode those axioms directly, plus the global
bound of half the Earth's circumference.
"""
use ExUnit.Case, async: true
use ExUnitProperties
alias Microwaveprop.Geo
@half_circumference_km 20_015.1
describe "haversine_km/4 metric axioms" do
property "distance from a point to itself is exactly 0" do
check all(
lat <- float(min: -89.9, max: 89.9),
lon <- float(min: -179.9, max: 179.9)
) do
assert Geo.haversine_km(lat, lon, lat, lon) == 0.0
end
end
property "is non-negative for every pair of points" do
check all(
lat1 <- float(min: -89.9, max: 89.9),
lon1 <- float(min: -179.9, max: 179.9),
lat2 <- float(min: -89.9, max: 89.9),
lon2 <- float(min: -179.9, max: 179.9)
) do
assert Geo.haversine_km(lat1, lon1, lat2, lon2) >= 0
end
end
property "is symmetric in its endpoints" do
check all(
lat1 <- float(min: -89.9, max: 89.9),
lon1 <- float(min: -179.9, max: 179.9),
lat2 <- float(min: -89.9, max: 89.9),
lon2 <- float(min: -179.9, max: 179.9)
) do
forward = Geo.haversine_km(lat1, lon1, lat2, lon2)
reverse = Geo.haversine_km(lat2, lon2, lat1, lon1)
assert_in_delta forward, reverse, 1.0e-9
end
end
property "satisfies the triangle inequality for any three points" do
check all(
lat1 <- float(min: -89.9, max: 89.9),
lon1 <- float(min: -179.9, max: 179.9),
lat2 <- float(min: -89.9, max: 89.9),
lon2 <- float(min: -179.9, max: 179.9),
lat3 <- float(min: -89.9, max: 89.9),
lon3 <- float(min: -179.9, max: 179.9)
) do
d_ab = Geo.haversine_km(lat1, lon1, lat2, lon2)
d_bc = Geo.haversine_km(lat2, lon2, lat3, lon3)
d_ac = Geo.haversine_km(lat1, lon1, lat3, lon3)
# Small epsilon for accumulated floating-point error.
assert d_ac <= d_ab + d_bc + 1.0e-6
end
end
property "is bounded above by half the Earth's circumference" do
check all(
lat1 <- float(min: -89.9, max: 89.9),
lon1 <- float(min: -179.9, max: 179.9),
lat2 <- float(min: -89.9, max: 89.9),
lon2 <- float(min: -179.9, max: 179.9)
) do
assert Geo.haversine_km(lat1, lon1, lat2, lon2) <= @half_circumference_km
end
end
end
describe "bearing_deg/4" do
property "result is always in [0, 360) for distinct points" do
check all(
lat1 <- float(min: -80.0, max: 80.0),
lon1 <- float(min: -179.0, max: 179.0),
# Ensure a nonzero displacement.
dlat <- float(min: 0.1, max: 5.0),
dlon <- float(min: 0.1, max: 5.0)
) do
b = Geo.bearing_deg(lat1, lon1, lat1 + dlat, lon1 + dlon)
assert b >= 0.0
assert b < 360.0
end
end
end
end

251
test/microwaveprop/terrain/terrain_analysis_property_test.exs Normal file
View file

@ -0,0 +1,251 @@
defmodule Microwaveprop.Terrain.TerrainAnalysisPropertyTest do
@moduledoc """
StreamData property tests for ITU-R P.526-16 diffraction math.
Each property encodes a single invariant of the pure-math core of
`Microwaveprop.Terrain.TerrainAnalysis`. The scenarios are
constructed so the expected bound is physically meaningful — e.g.
knife-edge loss is never negative regardless of ν, diffraction loss
vanishes on an unobstructed path, fresnel radius is symmetric in its
distance arguments, and the P.526 piecewise threshold at ν = −0.78
is a true zero floor.
"""
use ExUnit.Case, async: true
use ExUnitProperties
alias Microwaveprop.Terrain.TerrainAnalysis
describe "knife_edge_loss/1" do
property "is non-negative for every real ν" do
check all(nu <- float(min: -50.0, max: 50.0)) do
assert TerrainAnalysis.knife_edge_loss(nu) >= 0
end
end
property "is exactly 0 for ν ≤ −0.78 (P.526 piecewise branch)" do
check all(nu <- float(min: -50.0, max: -0.78)) do
assert TerrainAnalysis.knife_edge_loss(nu) == 0
end
end
property "is monotonically non-decreasing in ν above the cutoff" do
check all(
a <- float(min: -0.7, max: 20.0),
delta <- float(min: 0.0001, max: 5.0)
) do
loss_a = TerrainAnalysis.knife_edge_loss(a)
loss_b = TerrainAnalysis.knife_edge_loss(a + delta)
assert loss_b >= loss_a
end
end
end
describe "diffraction_param/4" do
property "sign of ν matches sign of h when distances and λ are positive" do
check all(
h <- float(min: -500.0, max: 500.0),
d1 <- float(min: 1.0, max: 100_000.0),
d2 <- float(min: 1.0, max: 100_000.0),
lambda <- float(min: 0.001, max: 1.0)
) do
nu = TerrainAnalysis.diffraction_param(h, d1, d2, lambda)
cond do
h > 0 -> assert nu > 0
h < 0 -> assert nu < 0
true -> assert nu == 0.0
end
end
end
property "degenerate inputs (d ≤ 0 or λ ≤ 0) return 0.0" do
check all(
h <- float(min: -500.0, max: 500.0),
bad <- float(min: -10.0, max: 0.0),
good <- float(min: 1.0, max: 1_000.0),
lambda <- float(min: 0.001, max: 1.0)
) do
assert TerrainAnalysis.diffraction_param(h, bad, good, lambda) == 0.0
assert TerrainAnalysis.diffraction_param(h, good, bad, lambda) == 0.0
assert TerrainAnalysis.diffraction_param(h, good, good, bad) == 0.0
end
end
end
describe "fresnel_radius/3" do
property "is non-negative and symmetric in d1 and d2" do
check all(
d1 <- float(min: 1.0, max: 100_000.0),
d2 <- float(min: 1.0, max: 100_000.0),
lambda <- float(min: 0.001, max: 1.0)
) do
r1 = TerrainAnalysis.fresnel_radius(d1, d2, lambda)
r2 = TerrainAnalysis.fresnel_radius(d2, d1, lambda)
assert r1 >= 0
assert_in_delta r1, r2, 1.0e-9
end
end
property "returns 0 when either distance is non-positive" do
check all(
bad <- float(min: -10.0, max: 0.0),
good <- float(min: 1.0, max: 1_000.0),
lambda <- float(min: 0.001, max: 1.0)
) do
assert TerrainAnalysis.fresnel_radius(bad, good, lambda) == 0
assert TerrainAnalysis.fresnel_radius(good, bad, lambda) == 0
end
end
end
describe "earth_bulge/3" do
property "is non-negative for fractions in [0, 1] and positive k" do
check all(
frac <- float(min: 0.0, max: 1.0),
dist_km <- float(min: 0.0, max: 500.0),
k <- float(min: 0.5, max: 10.0)
) do
assert TerrainAnalysis.earth_bulge(frac, dist_km, k) >= 0
end
end
property "is zero at both endpoints regardless of distance or k" do
check all(
dist_km <- float(min: 0.0, max: 500.0),
k <- float(min: 0.5, max: 10.0)
) do
assert TerrainAnalysis.earth_bulge(0.0, dist_km, k) == 0.0
assert TerrainAnalysis.earth_bulge(1.0, dist_km, k) == 0.0
end
end
property "is symmetric around the midpoint" do
check all(
frac <- float(min: 0.0, max: 0.5),
dist_km <- float(min: 1.0, max: 500.0),
k <- float(min: 0.5, max: 10.0)
) do
left = TerrainAnalysis.earth_bulge(frac, dist_km, k)
right = TerrainAnalysis.earth_bulge(1.0 - frac, dist_km, k)
assert_in_delta left, right, 1.0e-9
end
end
end
describe "k_factor/1" do
property "returns a positive value for any dN/dh" do
check all(dn_dh <- float(min: -500.0, max: 500.0)) do
assert TerrainAnalysis.k_factor(dn_dh) > 0
end
end
property "sub-refractive gradients (dN/dh > 0) yield k < 1" do
check all(dn_dh <- float(min: 1.0, max: 500.0)) do
assert TerrainAnalysis.k_factor(dn_dh) < 1.0
end
end
property "super-refractive gradients (dN/dh < 0, above ducting cap) yield k > 1" do
# Stay well above the -157 ducting threshold so the 100.0 cap doesn't kick in.
check all(dn_dh <- float(min: -100.0, max: -1.0)) do
assert TerrainAnalysis.k_factor(dn_dh) > 1.0
end
end
end
describe "analyse/4 invariants" do
# Flat sea-level terrain with both antennas raised high enough to
# clear the midpoint earth bulge *and* the first Fresnel zone: the
# beam must be unobstructed and produce zero diffraction loss.
# Antenna floor is chosen from the worst-case (max dist, max freq)
# so every generated case is guaranteed physically clear.
property "clear LoS (flat terrain, high antennas) produces 0 diffraction loss" do
max_dist_km = 50.0
min_freq_ghz = 5.0
# Midpoint earth bulge + first-Fresnel clearance with margin.
# bulge(0.5, 50km) ≈ 36.8m; r1(25km, 25km, 0.06m) ≈ 27.4m → ~64m.
ant_floor_m = 120.0
check all(
n_segs <- integer(4..20),
dist_km <- float(min: 5.0, max: max_dist_km),
freq_ghz <- float(min: min_freq_ghz, max: 50.0),
extra_h <- float(min: 0.0, max: 400.0)
) do
profile = flat_profile(n_segs, dist_km)
ant_h = ant_floor_m + extra_h
result =
TerrainAnalysis.analyse(profile, dist_km, freq_ghz, ant_ht_a: ant_h, ant_ht_b: ant_h)
assert result.diffraction_db == 0
assert result.obstructed_count == 0
assert result.verdict in ["CLEAR", "FRESNEL_MINOR"]
end
end
property "diffraction loss is always non-negative" do
check all(
n_segs <- integer(4..12),
dist_km <- float(min: 5.0, max: 100.0),
freq_ghz <- float(min: 1.0, max: 50.0),
peak_elev <- float(min: 0.0, max: 2_000.0)
) do
profile = ridge_profile(n_segs, dist_km, peak_elev)
result = TerrainAnalysis.analyse(profile, dist_km, freq_ghz)
assert result.diffraction_db >= 0
end
end
# Deygout degenerate case: when only one interior point rises above
# the direct ray, the three-edge method collapses to a single
# knife-edge and the total loss must equal the principal edge's loss.
property "single-obstacle path: total loss equals the principal knife-edge loss" do
check all(
dist_km <- float(min: 20.0, max: 80.0),
freq_ghz <- float(min: 2.0, max: 24.0),
peak_elev <- float(min: 200.0, max: 1_500.0)
) do
# 3-segment profile: endpoints at 0, single interior peak.
profile = [
%{lat: 32.9, lon: -97.0, d: 0.0, elev: 0.0, dist_km: 0.0},
%{lat: 33.0, lon: -97.0, d: 0.5, elev: peak_elev, dist_km: dist_km * 0.5},
%{lat: 33.1, lon: -97.0, d: 1.0, elev: 0.0, dist_km: dist_km}
]
result = TerrainAnalysis.analyse(profile, dist_km, freq_ghz)
# Hand-compute the knife-edge loss for the single interior point.
lambda_m = 0.3 / freq_ghz
d1_m = dist_km * 0.5 * 1000
d2_m = dist_km * 0.5 * 1000
# h = peak above direct ray = peak_elev + earth_bulge - 0 (endpoints)
bulge = TerrainAnalysis.earth_bulge(0.5, dist_km)
h = peak_elev + bulge
nu = TerrainAnalysis.diffraction_param(h, d1_m, d2_m, lambda_m)
expected = TerrainAnalysis.knife_edge_loss(nu)
assert_in_delta result.diffraction_db, expected, 0.001
end
end
end
# ── helpers ────────────────────────────────────────────────────
defp flat_profile(n_segs, dist_km) do
for i <- 0..n_segs do
f = i / n_segs
%{lat: 32.9 + f * 0.1, lon: -97.0, d: f, elev: 0.0, dist_km: f * dist_km}
end
end
defp ridge_profile(n_segs, dist_km, peak_elev) do
mid = div(n_segs, 2)
for i <- 0..n_segs do
f = i / n_segs
elev = if i == mid, do: peak_elev, else: 0.0
%{lat: 32.9 + f * 0.1, lon: -97.0, d: f, elev: elev, dist_km: f * dist_km}
end
end
end

119
test/microwaveprop_web/skew_t_property_test.exs Normal file
View file

@ -0,0 +1,119 @@
defmodule MicrowavepropWeb.SkewTPropertyTest do
@moduledoc """
StreamData property tests for `MicrowavepropWeb.SkewT`'s atmospheric
math helpers.
Invariants covered:
* Magnus saturation vapor pressure is positive and strictly
monotonic in temperature.
* `temperature_from_mixing_ratio/2` is the inverse of the
ws = 622·es / (p − es) forward relation.
* Along a dry adiabat, temperature decreases monotonically as
pressure falls (lift a parcel → it cools).
* At the 1000 mb reference pressure, dry-adiabat T equals θ − 273.15
(definition of potential temperature).
* `project/3` reproduces the chart corners and preserves pressure
ordering on the y-axis.
"""
use ExUnit.Case, async: true
use ExUnitProperties
alias MicrowavepropWeb.SkewT
describe "saturation_vapor_pressure/1" do
property "is strictly positive across the troposphere" do
check all(t_c <- float(min: -80.0, max: 50.0)) do
assert SkewT.saturation_vapor_pressure(t_c) > 0
end
end
property "is strictly monotonic in temperature" do
check all(
t_c <- float(min: -80.0, max: 50.0),
delta <- float(min: 0.001, max: 10.0)
) do
lo = SkewT.saturation_vapor_pressure(t_c)
hi = SkewT.saturation_vapor_pressure(t_c + delta)
assert hi > lo
end
end
end
describe "temperature_from_mixing_ratio/2" do
property "round-trips T → ws → T for realistic (T, p) pairs" do
check all(
t_c <- float(min: -30.0, max: 35.0),
p_mb <- float(min: 300.0, max: 1050.0)
) do
es = SkewT.saturation_vapor_pressure(t_c)
# Only test where ws is well-defined (p > es).
if p_mb > es + 1.0 do
ws_g_kg = 622.0 * es / (p_mb - es)
recovered = SkewT.temperature_from_mixing_ratio(ws_g_kg, p_mb)
assert_in_delta recovered, t_c, 0.01
end
end
end
end
describe "dry_adiabat_temperature/2" do
property "equals θ − 273.15 at the 1000 mb reference pressure" do
check all(theta_k <- float(min: 200.0, max: 400.0)) do
t_c = SkewT.dry_adiabat_temperature(theta_k, 1000.0)
assert_in_delta t_c, theta_k - 273.15, 1.0e-6
end
end
property "temperature decreases as pressure decreases along an adiabat" do
check all(
theta_k <- float(min: 250.0, max: 360.0),
p_hi <- float(min: 500.0, max: 1050.0),
drop <- float(min: 10.0, max: 400.0)
) do
p_lo = p_hi - drop
if p_lo > 50.0 do
t_hi = SkewT.dry_adiabat_temperature(theta_k, p_hi)
t_lo = SkewT.dry_adiabat_temperature(theta_k, p_lo)
assert t_hi > t_lo
end
end
end
end
describe "project/3" do
setup do
%{chart: SkewT.chart_config([])}
end
property "projects p_bot to the plot bottom and p_top to the plot top",
%{chart: chart} do
check all(t_c <- float(min: chart.t_min, max: chart.t_max)) do
{_x_bot, y_bot} = SkewT.project(t_c, chart.p_bot, chart)
{_x_top, y_top} = SkewT.project(t_c, chart.p_top, chart)
# SVG y grows downward: the bottom of the chart has a larger y.
assert y_bot > y_top
assert_in_delta y_bot, chart.padding_top + chart.plot_height, 1.0e-6
assert_in_delta y_top, chart.padding_top, 1.0e-6
end
end
property "pressure ordering is inverted on the y-axis (higher p → larger y)",
%{chart: chart} do
check all(
p_hi <- float(min: 200.0, max: 1050.0),
p_lo <- float(min: 100.0, max: 1000.0),
t_c <- float(min: -20.0, max: 20.0)
) do
if p_hi > p_lo + 1.0 do
{_x1, y_hi} = SkewT.project(t_c, p_hi, chart)
{_x2, y_lo} = SkewT.project(t_c, p_lo, chart)
assert y_hi > y_lo
end
end
end
end
end