Graham McIntire
51e390959c
277 @spec/@type annotations added to 58 files covering all public APIs: contexts (propagation, radio, weather, terrain, beacons, commercial), GRIB2 decoders, terrain analysis, duct detection, rain scatter, CSV/ADIF import, weather clients, and all Ecto schemas. Dialyzer passes with 0 errors.
227 lines
7.8 KiB
Elixir
227 lines
7.8 KiB
Elixir
defmodule Microwaveprop.Terrain.Viewshed do
|
||
@moduledoc false
|
||
|
||
alias Microwaveprop.Terrain.Srtm
|
||
alias Microwaveprop.Terrain.TerrainAnalysis
|
||
|
||
@type boundary_point :: %{bearing: non_neg_integer(), reach_km: float(), lat: float(), lon: float()}
|
||
|
||
@type viewshed_result :: %{
|
||
origin: %{lat: float(), lon: float()},
|
||
boundary: [boundary_point()]
|
||
}
|
||
|
||
@earth_radius_km 6371.0
|
||
@default_angular_step 2
|
||
@default_max_range_km 50
|
||
|
||
@doc "Compute destination lat/lon given origin, bearing (degrees), and distance (km)."
|
||
@spec destination_point(float(), float(), float(), float()) :: {float(), float()}
|
||
def destination_point(lat, lon, bearing_deg, dist_km) do
|
||
lat_rad = deg_to_rad(lat)
|
||
lon_rad = deg_to_rad(lon)
|
||
brg_rad = deg_to_rad(bearing_deg)
|
||
d = dist_km / @earth_radius_km
|
||
|
||
lat2 =
|
||
:math.asin(
|
||
:math.sin(lat_rad) * :math.cos(d) +
|
||
:math.cos(lat_rad) * :math.sin(d) * :math.cos(brg_rad)
|
||
)
|
||
|
||
lon2 =
|
||
lon_rad +
|
||
:math.atan2(
|
||
:math.sin(brg_rad) * :math.sin(d) * :math.cos(lat_rad),
|
||
:math.cos(d) - :math.sin(lat_rad) * :math.sin(lat2)
|
||
)
|
||
|
||
{rad_to_deg(lat2), rad_to_deg(lon2)}
|
||
end
|
||
|
||
@doc """
|
||
Find the max clear distance along a ray from TerrainAnalysis points.
|
||
Skips endpoints (first/last), returns the dist_km of the last clear
|
||
interior point before the first obstruction.
|
||
"""
|
||
@spec find_reach_km([TerrainAnalysis.analysis_point()], float()) :: float()
|
||
def find_reach_km(points, max_range_km) do
|
||
interior = Enum.slice(points, 1..-2//1)
|
||
|
||
case find_first_obstructed_index(interior) do
|
||
nil ->
|
||
max_range_km
|
||
|
||
0 ->
|
||
0.0
|
||
|
||
idx ->
|
||
Enum.at(interior, idx - 1).dist_km
|
||
end
|
||
end
|
||
|
||
@doc """
|
||
Compute effective reach accounting for both terrain and propagation conditions.
|
||
Higher propagation scores indicate ducting potential, which lets signals
|
||
propagate beyond terrain obstructions via atmospheric waveguide.
|
||
"""
|
||
@spec effective_reach_km(TerrainAnalysis.analysis_result(), float(), number()) :: float()
|
||
def effective_reach_km(analysis, max_range_km, score) do
|
||
case analysis.verdict do
|
||
"CLEAR" ->
|
||
max_range_km
|
||
|
||
"FRESNEL_MINOR" ->
|
||
max_range_km * 0.9
|
||
|
||
"FRESNEL_PARTIAL" ->
|
||
max_range_km * 0.7
|
||
|
||
"BLOCKED" ->
|
||
terrain_factor = terrain_reach_factor(analysis.diffraction_db)
|
||
ducting_factor = ducting_reach_factor(score)
|
||
max_range_km * max(terrain_factor, ducting_factor)
|
||
end
|
||
end
|
||
|
||
# Convert diffraction loss to a range reduction factor (0.0–1.0).
|
||
# Mild diffraction still allows significant propagation;
|
||
# heavy blockage attenuates sharply.
|
||
defp terrain_reach_factor(db) when db <= 3, do: 0.8
|
||
defp terrain_reach_factor(db) when db <= 6, do: 0.5
|
||
defp terrain_reach_factor(db) when db <= 12, do: 0.3
|
||
defp terrain_reach_factor(db) when db <= 20, do: 0.15
|
||
defp terrain_reach_factor(_db), do: 0.05
|
||
|
||
# Higher propagation scores mean stronger ducting potential —
|
||
# signals can ride atmospheric layers over terrain obstacles.
|
||
# Returns a floor factor so blocked paths still get range
|
||
# proportional to atmospheric conditions.
|
||
defp ducting_reach_factor(score) when score >= 80, do: 0.7
|
||
defp ducting_reach_factor(score) when score >= 65, do: 0.5
|
||
defp ducting_reach_factor(score) when score >= 50, do: 0.3
|
||
defp ducting_reach_factor(score) when score >= 33, do: 0.15
|
||
defp ducting_reach_factor(_score), do: 0.05
|
||
|
||
@doc """
|
||
Compute a terrain viewshed from a point. Returns a map with :origin
|
||
and :boundary (list of %{bearing, reach_km, lat, lon}).
|
||
|
||
Options:
|
||
- :freq_ghz — frequency for Fresnel zone calc (default 10.0)
|
||
- :max_range_km — max ray distance (default 50)
|
||
- :ant_height_m — antenna height at both ends (default 2.4)
|
||
- :angular_step — degrees between rays (default 2)
|
||
- :tiles_dir — SRTM tiles directory (default from config)
|
||
"""
|
||
@spec compute(float(), float(), keyword()) :: viewshed_result()
|
||
def compute(lat, lon, opts \\ []) do
|
||
freq_ghz = Keyword.get(opts, :freq_ghz, 10.0)
|
||
max_range_km = Keyword.get(opts, :max_range_km, @default_max_range_km)
|
||
ant_height_m = Keyword.get(opts, :ant_height_m, 2.4)
|
||
score = Keyword.get(opts, :score, 50)
|
||
angular_step = Keyword.get(opts, :angular_step, @default_angular_step)
|
||
tiles_dir = Keyword.get(opts, :tiles_dir, srtm_tiles_dir())
|
||
|
||
bearings = Enum.to_list(0..359//angular_step)
|
||
|
||
boundary =
|
||
bearings
|
||
|> Task.async_stream(
|
||
fn bearing ->
|
||
compute_ray(lat, lon, bearing, freq_ghz, max_range_km, ant_height_m, score, tiles_dir)
|
||
end,
|
||
max_concurrency: System.schedulers_online(),
|
||
timeout: 30_000,
|
||
on_timeout: :kill_task
|
||
)
|
||
|> Enum.map(fn
|
||
{:ok, result} -> result
|
||
{:exit, _} -> nil
|
||
end)
|
||
|> Enum.reject(&is_nil/1)
|
||
|> smooth_boundary(lat, lon)
|
||
|
||
%{origin: %{lat: lat, lon: lon}, boundary: boundary}
|
||
end
|
||
|
||
@doc """
|
||
Analyse a single ray's profile. Public for testing.
|
||
Returns %{reach_km: float, verdict: string}.
|
||
"""
|
||
@spec analyse_ray([TerrainAnalysis.elevation_point()], float(), float(), float(), float()) ::
|
||
%{reach_km: float(), verdict: String.t()}
|
||
def analyse_ray(profile, dist_km, freq_ghz, ant_ht_a_m, ant_ht_b_m) do
|
||
analysis = TerrainAnalysis.analyse(profile, dist_km, freq_ghz, ant_ht_a: ant_ht_a_m, ant_ht_b: ant_ht_b_m)
|
||
reach_km = find_reach_km(analysis.points, dist_km)
|
||
%{reach_km: reach_km, verdict: analysis.verdict}
|
||
end
|
||
|
||
defp compute_ray(origin_lat, origin_lon, bearing, freq_ghz, max_range_km, ant_height_m, score, tiles_dir) do
|
||
# Check terrain within the radio horizon where terrain features matter.
|
||
# Beyond this, propagation is atmospheric (ducting/scatter) and terrain
|
||
# doesn't block — earth curvature is handled by the atmosphere.
|
||
terrain_check_km = min(radio_horizon_km(ant_height_m) * 2.0, max_range_km)
|
||
{end_lat, end_lon} = destination_point(origin_lat, origin_lon, bearing, terrain_check_km)
|
||
|
||
{:ok, profile} =
|
||
Srtm.fetch_elevation_profile(origin_lat, origin_lon, end_lat, end_lon, tiles_dir, 64, download: true)
|
||
|
||
analysis =
|
||
TerrainAnalysis.analyse(profile, terrain_check_km, freq_ghz, ant_ht_a: ant_height_m, ant_ht_b: ant_height_m)
|
||
|
||
reach_km = effective_reach_km(analysis, max_range_km, score)
|
||
{reach_lat, reach_lon} = destination_point(origin_lat, origin_lon, bearing, reach_km)
|
||
|
||
%{
|
||
bearing: bearing,
|
||
reach_km: reach_km,
|
||
lat: reach_lat,
|
||
lon: reach_lon
|
||
}
|
||
end
|
||
|
||
# Radio horizon distance in km for a given antenna height (K=4/3 atmosphere)
|
||
defp radio_horizon_km(height_m) do
|
||
:math.sqrt(2 * (4 / 3) * @earth_radius_km * 1000 * height_m) / 1000
|
||
end
|
||
|
||
defp srtm_tiles_dir do
|
||
Application.get_env(:microwaveprop, :srtm_tiles_dir, Path.expand("~/srtm/tiles"))
|
||
end
|
||
|
||
# Smooth reach_km values with a circular moving average to remove spikes
|
||
# from SRTM elevation artifacts, then recompute boundary lat/lon.
|
||
defp smooth_boundary(points, _origin_lat, _origin_lon) when length(points) < 5, do: points
|
||
|
||
defp smooth_boundary(points, origin_lat, origin_lon) do
|
||
reaches = Enum.map(points, & &1.reach_km)
|
||
n = length(reaches)
|
||
arr = :array.from_list(reaches)
|
||
half = 2
|
||
|
||
smoothed_reaches =
|
||
Enum.map(0..(n - 1), fn i ->
|
||
window =
|
||
for offset <- -half..half do
|
||
:array.get(rem(i + offset + n, n), arr)
|
||
end
|
||
|
||
Enum.sum(window) / length(window)
|
||
end)
|
||
|
||
points
|
||
|> Enum.zip(smoothed_reaches)
|
||
|> Enum.map(fn {pt, reach} ->
|
||
{lat, lon} = destination_point(origin_lat, origin_lon, pt.bearing, reach)
|
||
%{pt | reach_km: reach, lat: lat, lon: lon}
|
||
end)
|
||
end
|
||
|
||
defp find_first_obstructed_index(interior) do
|
||
Enum.find_index(interior, & &1.obstructed)
|
||
end
|
||
|
||
defp deg_to_rad(deg), do: deg * :math.pi() / 180
|
||
defp rad_to_deg(rad), do: rad * 180 / :math.pi()
|
||
end
|