diff --git a/test/microwaveprop/geo_property_test.exs b/test/microwaveprop/geo_property_test.exs
new file mode 100644
index 00000000..e06e995d
--- /dev/null
+++ b/test/microwaveprop/geo_property_test.exs
@@ -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
diff --git a/test/microwaveprop/terrain/terrain_analysis_property_test.exs b/test/microwaveprop/terrain/terrain_analysis_property_test.exs
new file mode 100644
index 00000000..9d5f99b5
--- /dev/null
+++ b/test/microwaveprop/terrain/terrain_analysis_property_test.exs
@@ -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
diff --git a/test/microwaveprop_web/skew_t_property_test.exs b/test/microwaveprop_web/skew_t_property_test.exs
new file mode 100644
index 00000000..5b5e91c9
--- /dev/null
+++ b/test/microwaveprop_web/skew_t_property_test.exs
@@ -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