defmodule Microwaveprop.Propagation.PathCompute do
  @moduledoc """
  Pure(ish) path-calculator engine extracted from `MicrowavepropWeb.PathLive`.

  Given a source / destination / band / station-params tuple, runs the
  same pipeline `/path` does: terrain analysis (ITU-R P.526), HRRR
  profile lookup at 9 evenly-spaced points, sounding & ionosphere
  readouts, native HRRR duct info, composite scoring, loss + power
  budgets, and the 48-hour propagation forecast at the path midpoint.

  Used by both:
    * the live `MicrowavepropWeb.PathLive` page (live recompute)
    * the rover-planning `RoverPathProfileWorker` (background cache —
      result map is stored on the matching `RoverPlanning.Path`).
  """

  import Microwaveprop.Geo, only: [haversine_km: 4, bearing_deg: 4]

  alias Microwaveprop.Ionosphere
  alias Microwaveprop.Propagation
  alias Microwaveprop.Propagation.BandConfig
  alias Microwaveprop.Propagation.ProfilesFile
  alias Microwaveprop.Propagation.Scorer
  alias Microwaveprop.Propagation.SporadicE
  alias Microwaveprop.Terrain.ElevationClient
  alias Microwaveprop.Terrain.TerrainAnalysis
  alias Microwaveprop.Weather

  require Logger

  @type station_params :: %{
          required(:src_height_m) => float(),
          required(:dst_height_m) => float(),
          required(:tx_power_dbm) => float(),
          required(:src_gain_dbi) => float(),
          required(:dst_gain_dbi) => float(),
          optional(:src_height_ft) => float(),
          optional(:dst_height_ft) => float(),
          optional(:tx_power_mw) => float()
        }

  @type result :: %{
          source: map(),
          destination: map(),
          station_params: station_params(),
          band_config: map(),
          band_mhz: integer(),
          freq_ghz: float(),
          dist_km: float(),
          bearing: float(),
          terrain: map() | nil,
          conditions: map() | nil,
          scoring: map() | nil,
          loss_budget: map(),
          power_budget: map(),
          forecast: list(),
          hrrr_count: integer(),
          hrrr_points: list(),
          ionosphere: map() | nil,
          sounding: map() | nil
        }

  @type progress_fn :: (pos_integer(), pos_integer(), String.t() -> any())

  @stages [
    "Resolving locations",
    "Analyzing terrain",
    "Loading HRRR profile grid",
    "Fetching atmospheric data along path",
    "Loading sounding & duct data",
    "Scoring propagation conditions",
    "Computing link & power budget",
    "Loading 48-hour forecast",
    "Loading ionosphere data"
  ]
  @total_stages length(@stages)

  @doc "Total number of progress stages emitted by `compute/5`."
  @spec total_stages() :: integer()
  def total_stages, do: @total_stages

  @doc "Ordered stage labels emitted by `compute/5` (1-indexed)."
  @spec stage_labels() :: list(String.t())
  def stage_labels, do: @stages

  @doc """
  Run the full path-calculator pipeline. Returns `{:ok, result}` on
  success, `{:error, reason}` when either endpoint fails to resolve
  (callsign / grid / coords).

  Options:
    * `:on_progress` — `(step, total, label)` callback invoked at the
      start of each stage. Defaults to a no-op. Use it to push live
      progress updates to a `Phoenix.LiveView`.
  """
  @spec compute(String.t(), String.t(), integer(), station_params(), keyword()) ::
          {:ok, result()} | {:error, term()}
  def compute(source, dest, band_mhz, station_params, opts \\ []) do
    on_progress = Keyword.get(opts, :on_progress, fn _step, _total, _label -> :ok end)
    report = fn step -> on_progress.(step, @total_stages, Enum.at(@stages, step - 1)) end

    report.(1)

    with {:ok, src} <- resolve_location(source),
         {:ok, dst} <- resolve_location(dest) do
      band_config = BandConfig.get(band_mhz) || BandConfig.get(10_000)
      freq_ghz = band_mhz / 1000

      dist_km = haversine_km(src.lat, src.lon, dst.lat, dst.lon)
      bearing = bearing_deg(src.lat, src.lon, dst.lat, dst.lon)

      report.(2)
      terrain_result = compute_terrain(src, dst, dist_km, freq_ghz, station_params)

      now = DateTime.utc_now()
      midlat = (src.lat + dst.lat) / 2
      midlon = (src.lon + dst.lon) / 2

      report.(3)
      profile_valid_time = latest_profile_valid_time(now)
      profile_grid = profile_grid_for(profile_valid_time)
      sample_points = path_sample_points(src, dst, 9)

      {grid_hits, misses} = partition_grid_hits(sample_points, profile_grid, profile_valid_time)

      report.(4)
      fallback_hits = fallback_hits(misses, now)

      hrrr_points = Enum.reverse(grid_hits, fallback_hits)
      hrrr_profiles = Enum.map(hrrr_points, & &1.profile)

      report.(5)
      sounding = build_sounding_readout(midlat, midlon, now)

      native_duct =
        case Weather.nearest_native_duct_info(midlat, midlon, now) do
          {:ok, info} -> info
          {:error, :not_found} -> %{best_duct_band_ghz: nil, bulk_richardson: nil}
        end

      report.(6)

      {conditions, scoring} =
        build_scoring(hrrr_profiles, src, dst, now, band_config, native_duct)

      report.(7)
      loss_budget = compute_loss_budget(dist_km, freq_ghz, band_config, terrain_result, conditions)
      power_budget = compute_power_budget(station_params, loss_budget)

      report.(8)
      forecast = Propagation.point_forecast(band_mhz, midlat, midlon)

      report.(9)
      ionosphere = build_ionosphere_readout(band_mhz, midlat, midlon, dist_km)

      {:ok,
       %{
         source: src,
         destination: dst,
         station_params: station_params,
         band_config: band_config,
         band_mhz: band_mhz,
         freq_ghz: freq_ghz,
         dist_km: dist_km,
         bearing: bearing,
         terrain: terrain_result,
         conditions: conditions,
         scoring: scoring,
         loss_budget: loss_budget,
         power_budget: power_budget,
         forecast: forecast,
         hrrr_count: length(hrrr_profiles),
         hrrr_points: hrrr_points,
         ionosphere: ionosphere,
         sounding: sounding
       }}
    end
  end

  defp compute_terrain(src, dst, dist_km, freq_ghz, station_params) do
    case ElevationClient.fetch_elevation_profile(src.lat, src.lon, dst.lat, dst.lon, 64, download: true) do
      {:ok, profile} ->
        analysis =
          TerrainAnalysis.analyse(profile, dist_km, freq_ghz,
            ant_ht_a: station_params.src_height_m,
            ant_ht_b: station_params.dst_height_m
          )

        %{profile: profile, analysis: analysis}

      {:error, reason} ->
        Logger.warning(
          "PathCompute terrain profile load failed: src=#{src.lat},#{src.lon} dst=#{dst.lat},#{dst.lon} reason=#{inspect(reason)}"
        )

        nil
    end
  end

  defp fallback_hits(misses, now) do
    misses
    |> Enum.reverse()
    |> Enum.map(&fallback_hrrr_point(&1, now))
    |> Enum.filter(& &1)
  end

  @doc "Public for PathLive's `path_forecast_detail` event."
  @spec resolve_location(String.t()) :: {:ok, map()} | {:error, String.t()}
  def resolve_location(input) do
    case MicrowavepropWeb.LocationResolver.resolve(input) do
      :empty -> {:error, "Location is required"}
      other -> other
    end
  end

  @doc "Public for PathLive's `path_forecast_detail` event."
  @spec build_ionosphere_readout(integer(), float(), float(), float()) :: map() | nil
  def build_ionosphere_readout(band_mhz, midlat, midlon, dist_km) when band_mhz in [50, 144, 222, 432] do
    case Ionosphere.nearest_foes(midlat, midlon) do
      {:ok, obs} ->
        es_score = SporadicE.es_score(obs.fo_es_mhz, band_mhz, dist_km)
        muf = SporadicE.single_hop_muf(obs.fo_es_mhz, dist_km)

        %{
          station_code: obs.station_code,
          valid_time: obs.valid_time,
          fo_es_mhz: obs.fo_es_mhz,
          fo_f2_mhz: obs.fo_f2_mhz,
          mufd_mhz: obs.mufd_mhz,
          es_score: es_score,
          es_muf_mhz: muf,
          es_in_range?: dist_km >= 500 and dist_km <= 2500
        }

      {:error, _reason} ->
        nil
    end
  end

  def build_ionosphere_readout(_band_mhz, _lat, _lon, _dist), do: nil

  defp profile_grid_for(nil), do: nil

  defp profile_grid_for(valid_time) do
    case ProfilesFile.read(valid_time) do
      {:ok, grid} ->
        grid

      {:error, reason} ->
        Logger.warning(
          "PathCompute ProfilesFile.read failed: valid_time=#{inspect(valid_time)} reason=#{inspect(reason)}"
        )

        nil
    end
  end

  defp lookup_profile_grid(grid, lat, lon) do
    {snapped_lat, snapped_lon} = ProfilesFile.snap(lat, lon)
    Map.get(grid, {snapped_lat, snapped_lon})
  end

  defp partition_grid_hits(sample_points, profile_grid, profile_valid_time) do
    Enum.reduce(sample_points, {[], []}, fn pt, acc ->
      classify_grid_hit(pt, profile_grid, profile_valid_time, acc)
    end)
  end

  defp classify_grid_hit(pt, nil, _vt, {hits, mss}), do: {hits, [pt | mss]}

  defp classify_grid_hit({label, lat, lon} = pt, grid, vt, {hits, mss}) do
    case lookup_profile_grid(grid, lat, lon) do
      nil -> {hits, [pt | mss]}
      cell -> {[%{label: label, profile: profile_from_cell(cell, lat, lon, vt)} | hits], mss}
    end
  end

  defp fallback_hrrr_point({label, lat, lon}, now) do
    case Weather.find_nearest_hrrr(lat, lon, now) do
      nil -> nil
      profile -> %{label: label, profile: profile}
    end
  end

  defp latest_profile_valid_time(now) do
    case ProfilesFile.list_valid_times() do
      [] -> nil
      times -> pick_latest_at_or_before(times, now)
    end
  end

  defp pick_latest_at_or_before(times, now) do
    case Enum.filter(times, fn t -> DateTime.compare(t, now) != :gt end) do
      [] -> List.first(times)
      past -> Enum.max(past, DateTime)
    end
  end

  defp profile_from_cell(cell, lat, lon, valid_time) do
    cell
    |> Map.put_new(:lat, lat)
    |> Map.put_new(:lon, lon)
    |> Map.put_new(:valid_time, valid_time)
    |> Map.put_new(:min_refractivity_gradient, Map.get(cell, :native_min_gradient))
    |> Map.put_new(:surface_refractivity, nil)
    |> Map.put_new(:ducting_detected, Map.get(cell, :duct_count, 0) > 0)
    |> Map.put_new(:duct_characteristics, nil)
  end

  defp path_sample_points(src, dst, count) when count >= 2 do
    Enum.map(1..count, fn i ->
      t = (i - 1) / (count - 1)
      lat = src.lat + (dst.lat - src.lat) * t
      lon = src.lon + (dst.lon - src.lon) * t

      label =
        cond do
          i == 1 -> "Source"
          i == count -> "Destination"
          i * 2 == count + 1 -> "Midpoint"
          true -> "#{round(t * 100)}%"
        end

      {label, lat, lon}
    end)
  end

  defp build_sounding_readout(midlat, midlon, now) do
    case Weather.nearest_sounding_to(midlat, midlon, now) do
      {:ok, sounding} ->
        station = sounding.station

        distance_km =
          if station, do: haversine_km(midlat, midlon, station.lat, station.lon)

        %{
          station_code: station && station.station_code,
          station_name: station && station.name,
          observed_at: sounding.observed_at,
          ducting_detected: sounding.ducting_detected,
          min_refractivity_gradient: sounding.min_refractivity_gradient,
          boundary_layer_depth_m: sounding.boundary_layer_depth_m,
          precipitable_water_mm: sounding.precipitable_water_mm,
          distance_km: distance_km
        }

      {:error, :not_found} ->
        nil
    end
  end

  defp build_scoring([], _src, _dst, _now, _band_config, _native_duct), do: {nil, nil}

  defp build_scoring(profiles, src, dst, now, band_config, native_duct) do
    {temps, dewpoints, pressures, gradients, bl_depths, pwats} = collect_profile_fields(profiles)

    if temps == [] or dewpoints == [] do
      {nil, nil}
    else
      avg_temp_c = Enum.sum(temps) / length(temps)
      avg_dewpoint_c = Enum.sum(dewpoints) / length(dewpoints)

      conditions =
        build_conditions(
          avg_temp_c,
          avg_dewpoint_c,
          src,
          dst,
          now,
          {pressures, gradients, bl_depths, pwats},
          native_duct
        )

      scoring = Scorer.composite_score(conditions, band_config)
      {conditions, scoring}
    end
  end

  defp collect_profile_fields(profiles) do
    Enum.reduce(profiles, {[], [], [], [], [], []}, fn p, {ts, ds, ps, gs, bs, ws} ->
      {
        if(p.surface_temp_c == nil, do: ts, else: [p.surface_temp_c | ts]),
        if(p.surface_dewpoint_c == nil, do: ds, else: [p.surface_dewpoint_c | ds]),
        if(p.surface_pressure_mb == nil, do: ps, else: [p.surface_pressure_mb | ps]),
        if(p.min_refractivity_gradient == nil,
          do: gs,
          else: [p.min_refractivity_gradient | gs]
        ),
        if(p.hpbl_m == nil, do: bs, else: [p.hpbl_m | bs]),
        if(p.pwat_mm == nil, do: ws, else: [p.pwat_mm | ws])
      }
    end)
  end

  defp build_conditions(avg_temp_c, avg_dewpoint_c, src, dst, now, {pressures, gradients, bl_depths, pwats}, native_duct) do
    %{
      abs_humidity: Scorer.absolute_humidity(avg_temp_c, avg_dewpoint_c),
      temp_f: Scorer.c_to_f(avg_temp_c),
      dewpoint_f: Scorer.c_to_f(avg_dewpoint_c),
      temp_c: avg_temp_c,
      dewpoint_c: avg_dewpoint_c,
      wind_speed_kts: nil,
      sky_cover_pct: nil,
      utc_hour: now.hour,
      utc_minute: now.minute,
      month: now.month,
      latitude: (src.lat + dst.lat) / 2,
      longitude: (src.lon + dst.lon) / 2,
      pressure_mb: if(pressures != [], do: Enum.min(pressures)),
      prev_pressure_mb: nil,
      rain_rate_mmhr: 0.0,
      min_refractivity_gradient: if(gradients != [], do: Enum.min(gradients)),
      bl_depth_m: if(bl_depths != [], do: Enum.sum(bl_depths) / length(bl_depths)),
      pwat_mm: if(pwats != [], do: Enum.sum(pwats) / length(pwats)),
      best_duct_band_ghz: native_duct[:best_duct_band_ghz],
      bulk_richardson: native_duct[:bulk_richardson]
    }
  end

  @doc false
  @spec compute_loss_budget(float(), float(), map(), map() | nil, map() | nil) :: map()
  def compute_loss_budget(dist_km, freq_ghz, band_config, terrain_result, conditions) do
    freq_mhz = freq_ghz * 1000
    fspl = 20 * :math.log10(max(dist_km, 0.001)) + 20 * :math.log10(freq_mhz) + 32.44

    o2_loss = band_config.o2_db_km * dist_km
    h2o_coeff = band_config.h2o_coeff

    abs_humidity =
      if conditions do
        Map.get(conditions, :abs_humidity) || 7.5
      else
        7.5
      end

    h2o_loss = h2o_coeff * abs_humidity * dist_km

    rain_loss =
      if conditions && conditions.rain_rate_mmhr > 0 do
        gamma = band_config.rain_k * :math.pow(conditions.rain_rate_mmhr, band_config.rain_alpha)
        gamma * dist_km
      else
        0.0
      end

    diffraction_loss =
      if terrain_result do
        terrain_result.analysis.diffraction_db * 1.0
      else
        0.0
      end

    total = fspl + o2_loss + h2o_loss + rain_loss + diffraction_loss

    %{
      fspl: Float.round(fspl, 1),
      o2: Float.round(o2_loss, 2),
      h2o: Float.round(h2o_loss, 2),
      rain: Float.round(rain_loss, 2),
      diffraction: Float.round(diffraction_loss, 1),
      total: Float.round(total, 1)
    }
  end

  @doc false
  @spec compute_power_budget(station_params(), map()) :: map()
  def compute_power_budget(station_params, loss_budget) do
    tx_power_dbm = station_params.tx_power_dbm
    eirp_dbm = tx_power_dbm + station_params.src_gain_dbi
    rx_power_dbm = eirp_dbm - loss_budget.total + station_params.dst_gain_dbi

    rx_sensitivity_cw = -140.0
    rx_sensitivity_ssb = -130.0
    margin_cw = rx_power_dbm - rx_sensitivity_cw
    margin_ssb = rx_power_dbm - rx_sensitivity_ssb

    %{
      tx_power_dbm: Float.round(tx_power_dbm, 1),
      eirp_dbm: Float.round(eirp_dbm, 1),
      rx_power_dbm: Float.round(rx_power_dbm, 1),
      margin_cw: Float.round(margin_cw, 1),
      margin_ssb: Float.round(margin_ssb, 1)
    }
  end
end