Graham McIntire
9fea93f59c
Terrain elevations now include the earth bulge correction, so the profile visually humps up in the middle relative to the straight LOS beam — showing how the curved earth rises into the signal path. Uses actual k-factor from HRRR refractivity when available, otherwise standard 4/3. For a 300 km path at k=4/3, the midpoint bulge is ~2,650 m (8,700 ft).
237 lines
7.1 KiB
Elixir
237 lines
7.1 KiB
Elixir
defmodule Microwaveprop.Terrain.TerrainAnalysis do
|
||
@moduledoc """
|
||
Terrain diffraction analysis per ITU-R P.526-16.
|
||
|
||
Implements:
|
||
- Knife-edge diffraction loss (Eq. 31)
|
||
- Deygout 3-edge method for multiple obstacles (Section 6)
|
||
- Dynamic k-factor from refractivity gradient (Section 2)
|
||
"""
|
||
|
||
@earth_radius_m 6_371_000.0
|
||
@earth_radius_km 6_371.0
|
||
|
||
@doc """
|
||
Analyse an elevation profile for LOS clearance, Fresnel zone penetration,
|
||
and diffraction loss per ITU-R P.526-16.
|
||
|
||
## Options
|
||
* `:ant_ht_a` - antenna height A in meters (default 0.0)
|
||
* `:ant_ht_b` - antenna height B in meters (default 0.0)
|
||
* `:k_factor` - effective earth radius factor (default 4/3, compute with `k_factor/1`)
|
||
"""
|
||
def analyse(profile, dist_km, freq_ghz, opts \\ []) do
|
||
ant_ht_a = Keyword.get(opts, :ant_ht_a, 0.0)
|
||
ant_ht_b = Keyword.get(opts, :ant_ht_b, 0.0)
|
||
k = Keyword.get(opts, :k_factor, 4 / 3)
|
||
|
||
lambda_m = 0.3 / freq_ghz
|
||
n = length(profile) - 1
|
||
|
||
first = hd(profile)
|
||
last = List.last(profile)
|
||
elev1 = first.elev + ant_ht_a
|
||
elev2 = last.elev + ant_ht_b
|
||
|
||
# Build effective profile with earth curvature correction
|
||
points =
|
||
profile
|
||
|> Enum.with_index()
|
||
|> Enum.map(fn {p, i} ->
|
||
frac = i / n
|
||
d1_m = frac * dist_km * 1000
|
||
d2_m = (1 - frac) * dist_km * 1000
|
||
|
||
beam = beam_height(frac, elev1, elev2)
|
||
bulge = earth_bulge(frac, dist_km, k)
|
||
eff_terrain = p.elev + bulge
|
||
r1 = fresnel_radius(d1_m, d2_m, lambda_m)
|
||
clearance = beam - eff_terrain
|
||
f1_clear = if r1 > 0, do: clearance - r1, else: clearance
|
||
|
||
is_endpoint = i == 0 or i == n
|
||
|
||
%{
|
||
lat: p.lat,
|
||
lon: p.lon,
|
||
d: p.d,
|
||
elev: p.elev,
|
||
dist_km: p.dist_km,
|
||
idx: i,
|
||
beam: beam,
|
||
eff_terrain: eff_terrain,
|
||
bulge: bulge,
|
||
r1: r1,
|
||
d1_m: d1_m,
|
||
d2_m: d2_m,
|
||
clearance: clearance,
|
||
f1_clear: f1_clear,
|
||
obstructed: not is_endpoint and clearance < 0,
|
||
fresnel_penetrated: not is_endpoint and r1 > 0 and f1_clear < 0 and clearance >= 0
|
||
}
|
||
end)
|
||
|
||
interior = Enum.slice(points, 1..-2//1)
|
||
obstructed = Enum.filter(interior, & &1.obstructed)
|
||
fresnel_hit = Enum.filter(interior, & &1.fresnel_penetrated)
|
||
|
||
# Deygout diffraction: compute loss across all interior points
|
||
diffraction_db = deygout_diffraction(interior, lambda_m)
|
||
|
||
verdict = classify_verdict(diffraction_db, obstructed, fresnel_hit)
|
||
|
||
elev_vals = Enum.map(profile, & &1.elev)
|
||
clearance_vals = Enum.map(interior, & &1.clearance)
|
||
max_elevation_m = Enum.max(elev_vals, fn -> 0.0 end)
|
||
min_clearance_m = if clearance_vals == [], do: 999.0, else: Enum.min(clearance_vals)
|
||
|
||
%{
|
||
points: points,
|
||
diffraction_db: diffraction_db,
|
||
verdict: verdict,
|
||
k_factor: k,
|
||
max_elevation_m: max_elevation_m,
|
||
min_clearance_m: min_clearance_m,
|
||
obstructed_count: length(obstructed),
|
||
fresnel_hit_count: length(fresnel_hit)
|
||
}
|
||
end
|
||
|
||
@doc "First Fresnel zone radius at a point between two antennas."
|
||
def fresnel_radius(d1_m, d2_m, lambda_m) do
|
||
if d1_m <= 0 or d2_m <= 0 do
|
||
0
|
||
else
|
||
:math.sqrt(lambda_m * d1_m * d2_m / (d1_m + d2_m))
|
||
end
|
||
end
|
||
|
||
@doc "Earth bulge correction in meters at fractional distance along path."
|
||
def earth_bulge(frac, dist_km, k \\ 4 / 3) do
|
||
d1 = frac * dist_km * 1000
|
||
d2 = (1 - frac) * dist_km * 1000
|
||
d1 * d2 / (2 * k * @earth_radius_m)
|
||
end
|
||
|
||
@doc """
|
||
ITU-R P.526-16 Eq. 31 — knife-edge diffraction loss in dB.
|
||
|
||
J(ν) = 6.9 + 20·log10(√((ν−0.1)² + 1) + ν − 0.1) for ν > −0.78
|
||
J(ν) = 0 for ν ≤ −0.78
|
||
"""
|
||
def knife_edge_loss(v) when v <= -0.78, do: 0
|
||
|
||
def knife_edge_loss(v) do
|
||
inner = :math.sqrt((v - 0.1) * (v - 0.1) + 1) + v - 0.1
|
||
max(0, 6.9 + 20 * :math.log10(inner))
|
||
end
|
||
|
||
@doc """
|
||
ITU-R P.526-16 diffraction parameter ν.
|
||
|
||
ν = h · √(2·(d1+d2) / (λ·d1·d2))
|
||
|
||
where h is height above direct ray (positive = above/blocked, negative = below/clear).
|
||
"""
|
||
def diffraction_param(h, d1_m, d2_m, lambda_m) do
|
||
if d1_m <= 0 or d2_m <= 0 or lambda_m <= 0 do
|
||
0.0
|
||
else
|
||
h * :math.sqrt(2 * (d1_m + d2_m) / (lambda_m * d1_m * d2_m))
|
||
end
|
||
end
|
||
|
||
@doc """
|
||
Compute effective earth radius factor k from refractivity gradient dN/dh (N-units/km).
|
||
|
||
k = 1 / (1 + a·(dN/dh)·10⁻⁶)
|
||
|
||
Standard atmosphere: dN/dh ≈ −39 → k ≈ 4/3.
|
||
Returns 4/3 for nil input.
|
||
"""
|
||
def k_factor(nil), do: 4 / 3
|
||
|
||
def k_factor(dn_dh) do
|
||
denominator = 1 + @earth_radius_km * dn_dh * 1.0e-6
|
||
|
||
if denominator <= 0.01 do
|
||
# Near or beyond ducting — cap at very large k
|
||
100.0
|
||
else
|
||
1.0 / denominator
|
||
end
|
||
end
|
||
|
||
# ── Deygout 3-edge diffraction (P.526-16 Section 6) ─────────────
|
||
|
||
defp deygout_diffraction(interior, lambda_m) do
|
||
if interior == [] do
|
||
0
|
||
else
|
||
# Find principal edge: point with highest ν on full T→R path
|
||
principal = find_principal_edge(interior, lambda_m)
|
||
|
||
if principal == nil or principal.nu <= -0.78 do
|
||
0
|
||
else
|
||
loss_main = knife_edge_loss(principal.nu)
|
||
|
||
# Find subsidiary edges on each sub-path
|
||
{before_points, after_points} = split_at_edge(interior, principal.idx)
|
||
|
||
loss_t = subsidiary_loss(before_points, lambda_m)
|
||
loss_r = subsidiary_loss(after_points, lambda_m)
|
||
|
||
max(0, loss_main + loss_t + loss_r)
|
||
end
|
||
end
|
||
end
|
||
|
||
defp find_principal_edge(points, lambda_m) do
|
||
points
|
||
|> Enum.map(fn p ->
|
||
# h = height above direct ray = -clearance (clearance negative when above)
|
||
h = -p.clearance
|
||
nu = diffraction_param(h, p.d1_m, p.d2_m, lambda_m)
|
||
%{idx: p.idx, nu: nu, d1_m: p.d1_m, d2_m: p.d2_m}
|
||
end)
|
||
|> Enum.max_by(& &1.nu, fn -> nil end)
|
||
end
|
||
|
||
defp split_at_edge(points, edge_idx) do
|
||
before = Enum.filter(points, fn p -> p.idx < edge_idx end)
|
||
after_pts = Enum.filter(points, fn p -> p.idx > edge_idx end)
|
||
{before, after_pts}
|
||
end
|
||
|
||
defp subsidiary_loss([], _lambda_m), do: 0
|
||
|
||
defp subsidiary_loss(points, lambda_m) do
|
||
# Recompute ν relative to the sub-path endpoints
|
||
# For sub-path T→principal or principal→R, the beam line is different
|
||
# Use the existing d1_m/d2_m which are relative to the full path endpoints
|
||
# This is a simplified Deygout where we use the same reference line
|
||
sub_edge = find_principal_edge(points, lambda_m)
|
||
|
||
if sub_edge == nil or sub_edge.nu <= -0.78 do
|
||
0
|
||
else
|
||
knife_edge_loss(sub_edge.nu)
|
||
end
|
||
end
|
||
|
||
# ── Verdict classification ──────────────────────────────────────
|
||
|
||
defp classify_verdict(diffraction_db, obstructed, fresnel_hit) do
|
||
cond do
|
||
obstructed != [] -> "BLOCKED"
|
||
fresnel_hit != [] and diffraction_db > 3 -> "FRESNEL_PARTIAL"
|
||
fresnel_hit != [] -> "FRESNEL_MINOR"
|
||
true -> "CLEAR"
|
||
end
|
||
end
|
||
|
||
defp beam_height(frac, elev1, elev2) do
|
||
elev1 + frac * (elev2 - elev1)
|
||
end
|
||
end
|