feat(skewt): NOAA SPC line colors + sounding-parameter parity
Two changes that bring /skewt visually and analytically into parity with the NOAA SPC sounding viewer (https://www.spc.noaa.gov/exper/soundings/): NOAA-style line set ------------------- Updated the SkewtSvg style block to match the public SPC chart: * Temperature trace ........ pure red (#ff0000) at 2.5 px stroke * Dewpoint trace ........... pure green (#00aa00) at 2.5 px stroke * Isotherms ................ dashed teal (#14b8a6) * Dry adiabats ............. solid orange (#f59e0b) * Moist (saturated) adiabats dashed sky-blue (#38bdf8) * Mixing-ratio lines ....... dotted darker green (#16a34a), with labels at 700 mb showing the g/kg value, drawn only below 600 mb (matches the SPC chart's lower-band number row) * Parcel ascent ............ new dashed grey (#475569) line, drawn whenever SkewtParams.derive returns a non-empty parcel_trace The mixing-ratio set is the SPC standard (0.4 / 0.7 / 1 / 2 / 3 / 5 / 8 / 12 / 16 / 20 g/kg) instead of the prior arbitrary set. SkewtParams: SPC parcel-theory + indices module ----------------------------------------------- New `Microwaveprop.Weather.SkewtParams` derives the parameter set that ships in the SPC `*.txt` companion file: * LCL pressure / temperature / height (Bolton 1980 eq 22) * LFC pressure / height * EL pressure / height * SBCAPE, SBCIN (J/kg, virtual-temperature integration with surface parcel) * 3 km CAPE * Lifted Index (parcel ↔ environment T at 500 mb) * Freezing level (T = 0 °C height) and Wet-Bulb Zero * DCAPE (downdraft CAPE — picks the min-θₑ source level in 400–700 mb, descends moist-adiabatically to surface) * Layer lapse rates: sfc–3 km, 700–500 mb, 850–500 mb * `parcel_trace` — pressure/temperature trajectory used by SkewtSvg to draw the dashed parcel-ascent overlay Wind-derived indices (SRH, BWD, Bunkers motion, SCP, STP, SHIP) are deliberately *not* produced and the LiveView surfaces a one-line note explaining why: HRRR persists wind only at 10 m AGL, so per-pressure- level shear can't be computed from this data source. LiveView panel restructured into grouped sections (Surface, Convective parcel, Levels, Lapse rates, Moisture/stability, Refractivity) so the layout matches the SPC viewer's order. The Levels rows show pressure *and* height ("872 mb · 1,140 m") for quick cross-reference between the diagram axis and the readout. Tests: 8/8 new SkewtParamsTest cases (LCL Bolton check, saturated-air edge case, full field-set presence, summer profile produces SBCAPE > 0 and LI < 0, capped profile produces SBCAPE = 0 and LI > 0, freezing level interpolated correctly, sfc–3 km lapse rate, empty profile no crash). 25/25 across all skewt-related test files. mix test green (2,902 / 2,902 + 221 properties), mix credo --strict clean, mix assets.build clean.
This commit is contained in:
parent
a4cd7632eb
commit
be1174914f
4 changed files with 835 additions and 41 deletions
533
lib/microwaveprop/weather/skewt_params.ex
Normal file
533
lib/microwaveprop/weather/skewt_params.ex
Normal file
|
|
@ -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
|
||||
|
|
@ -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)}
|
||||
</div>
|
||||
<div class="bg-base-200 rounded-box p-4 text-sm space-y-3">
|
||||
<h3 class="font-semibold text-base">Derived parameters</h3>
|
||||
<div class="bg-base-200 rounded-box p-4 text-sm space-y-4">
|
||||
<h3 class="font-semibold text-base">Sounding parameters</h3>
|
||||
<%= if @derived do %>
|
||||
<dl class="grid grid-cols-2 gap-x-3 gap-y-1.5">
|
||||
<%= for {label, value} <- derived_rows(@derived) do %>
|
||||
<dt class="opacity-70">{label}</dt>
|
||||
<dd class="font-mono text-right">{value}</dd>
|
||||
<% end %>
|
||||
</dl>
|
||||
<%= for {section_title, rows} <- derived_sections(@derived) do %>
|
||||
<div class="space-y-1">
|
||||
<h4 class="font-semibold text-xs uppercase tracking-wide opacity-70">
|
||||
{section_title}
|
||||
</h4>
|
||||
<dl class="grid grid-cols-2 gap-x-3 gap-y-1">
|
||||
<%= for {label, value} <- rows do %>
|
||||
<dt class="opacity-70">{label}</dt>
|
||||
<dd class="font-mono text-right">{value}</dd>
|
||||
<% end %>
|
||||
</dl>
|
||||
</div>
|
||||
<% end %>
|
||||
|
||||
<%= if @derived[:ducting_detected] && @derived[:ducts] != [] do %>
|
||||
<div class="mt-3 pt-3 border-t border-base-300 space-y-1">
|
||||
<h4 class="font-semibold">Ducts</h4>
|
||||
<div class="pt-2 border-t border-base-300 space-y-1">
|
||||
<h4 class="font-semibold text-xs uppercase tracking-wide opacity-70">
|
||||
Ducts
|
||||
</h4>
|
||||
<%= for duct <- @derived.ducts do %>
|
||||
<div class="font-mono text-xs">
|
||||
base {round(duct["base"])} m → top {round(duct["top"])} m, ΔM {duct["strength"]}
|
||||
|
|
@ -288,6 +303,12 @@ defmodule MicrowavepropWeb.SkewtLive do
|
|||
<% end %>
|
||||
</div>
|
||||
<% end %>
|
||||
|
||||
<p class="pt-2 border-t border-base-300 text-xs opacity-60 leading-relaxed">
|
||||
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.
|
||||
</p>
|
||||
<% end %>
|
||||
</div>
|
||||
</div>
|
||||
|
|
@ -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)
|
||||
|
|
|
|||
|
|
@ -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|<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">|,
|
||||
|
|
@ -48,11 +54,37 @@ defmodule MicrowavepropWeb.SkewtSvg do
|
|||
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`.
|
||||
|
|
@ -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
|
||||
"""
|
||||
<defs>
|
||||
<style>
|
||||
.frame { fill: none; stroke: currentColor; stroke-width: 1.2; }
|
||||
.grid-major { fill: none; stroke: currentColor; stroke-opacity: 0.35; stroke-width: 0.6; }
|
||||
.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: #b91c1c; stroke-opacity: 0.30; stroke-width: 0.5; }
|
||||
.grid-dry { fill: none; stroke: #b45309; stroke-opacity: 0.30; stroke-width: 0.4; stroke-dasharray: 4 3; }
|
||||
.grid-moist { fill: none; stroke: #047857; stroke-opacity: 0.30; stroke-width: 0.4; stroke-dasharray: 1 4; }
|
||||
.grid-mix { fill: none; stroke: #6d28d9; stroke-opacity: 0.30; stroke-width: 0.4; stroke-dasharray: 2 2; }
|
||||
.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: #b91c1c; fill-opacity: 0.7; }
|
||||
.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: #dc2626; stroke-width: 2; stroke-linejoin: round; stroke-linecap: round; }
|
||||
.dew-line { fill: none; stroke: #16a34a; stroke-width: 2; stroke-linejoin: round; stroke-linecap: round; }
|
||||
.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: #dc2626; pointer-events: none; }
|
||||
.skewt-cursor-dot-d { fill: #16a34a; 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>
|
||||
|
|
@ -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|<text class="mix-label" x="#{r(x)}" y="#{r(y - 2)}" text-anchor="middle">#{label}</text>|
|
||||
else
|
||||
""
|
||||
end
|
||||
end
|
||||
|
||||
|
|
|
|||
118
test/microwaveprop/weather/skewt_params_test.exs
Normal file
118
test/microwaveprop/weather/skewt_params_test.exs
Normal file
|
|
@ -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
|
||||
Loading…
Add table
Reference in a new issue