sidereon

defmodule Sidereon.GNSS.BroadcastComparison do
  @moduledoc """
  Broadcast-ephemeris accuracy: compare a broadcast navigation product against a
  precise SP3 product over a window.

  This is the standard broadcast-orbit / clock accuracy check (the orbit and clock
  pieces of the signal-in-space range error, SISRE). For each satellite at each
  epoch it differences the broadcast-evaluated ECEF position and clock against the
  precise SP3 values, decomposes the position error into radial / along-track /
  cross-track (RAC) components, and summarizes the differences as RMS and maximum
  statistics per satellite and overall.

  Both products are evaluated through `Sidereon.GNSS.Ephemeris`, so the frame
  (ITRF/IGS ECEF, meters), the time scale (GPST), and the clock sign convention
  (positive = satellite clock ahead of system time) are exactly as documented
  there. Only epochs where **both** sources return a valid state contribute to the
  statistics; an epoch missing from either product is skipped, never extrapolated.

  ## RAC frame

  The radial/along-track/cross-track unit vectors are built from the **precise**
  state: radial along the position vector, cross-track along the orbital angular
  momentum `r x v`, and along-track completing the right-handed triad. The SP3
  position/clock interpolation does not expose a velocity, so the velocity is
  derived by a centered finite difference of the precise position
  (`(r(t+dt) - r(t-dt)) / 2dt`, falling back to a one-sided difference at a window
  edge). The position-difference vector `broadcast - precise` is projected onto
  this triad.

  ## Expected magnitudes

  GPS LNAV broadcast orbits differ from IGS precise orbits at roughly the 1-2 m
  RMS level (3D), dominated by the along-track and radial components; Galileo and
  BeiDou MEO are comparable. A result far outside this band (tens of meters or
  more) indicates a parse or evaluation defect rather than normal broadcast error.

  The broadcast models follow IS-GPS-200 (GPS LNAV), the Galileo OS-SIS-ICD, and
  the BeiDou BDS-SIS-ICD; the precise product is SP3-c / SP3-d (IGS).

  ## Example

      {:ok, broadcast} = Sidereon.GNSS.Broadcast.load("BRDC.rnx")
      {:ok, sp3} = Sidereon.GNSS.SP3.load("igs.sp3")

      report =
        Sidereon.GNSS.BroadcastComparison.compare(broadcast, sp3, ["G01", "E11"], %{
          from: ~N[2020-06-25 02:00:00],
          to: ~N[2020-06-25 04:00:00],
          step_s: 300
        })

      report.overall.orbit_3d_rms_m   # 3D orbit RMS over all satellites, meters
      report.per_satellite["G01"].radial_rms_m

  ## Broadcast and precise products on a sample day

  Comparing the GPS LNAV broadcast message against the GBM precise SP3 product for
  2020 day-of-year 177, all GPS satellites over a multi-hour window at a 15 min
  step, gives an overall 3D orbit RMS of about **1.5 m** (max ~4 m), split as
  roughly 1.1 m radial / 0.9 m along-track / 0.5 m cross-track, the expected GPS
  broadcast accuracy. The raw clock differences (`clock_rms_m`) are larger
  (several meters) because the broadcast and precise clocks are referenced to
  different time datums, which differ by a common per-epoch offset that drifts
  over the day. `clock_datum_removed_rms_m` removes that common offset (the
  per-epoch median over satellites) and reports the actual signal-in-space clock
  error, which is several times smaller.
  The same call works against a broadcast product by calling
  `compare_broadcast_vs_precise/4` directly with loaded products.

  The orbit/clock differencing, RAC decomposition, finite-difference velocity, and
  the RMS/median/datum statistics are computed in the Rust core; this module
  validates inputs, marshals the per-epoch evaluation keys across the NIF, and
  rebuilds the public report.
  """

  alias Sidereon.GNSS.Broadcast
  alias Sidereon.GNSS.Ephemeris
  alias Sidereon.GNSS.SP3
  alias Sidereon.GNSS.Time
  alias Sidereon.NIF

  defmodule Stats do
    @moduledoc """
    Orbit and clock difference statistics for one satellite (or the overall set).

    All values are meters except `count` (the number of compared epochs).
    `orbit_3d_rms_m` / `orbit_3d_max_m` are the Euclidean position-difference
    magnitudes. `radial_*`, `along_*`, `cross_*` are the RMS and max of the signed
    RAC components of the position difference (`broadcast - precise`).
    `clock_rms_m` / `clock_max_m` are the **raw** satellite-clock differences
    scaled to meters by the speed of light; they are `nil` when neither product
    carried a clock estimate for any compared epoch. They are dominated by the
    per-epoch common reference-clock offset between the two products' time
    datums. `clock_datum_removed_rms_m` / `clock_datum_removed_max_m` are the
    same differences after that per-epoch common offset (the median over all
    satellites at the epoch) is removed, the actual signal-in-space clock error
    (the SISRE clock term), typically several times smaller than the raw value.
    """
    @enforce_keys [
      :count,
      :orbit_3d_rms_m,
      :orbit_3d_max_m,
      :radial_rms_m,
      :radial_max_m,
      :along_rms_m,
      :along_max_m,
      :cross_rms_m,
      :cross_max_m,
      :clock_rms_m,
      :clock_max_m,
      :clock_datum_removed_rms_m,
      :clock_datum_removed_max_m
    ]
    defstruct [
      :count,
      :orbit_3d_rms_m,
      :orbit_3d_max_m,
      :radial_rms_m,
      :radial_max_m,
      :along_rms_m,
      :along_max_m,
      :cross_rms_m,
      :cross_max_m,
      :clock_rms_m,
      :clock_max_m,
      :clock_datum_removed_rms_m,
      :clock_datum_removed_max_m
    ]

    @type t :: %__MODULE__{
            count: non_neg_integer(),
            orbit_3d_rms_m: float() | nil,
            orbit_3d_max_m: float() | nil,
            radial_rms_m: float() | nil,
            radial_max_m: float() | nil,
            along_rms_m: float() | nil,
            along_max_m: float() | nil,
            cross_rms_m: float() | nil,
            cross_max_m: float() | nil,
            clock_rms_m: float() | nil,
            clock_max_m: float() | nil,
            clock_datum_removed_rms_m: float() | nil,
            clock_datum_removed_max_m: float() | nil
          }
  end

  defmodule Report do
    @moduledoc """
    The result of a broadcast and precise product comparison.

    `per_satellite` maps each satellite id to its `Stats`; `overall` aggregates
    every compared epoch across all satellites. `missing` lists `{satellite_id,
    count}` pairs counting epochs that were skipped because one or both products
    had no valid state there.
    """
    @enforce_keys [:per_satellite, :overall, :missing]
    defstruct [:per_satellite, :overall, :missing]

    @type t :: %__MODULE__{
            per_satellite: %{String.t() => Stats.t()},
            overall: Stats.t(),
            missing: [{String.t(), non_neg_integer()}]
          }
  end

  @doc """
  Compare a broadcast product against a precise SP3 product over `window`.

  `broadcast` is an `Sidereon.GNSS.Broadcast` handle and `precise` an
  `Sidereon.GNSS.SP3` handle for the same day; `sat_ids` is a list of canonical
  RINEX tokens; `window` is the `Sidereon.GNSS.Ephemeris` window map (`:from`,
  `:to`, `:step_s`). Returns an `Sidereon.GNSS.BroadcastComparison.Report`.

  The window should sit within both the SP3 file span and the broadcast records'
  fit intervals; epochs missing from either product are counted in
  `report.missing` and excluded from the statistics.
  """
  @spec compare(Broadcast.t(), SP3.t(), [String.t()], Ephemeris.window()) :: Report.t()
  def compare(%Broadcast{handle: broadcast}, %SP3{handle: precise}, sat_ids, %{from: from, to: to, step_s: step_s})
      when is_list(sat_ids) do
    # The interface picks the two start anchors; the core driver builds the
    # regular per-epoch grid and the velocity neighbours. The broadcast axis is a
    # continuous J2000 second (the value `Sidereon.GNSS.Broadcast.position`
    # marshals); the precise axis is a split Julian date advanced in lockstep.
    half = round(step_s / 2.0)
    {:ok, broadcast_t0_j2000_s} = Time.epoch_to_j2000_seconds_fractional(from)
    {:ok, broadcast_t1_j2000_s} = Time.epoch_to_j2000_seconds_fractional(to)
    {precise_start_jd_whole, precise_start_fraction} = Time.epoch_to_split_jd(from)

    {overall, per_satellite, missing} =
      NIF.broadcast_comparison(
        broadcast,
        precise,
        sat_ids,
        broadcast_t0_j2000_s,
        broadcast_t1_j2000_s,
        precise_start_jd_whole,
        precise_start_fraction,
        step_s / 1.0,
        half / 1.0
      )

    %Report{
      overall: decode_stats(overall),
      per_satellite: Map.new(per_satellite, fn {sat, stats} -> {sat, decode_stats(stats)} end),
      missing: Enum.map(missing, fn {sat, count} -> {sat, count} end)
    }
  end

  defp decode_stats(
         {count,
          [
            orbit_3d_rms_m,
            orbit_3d_max_m,
            radial_rms_m,
            radial_max_m,
            along_rms_m,
            along_max_m,
            cross_rms_m,
            cross_max_m,
            clock_rms_m,
            clock_max_m,
            clock_datum_removed_rms_m,
            clock_datum_removed_max_m
          ]}
       ) do
    %Stats{
      count: count,
      orbit_3d_rms_m: orbit_3d_rms_m,
      orbit_3d_max_m: orbit_3d_max_m,
      radial_rms_m: radial_rms_m,
      radial_max_m: radial_max_m,
      along_rms_m: along_rms_m,
      along_max_m: along_max_m,
      cross_rms_m: cross_rms_m,
      cross_max_m: cross_max_m,
      clock_rms_m: clock_rms_m,
      clock_max_m: clock_max_m,
      clock_datum_removed_rms_m: clock_datum_removed_rms_m,
      clock_datum_removed_max_m: clock_datum_removed_max_m
    }
  end
end