prop/lib/microwaveprop/terrain/terrain_analysis.ex

defmodule Microwaveprop.Terrain.TerrainAnalysis do
  @moduledoc false

  @earth_radius_m 6_371_000.0

  @doc """
  Analyse an elevation profile for LOS clearance, Fresnel zone penetration,
  and knife-edge diffraction loss.

  K=4/3 standard atmosphere. Antenna heights default to 0m (conservative).
  """
  def analyse(profile, dist_km, freq_ghz, ant_ht_a_m \\ 0.0, ant_ht_b_m \\ 0.0) do
    lambda_m = 0.3 / freq_ghz
    n = length(profile) - 1

    first = hd(profile)
    last = List.last(profile)
    elev1 = first.elev + ant_ht_a_m
    elev2 = last.elev + ant_ht_b_m

    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)
        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,
          beam: beam,
          eff_terrain: eff_terrain,
          bulge: bulge,
          r1: r1,
          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)

    {diffraction_db, verdict} = compute_verdict(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,
      max_elevation_m: max_elevation_m,
      min_clearance_m: min_clearance_m,
      obstructed_count: length(obstructed),
      fresnel_hit_count: length(fresnel_hit)
    }
  end

  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

  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

  def knife_edge_loss(v) when v <= -0.7787, do: 0

  def knife_edge_loss(v) when v <= 0 do
    max(0, -20 * :math.log10(0.5 - 0.62 * v))
  end

  def knife_edge_loss(v) when v <= 1 do
    max(0, -20 * :math.log10(0.5 * :math.exp(-0.95 * v)))
  end

  def knife_edge_loss(v) when v <= 2.4 do
    inner = max(0, 0.1184 - (0.38 - 0.1 * v) ** 2)
    max(0, -20 * :math.log10(0.4 - :math.sqrt(inner)))
  end

  def knife_edge_loss(v) do
    # Asymptotic: loss = 20·log10(v) + 13 dB
    max(0, 20 * :math.log10(v) + 13.0)
  end

  defp beam_height(frac, elev1, elev2) do
    elev1 + frac * (elev2 - elev1)
  end

  defp compute_verdict(obstructed, _fresnel_hit) when obstructed != [] do
    worst = Enum.min_by(obstructed, & &1.clearance)

    v_blocked =
      if worst.r1 > 0 do
        abs(worst.clearance) / worst.r1 + 0.5
      else
        2.0
      end

    {max(0, knife_edge_loss(v_blocked)), "BLOCKED"}
  end

  defp compute_verdict(_obstructed, fresnel_hit) when fresnel_hit != [] do
    worst = Enum.min_by(fresnel_hit, & &1.f1_clear)

    v_fresnel =
      if worst.r1 > 0 do
        -worst.f1_clear / worst.r1
      else
        0
      end

    db = max(0, knife_edge_loss(v_fresnel))
    verdict = if db > 3, do: "FRESNEL_PARTIAL", else: "FRESNEL_MINOR"
    {db, verdict}
  end

  defp compute_verdict(_obstructed, _fresnel_hit) do
    {0, "CLEAR"}
  end
end