space_dust

defmodule SpaceDust.Data.EarthOrientationParameters do
  @moduledoc """
  Parameters required to compute the CRS-TRS transformation matrix
  """

  @type t :: %__MODULE__{
          modifiedJulianDate: float() | nil,
          polarMotionX: float() | nil,
          polarMotionY: float() | nil,
          ut1UTC: float() | nil,
          lod: float() | nil,
          dEps: float() | nil,
          dPsi: float() | nil
        }

  defstruct [
    :modifiedJulianDate,
    :polarMotionX,
    :polarMotionY,
    :ut1UTC,
    :lod,
    :dEps,
    :dPsi
  ]
end

defmodule SpaceDust.Data.EOP do
  @moduledoc """
  Earth Orientation Parameters (EOP) data - ref https://hpiers.obspm.fr/eoppc/bul/bulb/explanatory.html

  Uses the IERS finals.all data which includes both historical observations (marked 'I')
  and predictions up to ~1 year into the future (marked 'P').
  """

  # finals.all contains historical + predicted EOP data (predictions ~1 year ahead)
  @eopDataUrl "https://datacenter.iers.org/data/latestVersion/finals.all.iau2000.txt"
  rootPath = Path.expand("../..", __DIR__)
  @dataPath rootPath <> "/data"
  @eopDataPath @dataPath <> "/eop_data.txt"

  case File.exists?(@dataPath) do
    true -> IO.puts("Data directory exists")
    false -> File.mkdir(@dataPath)
  end

  require Logger
  alias SpaceDust.Data.EarthOrientationParameters, as: EarthOrientationParameters
  alias SpaceDust.Math.Functions, as: Functions

  @spec pullEOPData() :: {:error, String.t()} | {:ok, list()}
  @doc "pull EOP data from the IERS - realistically only needed once a year, or at container start"
  def pullEOPData() do
    response = Req.get!(@eopDataUrl)

    case response.status do
      200 ->
        # split the response body into lines
        lines =
          String.split(response.body, "\n")
          |> Enum.map(&String.trim/1)
          # remove empty lines
          |> Enum.reject(&(&1 == ""))
          # remove comment lines
          |> Enum.reject(&String.starts_with?(&1, "#"))

        {:ok, lines}

      _ ->
        {:error, "Failed to retrieve EOP data from IERS with status #{response.status}"}
    end
  end

  @spec saveEopData(any()) :: {:error, atom()}
  @doc "save EOP data to a file"
  def saveEopData(lines) do
    # TODO: make a data directory if it doesn't exist
    case File.write(@eopDataPath, Enum.join(lines, "\n")) do
      :ok -> {:ok, @eopDataPath}
      {:error, reason} -> {:error, reason}
    end
  end

  @spec readSavedEopData() :: {:error, atom()} | {:ok, [binary()]}
  def readSavedEopData() do
    case File.read("eop_data.txt") do
      {:ok, data} ->
        {:ok, String.split(data, "\n")}

      {:error, :enoent} ->
        case pullEOPData() do
          {:ok, lines} ->
            saveEopData(lines)
            {:ok, lines}

          {:error, reason} ->
            {:error, reason}
        end

      {:error, reason} ->
        {:error, reason}
    end
  end

  def parseEopLine(line) do
    # finals.all format parsing
    # The date fields (year, month, day) at the start can vary in format,
    # but the MJD and subsequent fields are space-separated in a consistent pattern.
    #
    # Note: Sometimes the I/P flag is concatenated with negative values like "I-0.4575789"
    # We pre-process to add spaces before I/P flags to handle this.
    #
    # After MJD, the pattern is:
    #   [MJD] [flag] [PM-x] [PM-x err] [PM-y] [PM-y err] [flag] [UT1-UTC] [UT1-UTC err] [LOD] [LOD err] [flag] [dPSI] [dPSI err] [dEPS] [dEPS err] ...

    try do
      # Pre-process: Add space before I or P that's followed by a digit or minus sign
      # This handles cases like "I-0.4575789" -> "I -0.4575789"
      normalized_line = line
        |> String.replace(~r/([IP])(-?\d)/, "\\1 \\2")

      # Split by whitespace
      parts = String.split(normalized_line, ~r/\s+/, trim: true)

      # Find the MJD field (matches pattern like 41684.00 or 60681.00)
      # Must contain decimal point to distinguish from date fields like "41010"
      mjd_idx = Enum.find_index(parts, fn p ->
        String.contains?(p, ".") and
        case Float.parse(p) do
          {val, ""} -> val > 40000 and val < 70000  # MJD range: ~1970 to ~2050
          _ -> false
        end
      end)

      if mjd_idx == nil or length(parts) < mjd_idx + 9 do
        {:error, "Cannot find MJD or line too short"}
      else
        # Parse fields relative to MJD position
        mjd = parse_float(Enum.at(parts, mjd_idx))
        # Skip flag at mjd_idx + 1
        pm_x = parse_float(Enum.at(parts, mjd_idx + 2))
        # Skip PM-x error at mjd_idx + 3
        pm_y = parse_float(Enum.at(parts, mjd_idx + 4))
        # Skip PM-y error at mjd_idx + 5
        # Skip flag at mjd_idx + 6
        ut1_utc = parse_float(Enum.at(parts, mjd_idx + 7))
        # Skip UT1-UTC error at mjd_idx + 8

        # LOD is at mjd_idx + 9, but may be missing or be the nutation flag
        lod_raw = Enum.at(parts, mjd_idx + 9)
        {lod, next_offset} = case lod_raw do
          nil -> {0.0, 10}
          "I" -> {0.0, 9}   # This is actually the nutation flag
          "P" -> {0.0, 9}   # This is actually the nutation flag
          str ->
            case parse_float(str) do
              nil -> {0.0, 10}
              val -> {val / 1000.0, 11}  # Convert ms to s, skip LOD error
            end
        end

        # Find nutation values - look for next I/P flag after LOD
        nutation_flag_idx = mjd_idx + next_offset
        flag_val = Enum.at(parts, nutation_flag_idx)

        {dpsi, deps} = if flag_val in ["I", "P"] and length(parts) > nutation_flag_idx + 4 do
          dpsi_raw = parse_float(Enum.at(parts, nutation_flag_idx + 1))
          # Skip dPSI error at nutation_flag_idx + 2
          deps_raw = parse_float(Enum.at(parts, nutation_flag_idx + 3))

          dpsi = case dpsi_raw do
            nil -> 0.0
            val -> val / 1000.0  # Convert mas to arcsec
          end

          deps = case deps_raw do
            nil -> 0.0
            val -> val / 1000.0  # Convert mas to arcsec
          end

          {dpsi, deps}
        else
          {0.0, 0.0}
        end

        if mjd != nil and pm_x != nil and pm_y != nil and ut1_utc != nil do
          {:ok,
           %EarthOrientationParameters{
             modifiedJulianDate: mjd,
             polarMotionX: pm_x,
             polarMotionY: pm_y,
             ut1UTC: ut1_utc,
             dPsi: dpsi,
             dEps: deps,
             lod: lod
           }}
        else
          {:error, "Unable to parse required EOP fields"}
        end
      end
    rescue
      _ -> {:error, "Unable to parse EOP line"}
    end
  end

  # Helper to safely parse floats, returning nil for empty/invalid strings
  defp parse_float(""), do: nil
  defp parse_float(str) do
    case Float.parse(str) do
      {val, _} -> val
      :error -> nil
    end
  end

  @doc "parse an entire EOP data file"
  def parseEopData(lines) do
    eopData =
      lines
      |> Enum.map(&parseEopLine/1)
      |> Enum.filter(fn
        {:ok, _eop} -> true
        _ -> false
      end)
      |> Enum.map(fn {:ok, eop} -> eop end)

    # check if we have any EOP data
    case Enum.empty?(eopData) do
      true -> {:error, "No EOP data found"}
      false -> {:ok, eopData}
    end
  end

  @doc "interpolate between two EOP data points at a given MJD"
  def interpolateBetween(priorEop, futureEop, mjd) do
    # calculate the fraction of the time between the two EOP data points
    fraction =
      (mjd - priorEop.modifiedJulianDate) /
        (futureEop.modifiedJulianDate - priorEop.modifiedJulianDate)

    # interpolate the EOP data
    %EarthOrientationParameters{
      modifiedJulianDate: mjd,
      polarMotionX:
        Functions.linearInterpolate(priorEop.polarMotionX, futureEop.polarMotionX, fraction),
      polarMotionY:
        Functions.linearInterpolate(priorEop.polarMotionY, futureEop.polarMotionY, fraction),
      ut1UTC: Functions.linearInterpolate(priorEop.ut1UTC, futureEop.ut1UTC, fraction),
      dPsi: Functions.linearInterpolate(priorEop.dPsi, futureEop.dPsi, fraction),
      dEps: Functions.linearInterpolate(priorEop.dEps, futureEop.dEps, fraction),
      lod: Functions.linearInterpolate(priorEop.lod, futureEop.lod, fraction)
    }
  end

  @doc "get the EOP data at a given MJD, interpolating between the closest two data points"
  def getEopData(mjd) do
    case readSavedEopData() do
      {:ok, rawEopData} ->
        case parseEopData(rawEopData) do
          {:ok, eopData} ->
            # sort
            sortedData = Enum.sort_by(eopData, & &1.modifiedJulianDate)

            Enum.find(sortedData, fn eop -> eop.modifiedJulianDate >= mjd end)
            |> case do
              nil ->
                Logger.warning("MJD #{mjd} is after the last EOP data point")
                {:ok, Enum.at(sortedData, -1)}

              %EarthOrientationParameters{} = futureEop ->
                Enum.find(sortedData, fn eop -> eop.modifiedJulianDate < mjd end)
                |> case do
                  nil ->
                    Logger.warning("MJD #{mjd} is before the first EOP data point")
                    {:ok, Enum.at(sortedData, 0)}

                  %EarthOrientationParameters{} = priorEop ->
                    {:ok, interpolateBetween(priorEop, futureEop, mjd)}

                  _ ->
                    {:error, "Unable to find EOP data"}
                end

              _ ->
                {:error, "Unable to find EOP data"}
            end

          {:error, reason} ->
            Logger.error(reason)
            {:error, reason}
        end

      {:error, reason} ->
        IO.inspect(reason)
    end
  end
end