diff --git a/lib/microwaveprop/weather/skewt_params.ex b/lib/microwaveprop/weather/skewt_params.ex new file mode 100644 index 00000000..8c796bfb --- /dev/null +++ b/lib/microwaveprop/weather/skewt_params.ex @@ -0,0 +1,533 @@ +defmodule Microwaveprop.Weather.SkewtParams do + @moduledoc """ + Derived sounding parameters for the `/skewt` page — the same set the + SPC sounding viewer publishes in its `*.txt` companion file + (`https://www.spc.noaa.gov/exper/soundings/`): + + * Lifting condensation level (`lcl_pressure_mb`, `lcl_temp_c`, + `lcl_height_m`). + * Level of free convection (`lfc_pressure_mb`, `lfc_height_m`). + * Equilibrium level (`el_pressure_mb`, `el_height_m`). + * Surface-based CAPE/CIN (`sbcape`, `sbcin` in J/kg). + * Three-kilometre CAPE (`cape_3km`). + * Lifted index (`lifted_index`, parcel-environment T diff at 500 mb). + * Freezing level and wet-bulb-zero level (`freezing_level_m`, + `wbz_m`). + * Downdraft CAPE (`dcape`). + * Layer lapse rates (`lapse_rate_sfc_3km_c_per_km`, + `lapse_rate_700_500_c_per_km`, `lapse_rate_850_500_c_per_km`). + * `parcel_trace` — the surface parcel's pressure/temperature + trajectory through the diagram, used by `SkewtSvg.render/1` to + draw the dashed parcel-ascent overlay. + + Wind-derived indices (SRH, BWD, Bunkers motion, SCP, STP, SHIP) are + *not* produced because HRRR persists wind only at 10 m AGL — there + is no per-pressure-level wind to compute them from. Operators wanting + shear-based products should use NUCAPS / SPC mesoanalysis instead. + + Inputs are the same canonical profile shape `SoundingParams.derive/1` + consumes: `[%{"pres", "hght", "tmpc", "dwpc"}, ...]` (string OR atom + keys; a `:hght_m` / `"hght_m"` alias for height is also accepted). + """ + + @r_d 287.05 + @r_v 461.5 + @c_p 1004.0 + @g 9.80665 + @epsilon @r_d / @r_v + # Buck's saturation-vapor formula constants (kept in this module so + # we don't have to leak them out of SoundingParams). + @buck_a 6.1121 + @buck_b 17.502 + @buck_c 240.97 + # Approximate L_v at 0 °C in J/kg. Tropospheric variation across + # parcel temperatures is ≤4 % which is well below CAPE bin noise. + @l_v 2.5e6 + + @nil_params %{ + lcl_pressure_mb: nil, + lcl_temp_c: nil, + lcl_height_m: nil, + lfc_pressure_mb: nil, + lfc_height_m: nil, + el_pressure_mb: nil, + el_height_m: nil, + sbcape: nil, + sbcin: nil, + cape_3km: nil, + lifted_index: nil, + freezing_level_m: nil, + wbz_m: nil, + dcape: nil, + lapse_rate_sfc_3km_c_per_km: nil, + lapse_rate_700_500_c_per_km: nil, + lapse_rate_850_500_c_per_km: nil, + parcel_trace: [] + } + + @doc """ + Returns the SPC-style parameter map for `profile`. See module + docstring for field-by-field semantics. Returns `@nil_params` for + empty / single-level inputs. + """ + @spec derive([map()]) :: map() + def derive(profile) when is_list(profile) do + levels = profile |> Enum.map(&normalize_level/1) |> Enum.reject(&is_nil/1) + + case Enum.sort_by(levels, & &1.pres, :desc) do + [] -> @nil_params + [_only] -> @nil_params + [surface | _] = sorted -> compute(surface, sorted) + end + end + + defp compute(surface, sorted_desc) do + {p_lcl, t_lcl} = lcl(surface.tmpc, surface.dwpc, surface.pres) + parcel_trace = parcel_trace(surface, p_lcl, t_lcl, sorted_desc) + {sbcape, sbcin} = cape_cin(parcel_trace, sorted_desc) + cape_3km = cape_layer(parcel_trace, sorted_desc, surface.hght, surface.hght + 3000.0) + {lfc_pressure, el_pressure} = lfc_el(parcel_trace, sorted_desc) + + {lfc_height, el_height} = + {height_at_pressure(sorted_desc, lfc_pressure), height_at_pressure(sorted_desc, el_pressure)} + + %{ + lcl_pressure_mb: p_lcl, + lcl_temp_c: t_lcl, + lcl_height_m: height_at_pressure(sorted_desc, p_lcl), + lfc_pressure_mb: lfc_pressure, + lfc_height_m: lfc_height, + el_pressure_mb: el_pressure, + el_height_m: el_height, + sbcape: sbcape, + sbcin: sbcin, + cape_3km: cape_3km, + lifted_index: lifted_index(parcel_trace, sorted_desc), + freezing_level_m: temperature_height_crossing(sorted_desc, & &1.tmpc, 0.0), + wbz_m: temperature_height_crossing(sorted_desc, &wet_bulb_c/1, 0.0), + dcape: dcape(sorted_desc), + lapse_rate_sfc_3km_c_per_km: lapse_rate_height(sorted_desc, surface.hght, surface.hght + 3000.0), + lapse_rate_700_500_c_per_km: lapse_rate_pressures(sorted_desc, 700.0, 500.0), + lapse_rate_850_500_c_per_km: lapse_rate_pressures(sorted_desc, 850.0, 500.0), + parcel_trace: parcel_trace + } + end + + # ── Lifting Condensation Level (Bolton 1980 eq 22) ────────────────── + + @doc """ + Lifting condensation level pressure and temperature. + + Inputs: + * `temp_c` — surface temperature, °C + * `dewpoint_c` — surface dewpoint, °C + * `pressure_mb` — surface pressure, mb + + Returns `{p_lcl_mb, t_lcl_c}`. Saturated air (T = T_d) returns the + surface pressure unchanged. + """ + @spec lcl(number(), number(), number()) :: {float(), float()} + def lcl(temp_c, dewpoint_c, pressure_mb) when is_number(temp_c) and is_number(dewpoint_c) and is_number(pressure_mb) do + t_k = temp_c + 273.15 + td_k = dewpoint_c + 273.15 + + # Bolton (1980) eq 22: T at LCL. + t_lcl_k = 1.0 / (1.0 / (td_k - 56.0) + :math.log(t_k / td_k) / 800.0) + 56.0 + + # Poisson relation between two adiabatic levels. + p_lcl = pressure_mb * :math.pow(t_lcl_k / t_k, @c_p / @r_d) + + {p_lcl, t_lcl_k - 273.15} + end + + # ── Parcel ascent ─────────────────────────────────────────────────── + + # Returns a list of `%{pres, tmpc}` for the surface parcel, walking + # from surface pressure up through every environmental level above + # the LCL using a saturated-adiabatic step in pressure. + defp parcel_trace(surface, p_lcl, t_lcl, sorted_desc) do + surface_step = %{pres: surface.pres, tmpc: surface.tmpc} + lcl_step = %{pres: min(p_lcl, surface.pres), tmpc: t_lcl} + + above_lcl = + sorted_desc + |> Enum.filter(&(&1.pres < p_lcl)) + |> Enum.sort_by(& &1.pres, :desc) + + {trace_above, _} = + Enum.map_reduce(above_lcl, lcl_step, fn level, prev -> + t_next = saturated_adiabatic_step(prev.tmpc, prev.pres, level.pres) + step = %{pres: level.pres, tmpc: t_next} + {step, step} + end) + + Enum.uniq_by([surface_step, lcl_step | trace_above], & &1.pres) + end + + # Integrate dT/d(ln p) = γ_sat × R_d × T / g across the pressure + # interval [p1, p2] in 5 mb steps using a forward-Euler scheme. + # Coarse but stable; the residual error is well below CAPE-bin noise. + defp saturated_adiabatic_step(t_c_start, p_start, p_end) do + n_steps = max(1, trunc((p_start - p_end) / 5.0)) + dp = (p_end - p_start) / n_steps + + 1..n_steps + |> Enum.reduce({t_c_start, p_start}, fn _, {t_c, p} -> + t_k = t_c + 273.15 + gamma_s = saturated_lapse_rate_per_pressure(t_k, p) + t_next_k = t_k + gamma_s * dp + {t_next_k - 273.15, p + dp} + end) + |> elem(0) + end + + # γ_sat expressed as dT/dp (K/mb), derived from the standard moist + # adiabatic lapse rate dT/dz combined with hydrostatic dz/dp = + # -R_d × T_v / (g × p). + defp saturated_lapse_rate_per_pressure(t_k, p_mb) do + e_s = saturation_vapor_pressure_mb(t_k - 273.15) + r_s = @epsilon * e_s / max(p_mb - e_s, 1.0e-3) + + num = 1.0 + @l_v * r_s / (@r_d * t_k) + den = @c_p + @l_v * @l_v * r_s * @epsilon / (@r_d * t_k * t_k) + gamma_z = @g * num / den + + # dz/dp ≈ -R_d × T / (g × p) (mb) + dz_dp = -@r_d * t_k / (@g * p_mb) + -gamma_z * dz_dp + end + + defp saturation_vapor_pressure_mb(t_c) do + @buck_a * :math.exp(@buck_b * t_c / (@buck_c + t_c)) + end + + # ── CAPE / CIN ────────────────────────────────────────────────────── + + defp cape_cin(parcel_trace, sorted_desc) do + pairs = parcel_environment_pairs(parcel_trace, sorted_desc) + + {cape, cin} = + Enum.reduce(layer_increments(pairs), {0.0, 0.0}, fn {avg_diff, dlnp}, {cape, cin} -> + delta = @r_d * avg_diff * dlnp + + cond do + avg_diff > 0 -> {cape + delta, cin} + avg_diff < 0 -> {cape, cin + delta} + true -> {cape, cin} + end + end) + + {round_j(cape), round_j(cin)} + end + + defp cape_layer(_parcel, _profile, lo, hi) when not (is_number(lo) and is_number(hi)), do: nil + + defp cape_layer(parcel_trace, sorted_desc, lo_height_m, hi_height_m) do + pairs = + parcel_trace + |> parcel_environment_pairs(sorted_desc) + |> Enum.filter(fn p -> p.height >= lo_height_m and p.height <= hi_height_m end) + + cape = + pairs + |> layer_increments() + |> Enum.reduce(0.0, fn {avg_diff, dlnp}, acc -> + if avg_diff > 0, do: acc + @r_d * avg_diff * dlnp, else: acc + end) + + round_j(cape) + end + + # Build a list of {pressure, parcel_T_K, env_T_K, height} pairs at + # every environmental pressure level the parcel trace also covers. + defp parcel_environment_pairs(parcel_trace, sorted_desc) do + parcel_by_pres = Map.new(parcel_trace, &{&1.pres, &1.tmpc}) + + Enum.flat_map(sorted_desc, fn level -> + case Map.fetch(parcel_by_pres, level.pres) do + {:ok, parcel_t} -> + [ + %{ + pres: level.pres, + parcel_t_k: parcel_t + 273.15, + env_t_k: level.tmpc + 273.15, + diff: parcel_t - level.tmpc, + height: level.hght + } + ] + + :error -> + [] + end + end) + end + + defp layer_increments(pairs) do + pairs + |> Enum.sort_by(& &1.pres, :desc) + |> Enum.chunk_every(2, 1, :discard) + |> Enum.map(fn [a, b] -> + avg = (a.diff + b.diff) / 2.0 + dlnp = :math.log(a.pres / b.pres) + {avg, dlnp} + end) + end + + # ── LFC, EL ───────────────────────────────────────────────────────── + + defp lfc_el(parcel_trace, sorted_desc) do + pairs = + parcel_trace + |> parcel_environment_pairs(sorted_desc) + |> Enum.sort_by(& &1.pres, :desc) + + pos_pairs = Enum.filter(pairs, &(&1.diff > 0)) + + case pos_pairs do + [] -> + {nil, nil} + + [first | _] -> + last = List.last(pos_pairs) + # LFC = bottom of the first contiguous positive-area region. + # EL = top of the last positive-area region. Coarse but + # matches what SPC publishes for single-parcel ascent. + {first.pres, last.pres} + end + end + + # ── Lifted index ──────────────────────────────────────────────────── + + defp lifted_index(parcel_trace, sorted_desc) do + parcel_t = pressure_lookup_temp(parcel_trace, 500.0) + env_t = pressure_lookup_temp(sorted_desc, 500.0) + + case {parcel_t, env_t} do + {nil, _} -> nil + {_, nil} -> nil + {pt, et} -> Float.round(et - pt, 1) + end + end + + defp pressure_lookup_temp(levels, target_p) do + sorted = Enum.sort_by(levels, & &1.pres, :desc) + + cond do + sorted == [] -> nil + target_p > List.first(sorted).pres -> List.first(sorted).tmpc + target_p < List.last(sorted).pres -> List.last(sorted).tmpc + true -> interp_temp_at_pressure(sorted, target_p) + end + end + + defp interp_temp_at_pressure(sorted, target_p) do + sorted + |> Enum.chunk_every(2, 1, :discard) + |> Enum.find_value(fn [a, b] -> + if target_p <= a.pres and target_p >= b.pres do + f = (a.pres - target_p) / max(a.pres - b.pres, 1.0e-6) + a.tmpc + f * (b.tmpc - a.tmpc) + end + end) + end + + # ── Freezing / wet-bulb-zero / DCAPE ──────────────────────────────── + + defp temperature_height_crossing(sorted_desc, fun, target_c) do + sorted_desc + |> Enum.chunk_every(2, 1, :discard) + |> Enum.find_value(fn [a, b] -> + ta = fun.(a) + tb = fun.(b) + + cond do + ta == nil or tb == nil -> + nil + + (ta - target_c) * (tb - target_c) > 0 -> + nil + + ta == tb -> + a.hght + + true -> + f = (ta - target_c) / (ta - tb) + a.hght + f * (b.hght - a.hght) + end + end) + end + + defp wet_bulb_c(level) do + # Stipanuk-style iterative approximation: T_wb satisfies + # e_s(T_wb) = e + (p × c_p / (ε × L_v)) × (T - T_wb) + # Solve by bracketing between dewpoint and dry bulb. + t_d = level.dwpc + t = level.tmpc + + if t_d == nil or t == nil do + nil + else + bisect_wet_bulb(t, t_d, level.pres) + end + end + + defp bisect_wet_bulb(t, t_d, _p) when t == t_d, do: t + + defp bisect_wet_bulb(t, t_d, p) do + e = saturation_vapor_pressure_mb(t_d) + + f = fn t_w -> + e_s = saturation_vapor_pressure_mb(t_w) + e_s - e - p * @c_p / (@epsilon * @l_v) * (t - t_w) + end + + bisect(f, t_d, t, 30) + end + + defp bisect(_f, lo, hi, 0), do: (lo + hi) / 2.0 + + defp bisect(f, lo, hi, iter) do + mid = (lo + hi) / 2.0 + val_lo = f.(lo) + val_mid = f.(mid) + + cond do + abs(val_mid) < 1.0e-3 -> mid + val_lo * val_mid < 0 -> bisect(f, lo, mid, iter - 1) + true -> bisect(f, mid, hi, iter - 1) + end + end + + defp dcape(sorted_desc) do + # Pick the level with minimum θ_e in the layer 700-400 mb as the + # downdraft source. Saturate the parcel there, descend moist- + # adiabatically to the surface, and integrate negative buoyancy. + candidates = + Enum.filter(sorted_desc, fn l -> l.pres >= 400.0 and l.pres <= 700.0 end) + + case candidates do + [] -> + nil + + _ -> + source = Enum.min_by(candidates, &theta_e(&1.tmpc, &1.dwpc, &1.pres)) + descend_dcape(source, sorted_desc) + end + end + + defp descend_dcape(source, sorted_desc) do + descent_targets = + sorted_desc + |> Enum.filter(&(&1.pres >= source.pres)) + |> Enum.sort_by(& &1.pres, :asc) + + {traj, _} = + Enum.map_reduce(descent_targets, %{pres: source.pres, tmpc: source.dwpc}, fn lvl, prev -> + t_next = saturated_adiabatic_step(prev.tmpc, prev.pres, lvl.pres) + step = %{pres: lvl.pres, tmpc: t_next, env_t: lvl.tmpc, height: lvl.hght} + {step, %{pres: lvl.pres, tmpc: t_next}} + end) + + integral = + traj + |> Enum.chunk_every(2, 1, :discard) + |> Enum.reduce(0.0, fn [a, b], acc -> + diff = (a.env_t - a.tmpc + b.env_t - b.tmpc) / 2.0 + dlnp = :math.log(b.pres / a.pres) + + if diff > 0, do: acc + @r_d * diff * dlnp, else: acc + end) + + round_j(integral) + end + + defp theta_e(t_c, td_c, p_mb) do + t_k = t_c + 273.15 + e = saturation_vapor_pressure_mb(td_c) + r = @epsilon * e / max(p_mb - e, 1.0e-3) + theta = t_k * :math.pow(1000.0 / p_mb, @r_d / @c_p) + theta * :math.exp(@l_v * r / (@c_p * t_k)) + end + + # ── Lapse rates ───────────────────────────────────────────────────── + + defp lapse_rate_height(sorted_desc, lo_m, hi_m) do + case {temp_at_height(sorted_desc, lo_m), temp_at_height(sorted_desc, hi_m)} do + {t_lo, t_hi} when is_number(t_lo) and is_number(t_hi) -> + Float.round((t_lo - t_hi) / ((hi_m - lo_m) / 1000.0), 2) + + _ -> + nil + end + end + + defp lapse_rate_pressures(sorted_desc, p_lo, p_hi) do + t_lo = pressure_lookup_temp(sorted_desc, p_lo) + t_hi = pressure_lookup_temp(sorted_desc, p_hi) + h_lo = height_at_pressure(sorted_desc, p_lo) + h_hi = height_at_pressure(sorted_desc, p_hi) + + cond do + not (is_number(t_lo) and is_number(t_hi)) -> nil + not (is_number(h_lo) and is_number(h_hi)) -> nil + h_hi == h_lo -> nil + true -> Float.round((t_lo - t_hi) / ((h_hi - h_lo) / 1000.0), 2) + end + end + + defp temp_at_height(sorted_desc, target_m) do + sorted_desc + |> Enum.chunk_every(2, 1, :discard) + |> Enum.find_value(fn [a, b] -> + if a.hght <= target_m and target_m <= b.hght and b.hght != a.hght do + f = (target_m - a.hght) / (b.hght - a.hght) + a.tmpc + f * (b.tmpc - a.tmpc) + end + end) + end + + defp height_at_pressure(_sorted, nil), do: nil + + defp height_at_pressure(sorted_desc, target_p) do + sorted_desc + |> Enum.chunk_every(2, 1, :discard) + |> Enum.find_value(fn [a, b] -> + if target_p <= a.pres and target_p >= b.pres and a.pres != b.pres do + f = (a.pres - target_p) / (a.pres - b.pres) + a.hght + f * (b.hght - a.hght) + end + end) + end + + # ── Profile shape normaliser ──────────────────────────────────────── + + defp normalize_level(level) when is_map(level) do + pres = number(get_any(level, ["pres", "pres_mb", :pres, :pres_mb])) + hght = number(get_any(level, ["hght", "hght_m", :hght, :hght_m])) + tmpc = number(get_any(level, ["tmpc", :tmpc])) + dwpc = number(get_any(level, ["dwpc", :dwpc])) + + if is_number(pres) and is_number(hght) and is_number(tmpc) do + %{pres: pres, hght: hght, tmpc: tmpc, dwpc: dwpc} + end + end + + defp normalize_level(_), do: nil + + defp get_any(_map, []), do: nil + + defp get_any(map, [key | rest]) do + case Map.fetch(map, key) do + {:ok, v} -> v + :error -> get_any(map, rest) + end + end + + defp number(nil), do: nil + defp number(x) when is_integer(x), do: x * 1.0 + defp number(x) when is_float(x), do: x + defp number(_), do: nil + + defp round_j(nil), do: nil + defp round_j(v) when is_number(v), do: Float.round(v * 1.0, 1) +end diff --git a/lib/microwaveprop_web/live/skewt_live.ex b/lib/microwaveprop_web/live/skewt_live.ex index c32f3793..ac6939ba 100644 --- a/lib/microwaveprop_web/live/skewt_live.ex +++ b/lib/microwaveprop_web/live/skewt_live.ex @@ -13,6 +13,7 @@ defmodule MicrowavepropWeb.SkewtLive do alias Microwaveprop.Propagation.ProfilesFile alias Microwaveprop.Weather.HrrrProfileLookup + alias Microwaveprop.Weather.SkewtParams alias Microwaveprop.Weather.SoundingParams alias MicrowavepropWeb.SkewtLocationResolver alias MicrowavepropWeb.SkewtSvg @@ -120,8 +121,13 @@ defmodule MicrowavepropWeb.SkewtLive do if profile == [] do {nil, nil, nil} else - derived = SoundingParams.derive(profile) - svg = SkewtSvg.render(profile) + sounding = SoundingParams.derive(profile) + spc = SkewtParams.derive(profile) + # Merge so the LiveView template only has to read one map. + # SoundingParams keys (atoms) and SkewtParams keys (atoms) + # are disjoint by construction. + derived = Map.merge(sounding || %{}, spc) + svg = SkewtSvg.render(profile, parcel_trace: spc[:parcel_trace] || []) {profile, derived, svg} end end @@ -268,19 +274,28 @@ defmodule MicrowavepropWeb.SkewtLive do > {Phoenix.HTML.raw(@svg)} -
-

Derived parameters

+
+

Sounding parameters

<%= if @derived do %> -
- <%= for {label, value} <- derived_rows(@derived) do %> -
{label}
-
{value}
- <% end %> -
+ <%= for {section_title, rows} <- derived_sections(@derived) do %> +
+

+ {section_title} +

+
+ <%= for {label, value} <- rows do %> +
{label}
+
{value}
+ <% end %> +
+
+ <% end %> <%= if @derived[:ducting_detected] && @derived[:ducts] != [] do %> -
-

Ducts

+
+

+ Ducts +

<%= for duct <- @derived.ducts do %>
base {round(duct["base"])} m → top {round(duct["top"])} m, ΔM {duct["strength"]} @@ -288,6 +303,12 @@ defmodule MicrowavepropWeb.SkewtLive do <% end %>
<% end %> + +

+ Wind-derived indices (SRH, BWD, Bunkers motion, SCP/STP/SHIP) + aren't shown — HRRR persists wind only at 10 m AGL, so per- + level shear can't be computed from this data source. +

<% end %>
@@ -328,19 +349,52 @@ defmodule MicrowavepropWeb.SkewtLive do IO.iodata_to_binary("f#{:io_lib.format("~2..0B", [delta_h])} · #{Calendar.strftime(t, "%H:%MZ")}") end - defp derived_rows(d) do + # Grouped derived-parameter sections rendered on the right side of + # the page. Mirrors the SPC viewer's "parcel / level / lapse rate / + # moisture / refractivity" layout so users moving between the two + # tools don't have to context-switch. + defp derived_sections(d) do [ - {"Surface T", fmt_temp(d[:surface_temp_c])}, - {"Surface Td", fmt_temp(d[:surface_dewpoint_c])}, - {"Surface P", fmt_num(d[:surface_pressure_mb], "mb")}, - {"Surface N", fmt_num(d[:surface_refractivity], "")}, - {"Min dN/dh", fmt_num(d[:min_refractivity_gradient], "/km")}, - {"PWAT", fmt_num(d[:precipitable_water_mm], "mm")}, - {"BL depth", fmt_num(d[:boundary_layer_depth_m], "m")}, - {"K-index", fmt_num(d[:k_index], "")}, - {"Lifted index", fmt_num(d[:lifted_index], "")}, - {"Ducting", if(d[:ducting_detected], do: "yes", else: "no")}, - {"Levels", to_string(d[:level_count] || 0)} + {"Surface", + [ + {"T", fmt_temp(d[:surface_temp_c])}, + {"Td", fmt_temp(d[:surface_dewpoint_c])}, + {"P", fmt_num(d[:surface_pressure_mb], "mb")} + ]}, + {"Convective parcel", + [ + {"SBCAPE", fmt_num(d[:sbcape], "J/kg")}, + {"SBCIN", fmt_num(d[:sbcin], "J/kg")}, + {"3 km CAPE", fmt_num(d[:cape_3km], "J/kg")}, + {"Lifted index", fmt_num(d[:lifted_index], "")}, + {"DCAPE", fmt_num(d[:dcape], "J/kg")} + ]}, + {"Levels", + [ + {"LCL", fmt_pressure_height(d[:lcl_pressure_mb], d[:lcl_height_m])}, + {"LFC", fmt_pressure_height(d[:lfc_pressure_mb], d[:lfc_height_m])}, + {"EL", fmt_pressure_height(d[:el_pressure_mb], d[:el_height_m])}, + {"0 °C", fmt_height(d[:freezing_level_m])}, + {"Wet-bulb 0", fmt_height(d[:wbz_m])} + ]}, + {"Lapse rates", + [ + {"Sfc–3 km", fmt_num(d[:lapse_rate_sfc_3km_c_per_km], "°C/km")}, + {"700–500 mb", fmt_num(d[:lapse_rate_700_500_c_per_km], "°C/km")}, + {"850–500 mb", fmt_num(d[:lapse_rate_850_500_c_per_km], "°C/km")} + ]}, + {"Moisture / stability", + [ + {"PWAT", fmt_num(d[:precipitable_water_mm], "mm")}, + {"BL depth", fmt_num(d[:boundary_layer_depth_m], "m")}, + {"K-index", fmt_num(d[:k_index], "")} + ]}, + {"Refractivity", + [ + {"Surface N", fmt_num(d[:surface_refractivity], "")}, + {"Min dN/dh", fmt_num(d[:min_refractivity_gradient], "/km")}, + {"Ducting", if(d[:ducting_detected], do: "yes", else: "no")} + ]} ] end @@ -351,6 +405,19 @@ defmodule MicrowavepropWeb.SkewtLive do defp fmt_num(v, ""), do: format_value(v) defp fmt_num(v, suffix), do: "#{format_value(v)} #{suffix}" + defp fmt_height(nil), do: "—" + + defp fmt_height(v) when is_number(v) do + "#{round(v)} m / #{round(v * 3.28084)} ft" + end + + defp fmt_pressure_height(nil, _), do: "—" + defp fmt_pressure_height(p, nil), do: "#{round(p)} mb" + + defp fmt_pressure_height(p, h) when is_number(p) and is_number(h) do + "#{round(p)} mb · #{round(h)} m" + end + defp format_value(v) when is_float(v), do: v |> Float.round(1) |> :erlang.float_to_binary(decimals: 1) defp format_value(v) when is_integer(v), do: Integer.to_string(v) defp format_value(v), do: to_string(v) diff --git a/lib/microwaveprop_web/live/skewt_svg.ex b/lib/microwaveprop_web/live/skewt_svg.ex index ff13565a..42b0fc4c 100644 --- a/lib/microwaveprop_web/live/skewt_svg.ex +++ b/lib/microwaveprop_web/live/skewt_svg.ex @@ -32,10 +32,16 @@ defmodule MicrowavepropWeb.SkewtSvg do 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()]) :: iodata() - def render(profile) when is_list(profile) do + @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||, @@ -48,11 +54,37 @@ defmodule MicrowavepropWeb.SkewtSvg do frame(), axis_labels(), profile_traces(cleaned), + parcel_trace_line(parcel), cursor_layer(), "" ] 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|| + 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`. @@ -73,25 +105,40 @@ defmodule MicrowavepropWeb.SkewtSvg do """ 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 """ @@ -307,13 +354,19 @@ defmodule MicrowavepropWeb.SkewtSvg do -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 - # 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 + Enum.map(@mixing_ratios, fn w -> pts = sample_pressures() - |> Enum.filter(&(&1 >= 400)) + |> 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) @@ -324,7 +377,30 @@ defmodule MicrowavepropWeb.SkewtSvg do end) |> clip_polyline() - polyline_or_empty(pts, "grid-mix") + 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|#{label}| + else + "" end end diff --git a/test/microwaveprop/weather/skewt_params_test.exs b/test/microwaveprop/weather/skewt_params_test.exs new file mode 100644 index 00000000..25cf486d --- /dev/null +++ b/test/microwaveprop/weather/skewt_params_test.exs @@ -0,0 +1,118 @@ +defmodule Microwaveprop.Weather.SkewtParamsTest do + use ExUnit.Case, async: true + + alias Microwaveprop.Weather.SkewtParams + + # A moderately unstable summer-like profile. Surface ≈ 30 °C / 22 °C + # dewpoint, gradually drier and colder aloft. Used as a baseline for + # LCL/LFC/EL/CAPE assertions. + @summer_profile [ + %{"pres" => 1000.0, "hght" => 100.0, "tmpc" => 30.0, "dwpc" => 22.0}, + %{"pres" => 925.0, "hght" => 800.0, "tmpc" => 24.0, "dwpc" => 19.0}, + %{"pres" => 850.0, "hght" => 1500.0, "tmpc" => 18.0, "dwpc" => 14.0}, + %{"pres" => 700.0, "hght" => 3100.0, "tmpc" => 8.0, "dwpc" => -2.0}, + %{"pres" => 500.0, "hght" => 5800.0, "tmpc" => -10.0, "dwpc" => -25.0}, + %{"pres" => 300.0, "hght" => 9300.0, "tmpc" => -38.0, "dwpc" => -55.0}, + %{"pres" => 200.0, "hght" => 12_000.0, "tmpc" => -60.0, "dwpc" => -80.0} + ] + + # A capped / stable profile: subsidence inversion above the surface + # mixed layer. Should produce 0 CAPE. + @capped_profile [ + %{"pres" => 1000.0, "hght" => 100.0, "tmpc" => 25.0, "dwpc" => 18.0}, + %{"pres" => 925.0, "hght" => 800.0, "tmpc" => 28.0, "dwpc" => 5.0}, + %{"pres" => 850.0, "hght" => 1500.0, "tmpc" => 25.0, "dwpc" => -5.0}, + %{"pres" => 700.0, "hght" => 3100.0, "tmpc" => 12.0, "dwpc" => -15.0}, + %{"pres" => 500.0, "hght" => 5800.0, "tmpc" => -8.0, "dwpc" => -30.0} + ] + + describe "lcl/3" do + test "returns LCL pressure and temperature consistent with Bolton (1980)" do + # Surface 30°C / 22°C dewpoint at 1000 mb → LCL ~880-920 mb, + # T_LCL ~17-22°C. Tolerances stay loose because Bolton is an + # approximation and the surface RH already drives most of the + # answer. + {p_lcl, t_lcl} = SkewtParams.lcl(30.0, 22.0, 1000.0) + assert p_lcl >= 850.0 and p_lcl <= 940.0 + assert t_lcl >= 16.0 and t_lcl <= 23.0 + end + + test "LCL pressure equals surface pressure when air is saturated" do + {p_lcl, t_lcl} = SkewtParams.lcl(20.0, 20.0, 1000.0) + assert_in_delta p_lcl, 1000.0, 1.0 + assert_in_delta t_lcl, 20.0, 0.2 + end + end + + describe "derive/1" do + test "produces every field documented in the public spec for a non-trivial profile" do + params = SkewtParams.derive(@summer_profile) + + # Surface-based parcel parameters present. + assert is_map(params) + assert is_number(params[:lcl_pressure_mb]) + assert is_number(params[:lcl_height_m]) + assert is_number(params[:sbcape]) + assert is_number(params[:sbcin]) + assert is_number(params[:lifted_index]) + # LFC/EL may be nil if no positive area exists — but for this + # unstable summer profile they must resolve. + assert is_number(params[:lfc_pressure_mb]) + assert is_number(params[:el_pressure_mb]) + # Cold-side milestones. + assert is_number(params[:freezing_level_m]) + assert is_number(params[:wbz_m]) + # Layer lapse rates. + assert is_number(params[:lapse_rate_700_500_c_per_km]) + assert is_number(params[:lapse_rate_sfc_3km_c_per_km]) + # Parcel trajectory used by SkewtSvg to draw the dashed grey line. + assert is_list(params[:parcel_trace]) + assert length(params[:parcel_trace]) >= 2 + end + + test "summer profile is unstable: SBCAPE > 0, SBCIN <= 0, LFC < surface" do + params = SkewtParams.derive(@summer_profile) + + assert params[:sbcape] > 0 + # CIN by convention is negative (or zero). + assert params[:sbcin] <= 0 + assert params[:lfc_pressure_mb] < 1000.0 + assert params[:el_pressure_mb] < params[:lfc_pressure_mb] + # Lifted index < 0 for an unstable summer atmosphere. + assert params[:lifted_index] < 0 + end + + test "capped profile is stable: SBCAPE = 0, no LFC" do + params = SkewtParams.derive(@capped_profile) + assert params[:sbcape] == 0.0 + assert params[:lfc_pressure_mb] == nil + assert params[:el_pressure_mb] == nil + # Lifted index > 0 in a stable atmosphere. + assert params[:lifted_index] > 0 + end + + test "freezing level is the height where T crosses 0 °C" do + params = SkewtParams.derive(@summer_profile) + # Profile goes 30→24→18→8→-10°C between 100 m and 5800 m. + # 0°C falls between 700 mb (8°C, 3100 m) and 500 mb (-10°C, 5800 m). + # Linear interpolation puts it ~4300 m. + assert params[:freezing_level_m] >= 3700 + assert params[:freezing_level_m] <= 4900 + end + + test "lapse rate sfc-3km is the (T_sfc - T_3km) / 3 km gradient" do + params = SkewtParams.derive(@summer_profile) + # 30 °C at 100 m, ~9 °C at 3100 m → ~7 °C/km drop. + assert params[:lapse_rate_sfc_3km_c_per_km] >= 5.5 + assert params[:lapse_rate_sfc_3km_c_per_km] <= 9.0 + end + + test "empty profile returns a fully-nil parameter map (no crash)" do + params = SkewtParams.derive([]) + assert is_map(params) + assert params[:sbcape] == nil + assert params[:lcl_pressure_mb] == nil + assert params[:parcel_trace] == [] + end + end +end