prop/lib/microwaveprop_web/skew_t.ex

defmodule MicrowavepropWeb.SkewT do
  @moduledoc """
  Skew-T log-P diagram renderer. Builds the full set of SVG
  primitives (isobars, isotherms, dry adiabats, saturation mixing
  ratio lines, and temperature/dewpoint traces) for a single HRRR
  or NARR pressure-level profile so the contact detail template can
  draw it declaratively.

  ## Meteorology

  * Magnus formula for saturation vapor pressure:
    `es(T) = 6.112 * exp(17.67 * T / (T + 243.5))` with T in °C.
  * Saturation mixing ratio: `ws = 622 * es / (p - es)` in g/kg.
  * Dry adiabat (constant potential temperature θ, K):
    `T(p) = θ * (p / 1000)^(R/cp)` with `R/cp ≈ 0.2854`.

  ## Skew transform

  Pressure is on a log axis (top to bottom, larger p at the bottom).
  Isotherms are skewed 45° by shifting each (T, p) point right as
  pressure decreases, so a vertical line on the chart corresponds
  to a temperature that rises with height — a hallmark of the
  classic skew-T log-P diagram.
  """

  @dew_point_line_color "#16a34a"
  @temp_line_color "#dc2626"

  @dry_adiabat_thetas_k 250..360//10
  @mixing_ratios_g_kg [0.4, 1.0, 2.0, 4.0, 8.0, 12.0, 16.0, 20.0]

  @type projected :: {x :: float(), y :: float()}

  @type chart :: %{
          width: number(),
          height: number(),
          padding_top: number(),
          padding_bottom: number(),
          padding_left: number(),
          padding_right: number(),
          plot_width: number(),
          plot_height: number(),
          t_min: number(),
          t_max: number(),
          p_bot: number(),
          p_top: number(),
          skew: number()
        }

  @default_opts [
    width: 560,
    height: 620,
    padding_top: 24,
    padding_bottom: 36,
    padding_left: 44,
    padding_right: 20,
    t_min: -40.0,
    t_max: 40.0,
    p_bot: 1050.0,
    p_top: 100.0,
    skew: 0.55
  ]

  @doc """
  Build a chart config from overrideable keyword options. Returned
  struct has the plot dimensions pre-computed.
  """
  @spec chart_config(keyword()) :: chart()
  def chart_config(opts \\ []) do
    opts = Keyword.merge(@default_opts, opts)

    width = Keyword.fetch!(opts, :width)
    height = Keyword.fetch!(opts, :height)
    padding_top = Keyword.fetch!(opts, :padding_top)
    padding_bottom = Keyword.fetch!(opts, :padding_bottom)
    padding_left = Keyword.fetch!(opts, :padding_left)
    padding_right = Keyword.fetch!(opts, :padding_right)

    plot_width = width - padding_left - padding_right
    plot_height = height - padding_top - padding_bottom

    %{
      width: width,
      height: height,
      padding_top: padding_top,
      padding_bottom: padding_bottom,
      padding_left: padding_left,
      padding_right: padding_right,
      plot_width: plot_width,
      plot_height: plot_height,
      t_min: Keyword.fetch!(opts, :t_min),
      t_max: Keyword.fetch!(opts, :t_max),
      p_bot: Keyword.fetch!(opts, :p_bot),
      p_top: Keyword.fetch!(opts, :p_top),
      skew: Keyword.fetch!(opts, :skew)
    }
  end

  @doc """
  Build the skew-T primitives for `profile` — a list of maps with
  string keys `"pres"`, `"tmpc"`, `"dwpc"`, and optionally `"hght"`.
  Returns `nil` when the profile is missing or empty.
  """
  @spec build([map()] | nil, keyword()) :: map() | nil
  def build(profile, opts \\ [])
  def build(nil, _), do: nil
  def build([], _), do: nil

  def build(profile, opts) when is_list(profile) do
    chart = chart_config(opts)

    t_points =
      profile
      |> Enum.filter(fn l -> is_number(l["tmpc"]) and is_number(l["pres"]) end)
      |> Enum.sort_by(&(-&1["pres"]))
      |> Enum.map(fn l -> project(l["tmpc"], l["pres"], chart) end)

    td_points =
      profile
      |> Enum.filter(fn l -> is_number(l["dwpc"]) and is_number(l["pres"]) end)
      |> Enum.sort_by(&(-&1["pres"]))
      |> Enum.map(fn l -> project(l["dwpc"], l["pres"], chart) end)

    %{
      width: chart.width,
      height: chart.height,
      padding_top: chart.padding_top,
      padding_bottom: chart.padding_bottom,
      padding_left: chart.padding_left,
      padding_right: chart.padding_right,
      plot_width: chart.plot_width,
      plot_height: chart.plot_height,
      isobars: isobars(chart),
      isotherms: isotherms(chart),
      dry_adiabats: dry_adiabats(chart),
      mixing_ratio_lines: mixing_ratio_lines(chart),
      t_points: t_points,
      td_points: td_points,
      t_path: points_to_path(t_points),
      td_path: points_to_path(td_points),
      temp_color: @temp_line_color,
      dewpoint_color: @dew_point_line_color
    }
  end

  # ── Meteorology ──────────────────────────────────────────────

  @doc """
  Saturation vapor pressure (hPa / mb) for temperature `t_c` in °C.
  Magnus formula with the Tetens coefficients used by IFS/GFS.
  """
  @spec saturation_vapor_pressure(number()) :: float()
  def saturation_vapor_pressure(t_c) when is_number(t_c) do
    6.112 * :math.exp(17.67 * t_c / (t_c + 243.5))
  end

  @doc """
  Invert the saturation-mixing-ratio relation: given a mixing
  ratio `ws_g_kg` and pressure `p_mb`, return the temperature
  (°C) whose saturation vapor pressure yields that mixing ratio.
  """
  @spec temperature_from_mixing_ratio(number(), number()) :: float()
  def temperature_from_mixing_ratio(ws_g_kg, p_mb) when is_number(ws_g_kg) and is_number(p_mb) and ws_g_kg > 0 do
    es = ws_g_kg * p_mb / (622.0 + ws_g_kg)
    ln = :math.log(es / 6.112)
    243.5 * ln / (17.67 - ln)
  end

  @doc """
  Temperature (°C) along the dry adiabat with potential temperature
  `theta_k` (K) at pressure `p_mb` (mb).
  """
  @spec dry_adiabat_temperature(number(), number()) :: float()
  def dry_adiabat_temperature(theta_k, p_mb) when is_number(theta_k) and is_number(p_mb) do
    theta_k * :math.pow(p_mb / 1000.0, 0.2854) - 273.15
  end

  # ── Projection ───────────────────────────────────────────────

  @doc """
  Project a `(t_c, p_mb)` point onto chart pixel coordinates.
  """
  @spec project(number(), number(), chart()) :: projected()
  def project(t_c, p_mb, chart) when is_number(t_c) and is_number(p_mb) do
    p_ratio = pressure_ratio(p_mb, chart)

    y = chart.padding_top + (1.0 - p_ratio) * chart.plot_height

    x_base =
      chart.padding_left +
        (t_c - chart.t_min) / (chart.t_max - chart.t_min) * chart.plot_width

    x = x_base + chart.skew * p_ratio * chart.plot_width
    {x, y}
  end

  defp pressure_ratio(p_mb, chart) do
    # 0 at p_bot (chart bottom), 1 at p_top (chart top). Log-P scale.
    :math.log(chart.p_bot / p_mb) / :math.log(chart.p_bot / chart.p_top)
  end

  # ── Background curves ────────────────────────────────────────

  @isobar_levels_mb [1000, 850, 700, 500, 400, 300, 250, 200, 150, 100]

  defp isobars(chart) do
    for p <- @isobar_levels_mb, p <= chart.p_bot and p >= chart.p_top do
      {_x, y} = project(0.0, p * 1.0, chart)

      %{
        y: y,
        x_start: chart.padding_left,
        x_end: chart.padding_left + chart.plot_width,
        label: "#{p}",
        label_y: y + 3
      }
    end
  end

  defp isotherms(chart) do
    # Draw isotherms every 10°C across the data area. Each one is a
    # straight line between its (t, p_bot) and (t, p_top) projections,
    # clipped to the plot rectangle by the svg viewport.
    for t_c <- -120..40//10 do
      {x_bot, y_bot} = project(t_c * 1.0, chart.p_bot, chart)
      {x_top, y_top} = project(t_c * 1.0, chart.p_top, chart)

      label_x =
        if t_c >= chart.t_min and t_c <= chart.t_max do
          {x, _} = project(t_c * 1.0, chart.p_bot, chart)
          x
        end

      %{
        x1: x_bot,
        y1: y_bot,
        x2: x_top,
        y2: y_top,
        t_c: t_c,
        label_x: label_x,
        label_y: chart.padding_top + chart.plot_height + 14
      }
    end
  end

  defp dry_adiabats(chart) do
    pressures = log_pressure_samples(chart.p_bot, chart.p_top, 24)

    for theta <- @dry_adiabat_thetas_k do
      points =
        Enum.map(pressures, fn p ->
          t = dry_adiabat_temperature(theta * 1.0, p)
          project(t, p, chart)
        end)

      %{theta_k: theta, path: points_to_path(points)}
    end
  end

  defp mixing_ratio_lines(chart) do
    # Mixing ratio isohumes only make sense in the lower troposphere.
    pressures = log_pressure_samples(chart.p_bot, 400.0, 18)

    for ws <- @mixing_ratios_g_kg do
      points =
        Enum.map(pressures, fn p ->
          t = temperature_from_mixing_ratio(ws, p)
          project(t, p, chart)
        end)

      %{ws_g_kg: ws, path: points_to_path(points)}
    end
  end

  defp log_pressure_samples(p_from, p_to, count) do
    log_from = :math.log(p_from)
    log_to = :math.log(p_to)
    step = (log_to - log_from) / (count - 1)

    for i <- 0..(count - 1) do
      :math.exp(log_from + i * step)
    end
  end

  defp points_to_path([]), do: ""

  defp points_to_path([{x0, y0} | rest]) do
    initial = "M#{format_coord(x0)},#{format_coord(y0)}"

    body =
      Enum.map_join(rest, " ", fn {x, y} -> "L#{format_coord(x)},#{format_coord(y)}" end)

    if body == "", do: initial, else: initial <> " " <> body
  end

  defp format_coord(n) when is_number(n), do: :erlang.float_to_binary(n * 1.0, decimals: 2)
end