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`. """ @spec render([map()]) :: iodata() def render(profile) when is_list(profile) do cleaned = clean_profile(profile) [ ~s||, style_block(), isobars(), isotherms(), dry_adiabats(), moist_adiabats(), mixing_ratio_lines(), frame(), axis_labels(), profile_traces(cleaned), cursor_layer(), "" ] 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 """ """ end defp style_block do """ """ 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|| 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|#{p}| end temps = -40..40//10 temp_labels = for t <- temps do x = temperature_x(t, @bottom) if x >= @left and x <= @right do ~s|#{t}| else "" end end [ pres_labels, temp_labels, ~s|Pressure (mb)|, ~s|Temperature (°C)| ] 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|| 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, 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|#{t}°| 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 defp mixing_ratio_lines do # Constant-mixing-ratio lines (g/kg). Each line is the locus of # (T_d, p) for a given saturation mixing ratio. for w <- [0.4, 1, 2, 4, 7, 10, 16, 24, 32] do pts = sample_pressures() |> Enum.filter(&(&1 >= 400)) |> 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() polyline_or_empty(pts, "grid-mix") 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|| 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