prop/lib/microwaveprop_web/live/skewt_svg.ex

defmodule MicrowavepropWeb.SkewtSvg do
  @moduledoc """
  Render an HRRR vertical profile as a Skew-T-Log-P diagram in SVG.

  HRRR's persisted profile stores `pres`, `hght`, `tmpc`, `dwpc` per
  pressure level — same canonical shape `Microwaveprop.Weather.SoundingParams`
  consumes. This renderer ignores anything else.

  Wind barbs are deliberately omitted: HRRR persists wind only at 10 m
  AGL, not on the pressure-level profile, so a per-level barb column
  would be misleading.
  """

  # Plot canvas (SVG user units).
  @width 720
  @height 640
  @left 60
  @right 660
  @top 30
  @bottom 580

  @p_bottom 1050.0
  @p_top 100.0
  @t_min -40.0
  @t_max 40.0

  # Skew factor — pixels of horizontal offset added per pixel of
  # height. 1.0 puts isotherms at 45°.
  @skew 1.0

  @doc """
  Render a Skew-T diagram for `profile` (a list of `%{"pres", "hght",
  "tmpc", "dwpc"}` maps) and return an inline SVG string suitable for
  `Phoenix.HTML.raw/1`.

  When `:parcel_trace` is supplied (typically from
  `Microwaveprop.Weather.SkewtParams.derive/1`) a dashed grey parcel-
  ascent line is overlaid on the temperature trace, matching the
  surface-parcel rendition the SPC sounding viewer publishes.
  """
  @spec render([map()], keyword()) :: iodata()
  def render(profile, opts \\ []) when is_list(profile) do
    cleaned = clean_profile(profile)
    parcel = Keyword.get(opts, :parcel_trace, [])

    [
      ~s|<svg viewBox="0 0 #{@width} #{@height}" xmlns="http://www.w3.org/2000/svg" class="w-full h-auto" role="img" aria-label="Skew-T-Log-P diagram">|,
      style_block(),
      isobars(),
      isotherms(),
      dry_adiabats(),
      moist_adiabats(),
      mixing_ratio_lines(),
      frame(),
      axis_labels(),
      profile_traces(cleaned),
      parcel_trace_line(parcel),
      cursor_layer(),
      "</svg>"
    ]
  end

  # Surface-parcel ascent line (dashed grey), matching the SPC viewer's
  # convention. The trace is a list of `%{pres, tmpc}` from surface up
  # through every level above the LCL.
  defp parcel_trace_line([]), do: ""

  defp parcel_trace_line(trace) do
    pts =
      trace
      |> Enum.sort_by(& &1.pres, :desc)
      |> Enum.map(fn %{pres: p, tmpc: t} ->
        y = pressure_y(p)
        {temperature_x(t, y), y}
      end)
      |> clip_polyline()

    case pts do
      [] -> ""
      _ -> ~s|<polyline class="parcel-line" points="#{points_str(pts)}" />|
    end
  end

  defp points_str(points) do
    Enum.map_join(points, " ", fn {x, y} -> "#{r(x)},#{r(y)}" end)
  end

  # Crosshair layer the JS hook (`assets/js/skewt_hook.js`) drives.
  # Initially `display: none`; the hook flips display + updates element
  # coordinates on `mousemove`.
  defp cursor_layer do
    """
    <g class="skewt-cursor" data-cursor style="display:none">
      <line class="skewt-cursor-line" x1="#{@left}" x2="#{@right}" y1="0" y2="0"/>
      <circle class="skewt-cursor-dot-t" data-cursor-dot-t cx="0" cy="0" r="3"/>
      <circle class="skewt-cursor-dot-d" data-cursor-dot-d cx="0" cy="0" r="3"/>
      <rect class="skewt-cursor-readout-bg" data-cursor-readout-bg x="0" y="0" width="148" height="62" rx="3"/>
      <text class="skewt-cursor-readout" data-cursor-readout x="0" y="0">
        <tspan data-cursor-line-1 x="0" dy="14"></tspan>
        <tspan data-cursor-line-2 x="0" dy="14"></tspan>
        <tspan data-cursor-line-3 x="0" dy="14"></tspan>
        <tspan data-cursor-line-4 x="0" dy="14"></tspan>
      </text>
    </g>
    """
  end

  # NOAA SPC Skew-T-Log-P colour set, copied from the SPC sounding viewer
  # (https://www.spc.noaa.gov/exper/soundings/) so the chart reads the
  # same way as the public reference product:
  #
  #   * Temperature trace .................. solid red (#ff0000)
  #   * Dewpoint trace ..................... solid green (#00aa00)
  #   * Isotherms .......................... dashed teal (#22d3a4-ish)
  #   * Dry adiabats ....................... solid orange (#f59e0b)
  #   * Moist (saturated) adiabats ......... dashed sky-blue (#38bdf8)
  #   * Mixing-ratio lines ................. dotted darker green (#16a34a),
  #                                          values 0.4–20 g/kg, drawn only
  #                                          in the lower atmosphere
  #   * Isobars ............................ light grey horizontal lines
  defp style_block do
    """
    <defs>
      <style>
        .frame { fill: none; stroke: currentColor; stroke-width: 1.2; }
        .grid-major { fill: none; stroke: currentColor; stroke-opacity: 0.45; stroke-width: 0.6; }
        .grid-minor { fill: none; stroke: currentColor; stroke-opacity: 0.18; stroke-width: 0.4; }
        .grid-iso { fill: none; stroke: #14b8a6; stroke-opacity: 0.55; stroke-width: 0.5; stroke-dasharray: 4 3; }
        .grid-dry { fill: none; stroke: #f59e0b; stroke-opacity: 0.65; stroke-width: 0.5; }
        .grid-moist { fill: none; stroke: #38bdf8; stroke-opacity: 0.55; stroke-width: 0.5; stroke-dasharray: 4 3; }
        .grid-mix { fill: none; stroke: #16a34a; stroke-opacity: 0.55; stroke-width: 0.5; stroke-dasharray: 1 3; }
        .axis { font: 10px ui-sans-serif, system-ui, sans-serif; fill: currentColor; }
        .iso-label { font: 9px ui-sans-serif, system-ui, sans-serif; fill: #14b8a6; fill-opacity: 0.85; }
        .mix-label { font: 9px ui-sans-serif, system-ui, sans-serif; fill: #16a34a; fill-opacity: 0.85; }
        .pres-label { font: 9px ui-sans-serif, system-ui, sans-serif; fill: currentColor; fill-opacity: 0.7; }
        .temp-line { fill: none; stroke: #ff0000; stroke-width: 2.5; stroke-linejoin: round; stroke-linecap: round; }
        .dew-line  { fill: none; stroke: #00aa00; stroke-width: 2.5; stroke-linejoin: round; stroke-linecap: round; }
        .parcel-line { fill: none; stroke: #475569; stroke-width: 1.5; stroke-dasharray: 5 4; stroke-linejoin: round; stroke-linecap: round; }
        .skewt-cursor-line { stroke: currentColor; stroke-opacity: 0.5; stroke-width: 0.8; stroke-dasharray: 3 3; pointer-events: none; }
        .skewt-cursor-dot-t { fill: #ff0000; pointer-events: none; }
        .skewt-cursor-dot-d { fill: #00aa00; pointer-events: none; }
        .skewt-cursor-readout-bg { fill: rgba(15,23,42,0.85); pointer-events: none; }
        .skewt-cursor-readout { font: 11px ui-monospace, SFMono-Regular, monospace; fill: #f8fafc; pointer-events: none; }
      </style>
    </defs>
    """
  end

  # ── Geometry exposed to the hover hook ──────────────────────────────

  @doc """
  Plot geometry constants the JS hover hook needs to invert pixel-Y
  back to pressure and to project temperature values onto the skewed
  diagram. Returned as a plain map so it serializes to JSON via
  `Jason.encode!/1` straight onto a `data-` attribute.
  """
  @spec geometry() :: map()
  def geometry do
    %{
      width: @width,
      height: @height,
      left: @left,
      right: @right,
      top: @top,
      bottom: @bottom,
      p_bottom: @p_bottom,
      p_top: @p_top,
      t_min: @t_min,
      t_max: @t_max,
      skew: @skew
    }
  end

  # ── Coordinate transforms ───────────────────────────────────────────

  @doc false
  def pressure_y(p) when is_number(p) do
    # High pressure → bottom of plot, low pressure → top.
    @bottom -
      (@bottom - @top) * (:math.log(@p_bottom) - :math.log(p)) /
        (:math.log(@p_bottom) - :math.log(@p_top))
  end

  @doc false
  def temperature_x(t, y) when is_number(t) and is_number(y) do
    plot_w = @right - @left
    base_x = @left + (t - @t_min) / (@t_max - @t_min) * plot_w
    base_x + @skew * (@bottom - y)
  end

  defp clamp_x(x), do: x |> max(@left) |> min(@right)

  defp clip_segment({x1, y1}, {x2, y2}) do
    cond do
      x1 < @left and x2 < @left -> nil
      x1 > @right and x2 > @right -> nil
      true -> {{clamp_x(x1), y1}, {clamp_x(x2), y2}}
    end
  end

  # ── Frame & axes ────────────────────────────────────────────────────

  defp frame do
    ~s|<rect class="frame" x="#{@left}" y="#{@top}" width="#{@right - @left}" height="#{@bottom - @top}" />|
  end

  defp axis_labels do
    pressures = [1000, 850, 700, 500, 400, 300, 200, 100]

    pres_labels =
      for p <- pressures do
        y = pressure_y(p)
        ~s|<text class="pres-label" x="#{@left - 6}" y="#{y + 3}" text-anchor="end">#{p}</text>|
      end

    temps = -40..40//10

    temp_labels =
      for t <- temps do
        x = temperature_x(t, @bottom)

        if x >= @left and x <= @right do
          ~s|<text class="axis" x="#{x}" y="#{@bottom + 14}" text-anchor="middle">#{t}</text>|
        else
          ""
        end
      end

    [
      pres_labels,
      temp_labels,
      ~s|<text class="axis" x="#{@left - 38}" y="#{(@top + @bottom) / 2}" text-anchor="middle" transform="rotate(-90 #{@left - 38} #{(@top + @bottom) / 2})">Pressure (mb)</text>|,
      ~s|<text class="axis" x="#{(@left + @right) / 2}" y="#{@bottom + 30}" text-anchor="middle">Temperature (°C)</text>|
    ]
  end

  # ── Background grid ─────────────────────────────────────────────────

  defp isobars do
    pressures = [1000, 925, 850, 700, 500, 400, 300, 250, 200, 150, 100]

    for p <- pressures do
      y = pressure_y(p)
      cls = if p in [1000, 850, 700, 500, 300, 200, 100], do: "grid-major", else: "grid-minor"
      ~s|<line class="#{cls}" x1="#{@left}" y1="#{y}" x2="#{@right}" y2="#{y}" />|
    end
  end

  defp isotherms do
    # Skewed isotherms every 10°C across the plot. Each line goes from
    # (T at bottom) to (T at top), with skew already baked into
    # temperature_x. Clip to the plot box.
    for t <- -100..40//10 do
      render_isotherm(t)
    end
  end

  defp render_isotherm(t) do
    p1 = {temperature_x(t, @bottom), @bottom}
    p2 = {temperature_x(t, @top), @top}

    case clip_segment(p1, p2) do
      nil -> ""
      clipped -> isotherm_with_label(t, clipped)
    end
  end

  defp isotherm_with_label(t, {{x1, y1}, {x2, y2}}) do
    line = ~s|<line class="grid-iso" x1="#{r(x1)}" y1="#{r(y1)}" x2="#{r(x2)}" y2="#{r(y2)}" />|
    [line, isotherm_label(t, x1)]
  end

  defp isotherm_label(t, x1) when t in [-40, -20, 0, 20, 40] do
    if x1 > @left + 2 and x1 < @right - 12 do
      ~s|<text class="iso-label" x="#{r(x1) + 2}" y="#{@bottom - 2}">#{t}°</text>|
    else
      ""
    end
  end

  defp isotherm_label(_t, _x1), do: ""

  defp dry_adiabats do
    # Potential-temperature curves: T_K(p) = θ * (p / 1000)^(R/cp).
    # Sample every theta = -30..150°C every 10° at 50-mb pressure steps.
    for theta_c <- -30..150//10 do
      theta_k = theta_c + 273.15

      pts =
        sample_pressures()
        |> Enum.map(fn p ->
          t_k = theta_k * :math.pow(p / 1000.0, 0.2854)
          t_c = t_k - 273.15
          y = pressure_y(p)
          x = temperature_x(t_c, y)
          {x, y}
        end)
        |> clip_polyline()

      polyline_or_empty(pts, "grid-dry")
    end
  end

  defp moist_adiabats do
    # Saturated pseudo-adiabats — quick approximation by integrating
    # dT/dp via the moist-adiabatic lapse rate. Good enough to put
    # green dashes in the right place; not used for parcel ascent.
    for theta_e_c <- 0..40//4 do
      pts =
        (theta_e_c + 273.15)
        |> moist_adiabat_curve()
        |> Enum.map(fn {p, t_c} ->
          y = pressure_y(p)
          x = temperature_x(t_c, y)
          {x, y}
        end)
        |> clip_polyline()

      polyline_or_empty(pts, "grid-moist")
    end
  end

  defp moist_adiabat_curve(theta_e_k) do
    # Integrate from 1000 mb upward using dT/dz = Γ_m and dp/dz from
    # hydrostatic. Cheap approximation good enough for plotting.
    starting_t = theta_e_k - 273.15

    sample_pressures()
    |> Enum.reduce({[], starting_t, 1000.0}, fn p, {acc, t, prev_p} ->
      dp = p - prev_p
      # Δz from hydrostatic with mean T_K ≈ t + 273
      mean_t_k = max(t + 273.15, 200.0)
      dz = -dp * 287.0 * mean_t_k / (9.81 * p)
      gamma_m = moist_lapse_c_per_m(t, p)
      new_t = t + dz * gamma_m
      {[{p, new_t} | acc], new_t, p}
    end)
    |> elem(0)
    |> Enum.reverse()
  end

  defp moist_lapse_c_per_m(t_c, p_mb) do
    t_k = t_c + 273.15
    # Saturation mixing ratio (g/g) via Buck's eq.
    es = 6.1121 * :math.exp((18.678 - t_c / 234.5) * (t_c / (257.14 + t_c)))
    qs = 0.622 * es / (p_mb - es)
    lv = 2.5e6
    cp = 1004.0
    rd = 287.0
    rv = 461.0
    g = 9.81
    num = g * (1.0 + lv * qs / (rd * t_k))
    den = cp + lv * lv * qs / (rv * t_k * t_k)
    -num / den
  end

  # NOAA SPC mixing-ratio set, in g/kg. Same values quoted along the
  # bottom of the SPC sounding charts (0.4 / 0.7 / 1 / 2 / 3 / 5 / 8 /
  # 12 / 16 / 20). Lines drawn from the surface up to 600 mb only —
  # above that the saturation values get unhelpfully close together.
  @mixing_ratios [0.4, 0.7, 1, 2, 3, 5, 8, 12, 16, 20]
  @mix_top_pres 600.0
  @mix_label_pres 700.0

  defp mixing_ratio_lines do
    Enum.map(@mixing_ratios, fn w ->
      pts =
        sample_pressures()
        |> Enum.filter(&(&1 >= @mix_top_pres))
        |> Enum.map(fn p ->
          # Solve for T from w = 0.622 * e_s(T) / (p - e_s).
          e = w / 1000.0 * p / (0.622 + w / 1000.0)
          t_d = dewpoint_from_vapor(e)
          y = pressure_y(p)
          x = temperature_x(t_d, y)
          {x, y}
        end)
        |> clip_polyline()

      label = mix_label(w)
      [polyline_or_empty(pts, "grid-mix"), label]
    end)
  end

  # Anchor each mixing-ratio line with its g/kg value at ~700 mb,
  # mirroring SPC's row of green numbers along the diagram's lower band.
  defp mix_label(w) do
    e = w / 1000.0 * @mix_label_pres / (0.622 + w / 1000.0)
    t_d = dewpoint_from_vapor(e)
    y = pressure_y(@mix_label_pres)
    x = temperature_x(t_d, y)

    if x >= @left and x <= @right do
      label =
        if w >= 1 do
          w |> trunc() |> Integer.to_string()
        else
          :erlang.float_to_binary(w * 1.0, decimals: 1)
        end

      ~s|<text class="mix-label" x="#{r(x)}" y="#{r(y - 2)}" text-anchor="middle">#{label}</text>|
    else
      ""
    end
  end

  # Inverse Buck — find T given vapor pressure e (hPa).
  defp dewpoint_from_vapor(e) when e <= 0, do: -80.0

  defp dewpoint_from_vapor(e) do
    # Magnus form: T = 243.5 * ln(e/6.112) / (17.67 - ln(e/6.112))
    ratio = :math.log(e / 6.112)
    243.5 * ratio / (17.67 - ratio)
  end

  defp sample_pressures do
    # 1050 down to 100 in 25-mb steps.
    100..1050//25 |> Enum.to_list() |> Enum.reverse()
  end

  defp clip_polyline(points) do
    Enum.filter(points, fn {x, y} -> x >= @left - 200 and x <= @right + 200 and y >= @top and y <= @bottom end)
  end

  defp polyline_or_empty([], _cls), do: ""

  defp polyline_or_empty(points, cls) do
    pts =
      Enum.map_join(points, " ", fn {x, y} -> "#{r(x)},#{r(y)}" end)

    ~s|<polyline class="#{cls}" points="#{pts}" />|
  end

  # ── Profile traces ──────────────────────────────────────────────────

  defp profile_traces([]), do: ""

  defp profile_traces(profile) do
    sorted = Enum.sort_by(profile, & &1.pres, :desc)

    temp_pts =
      sorted
      |> Enum.map(fn p ->
        y = pressure_y(p.pres)
        x = temperature_x(p.tmpc, y)
        {x, y}
      end)
      |> Enum.filter(fn {_x, y} -> y >= @top and y <= @bottom end)

    dew_pts =
      sorted
      |> Enum.filter(& &1.dwpc)
      |> Enum.map(fn p ->
        y = pressure_y(p.pres)
        x = temperature_x(p.dwpc, y)
        {x, y}
      end)
      |> Enum.filter(fn {_x, y} -> y >= @top and y <= @bottom end)

    [
      polyline_or_empty(dew_pts, "dew-line"),
      polyline_or_empty(temp_pts, "temp-line")
    ]
  end

  defp clean_profile(profile) do
    profile
    |> Enum.map(&normalize_level/1)
    |> Enum.filter(&(&1.pres && &1.tmpc))
  end

  defp normalize_level(%{} = level) do
    %{
      pres: pick_num(level, ["pres", "pres_mb", :pres, :pres_mb]),
      hght: pick_num(level, ["hght", "hght_m", :hght, :hght_m]),
      tmpc: pick_num(level, ["tmpc", :tmpc]),
      dwpc: pick_num(level, ["dwpc", :dwpc])
    }
  end

  defp pick_num(map, keys) do
    keys
    |> Enum.find_value(fn k -> Map.get(map, k) end)
    |> coerce_num()
  end

  defp coerce_num(n) when is_number(n), do: n * 1.0
  defp coerce_num(_), do: nil

  defp r(n) when is_float(n), do: Float.round(n, 1)
  defp r(n), do: n
end