fix: April 2026 codebase review — address 13 bugs across propagation chain

Each fix is covered by a regression test that fails on `main` and
passes on this commit.

Round 1 (initial review):

* propagation: thread `latitude` into the conditions map so
  `score_season/4` actually picks up regional multipliers
* hrrr_client / fetcher.rs: `nearest_hrrr_hour` rounds DOWN, never at
  a future cycle that NOAA hasn't published yet
* radio: spherical-vector great-circle midpoint replaces the
  arithmetic mean — anti-meridian paths no longer fold to Greenwich
* weather: `reconcile_weather_statuses` scales the longitude band by
  `1 / cos(lat)` so the bbox stays ~150 km wide at every latitude
* radio/maidenhead: clamp 90°/180° below the field-bucket overflow so
  `from_latlon` never emits invalid characters like 'S'
* prop_grid_rs/pipeline: merge HRRR + NEXRAD-derived rain rates and
  read `best_duct_freq_ghz` into `best_duct_band_ghz` so the Native
  Duct Boost actually fires
* propagation/region (Elixir + Rust): inclusive upper bounds so points
  exactly at lat_max get the regional multiplier
* weather/sounding_params (Elixir + Rust): drop the 10 m gradient
  floor so HRRR's thin near-surface layers stop hiding sharp ducts
* weather/sounding_params: when the profile ends inside a duct,
  finalize it with the highest sample as the top instead of throwing
  it away (Rust port already correct)

Round 2 (post-fix sweep):

* radio + commercial: single canonical haversine in Radio (atan2
  form); Commercial delegates instead of carrying a second copy that
  could disagree at threshold distances
* prop_grid_rs/profiles_file: `snap_coords` matches Elixir's
  step-aware snap (`round(coord/0.125) * 0.125`, then 3-dp round) so
  Rust-keyed and Elixir-keyed profile maps land on the same cell
* weather/grib2/wgrib2: `parse_lon_val_segment` uses `Float.parse`
  uniformly — wgrib2 dropping the trailing `.0` from a longitude no
  longer crashes the whole chain step
This commit is contained in:
Graham McIntire 2026-04-25 10:52:42 -05:00
parent 15050bd2d8
commit 6aa91e7656
No known key found for this signature in database
GPG key ID: F4ABF488E6029E59
24 changed files with 430 additions and 88 deletions

27
bugs.md Normal file
View file

@ -0,0 +1,27 @@
# Identified Bugs and Observations (Post-Fix Review)
This document reflects the state of the codebase after the fixes for the initial 15 bugs were applied.
All three items below were verified, covered by failing tests, and fixed on 2026-04-25.
## 1. Haversine Implementation Inconsistency — FIXED 2026-04-25
**Files:** `lib/microwaveprop/radio.ex` vs `lib/microwaveprop/commercial.ex`
**Description:**
- `Radio.haversine_km/4` uses the `asin` form: `2 * r * :math.asin(:math.sqrt(a))`
- `Commercial.haversine_km/4` uses the `atan2` form: `2 * r * :math.atan2(:math.sqrt(a), :math.sqrt(1 - a))`
**Impact:** While functionally similar, `atan2` is the industry standard for floating-point stability at large distances. Having two different implementations for the same physical calculation can lead to tiny discrepancies (sub-millimeter) that might cause a contact or commercial link to intermittently cross a distance threshold (e.g., 75 km) depending on which module performs the check.
**Fix:** Switched `Radio.haversine_km/4` to the `atan2` form and made `Commercial`'s private wrapper delegate to it. One canonical implementation, no per-module drift. Numerical-stability regression test added in `test/microwaveprop/radio_test.exs`.
## 2. Grid Snapping Mathematical Drift — FIXED 2026-04-25
**Files:** `lib/microwaveprop/propagation/profiles_file.ex` vs `rust/prop_grid_rs/src/profiles_file.rs`
**Description:**
- Elixir `snap/2`: `Float.round(Float.round(lat / step) * step, 3)`
- Rust `snap_coords`: `(lat * 1000.0).round() / 1000.0`
**Impact:** The Elixir version performs a division and multiplication by the grid step (`0.125`) before rounding. The Rust version rounds the raw coordinate to the third decimal place. For coordinates that land exactly on a $0.0625$ boundary (the halfway point between cells), these two methods may snap to different cells, causing the UI to report "No Data" or show the wrong factor breakdown for a clicked point.
**Fix:** Rust `snap_coords` now mirrors the Elixir step-aware snap (`round(coord / 0.125) * 0.125`, then 3-decimal round). Cross-stack parity test in `rust/prop_grid_rs/src/profiles_file.rs` enumerates known-disagreeing inputs.
## 3. Potential Exception in Wgrib2 Parsing (Theoretical) — FIXED 2026-04-25
**File:** `lib/microwaveprop/weather/grib2/wgrib2.ex`
**Description:** `parse_lon_val_segment/1` uses `String.to_float(lon_str)` for the longitude, but `Float.parse` for latitude and value.
**Impact:** If `wgrib2` were to ever return an integer string for longitude (e.g. `"235"`) instead of a float (`"235.0"`), `String.to_float` would raise an `ArgumentError`. While `wgrib2` typically outputs floats, using `Float.parse` would be more consistent with the other fields in the same line.
**Fix:** All three numeric fields now go through `Float.parse/1` (in a `with` chain so a malformed segment returns `nil` instead of bringing down the chain step). Function exposed as `@doc false` and covered by direct unit tests for both decimal-formatted and integer-formatted longitude strings.

View file

@ -166,20 +166,11 @@ defmodule Microwaveprop.Commercial do
defp link_endpoint(_), do: nil defp link_endpoint(_), do: nil
defp haversine_km(lat1, lon1, lat2, lon2) do # Single canonical haversine implementation lives in `Microwaveprop.Radio`.
r = 6371.0 # Two formulae for the same physical distance let identical coordinates
phi1 = lat1 * :math.pi() / 180 # land on opposite sides of a 75 km link-radius threshold depending on
phi2 = lat2 * :math.pi() / 180 # which module asked.
dphi = (lat2 - lat1) * :math.pi() / 180 defp haversine_km(lat1, lon1, lat2, lon2), do: Microwaveprop.Radio.haversine_km(lat1, lon1, lat2, lon2)
dlambda = (lon2 - lon1) * :math.pi() / 180
a =
:math.sin(dphi / 2) * :math.sin(dphi / 2) +
:math.cos(phi1) * :math.cos(phi2) * :math.sin(dlambda / 2) * :math.sin(dlambda / 2)
c = 2 * :math.atan2(:math.sqrt(a), :math.sqrt(1 - a))
r * c
end
defp average([]), do: 0.0 defp average([]), do: 0.0
defp average(list), do: Enum.sum(list) / length(list) defp average(list), do: Enum.sum(list) / length(list)

View file

@ -76,7 +76,7 @@ defmodule Microwaveprop.Propagation do
score_with_algorithm(hrrr_profile, valid_time, temp_c, dewpoint_c, derived, latitude, longitude) score_with_algorithm(hrrr_profile, valid_time, temp_c, dewpoint_c, derived, latitude, longitude)
end end
defp score_with_algorithm(hrrr_profile, valid_time, temp_c, dewpoint_c, derived, _latitude, longitude) do defp score_with_algorithm(hrrr_profile, valid_time, temp_c, dewpoint_c, derived, latitude, longitude) do
temp_f = Scorer.c_to_f(temp_c) temp_f = Scorer.c_to_f(temp_c)
dewpoint_f = Scorer.c_to_f(dewpoint_c) dewpoint_f = Scorer.c_to_f(dewpoint_c)
@ -89,6 +89,7 @@ defmodule Microwaveprop.Propagation do
utc_hour: valid_time.hour, utc_hour: valid_time.hour,
utc_minute: valid_time.minute, utc_minute: valid_time.minute,
month: valid_time.month, month: valid_time.month,
latitude: latitude,
longitude: longitude, longitude: longitude,
pressure_mb: hrrr_profile.surface_pressure_mb, pressure_mb: hrrr_profile.surface_pressure_mb,
prev_pressure_mb: nil, prev_pressure_mb: nil,

View file

@ -42,7 +42,7 @@ defmodule Microwaveprop.Propagation.Region do
@spec for_point(float, float) :: atom @spec for_point(float, float) :: atom
def for_point(lat, lon) do def for_point(lat, lon) do
Enum.find_value(@regions, :other, fn {name, {lat_min, lat_max}, {lon_min, lon_max}} -> Enum.find_value(@regions, :other, fn {name, {lat_min, lat_max}, {lon_min, lon_max}} ->
if lat >= lat_min and lat < lat_max and lon >= lon_min and lon < lon_max do if lat >= lat_min and lat <= lat_max and lon >= lon_min and lon <= lon_max do
name name
end end
end) end)

View file

@ -368,14 +368,43 @@ defmodule Microwaveprop.Radio do
defp extract_latlon(_), do: nil defp extract_latlon(_), do: nil
defp build_path_points({lat1, lon1}, {lat2, lon2}) do defp build_path_points({lat1, lon1}, {lat2, lon2}) do
mid_lat = (lat1 + lat2) / 2 {mid_lat, mid_lon} = great_circle_midpoint(lat1, lon1, lat2, lon2)
mid_lon = (lon1 + lon2) / 2
[{lat1, lon1}, {mid_lat, mid_lon}, {lat2, lon2}] [{lat1, lon1}, {mid_lat, mid_lon}, {lat2, lon2}]
end end
defp build_path_points({lat1, lon1}, _), do: [{lat1, lon1}] defp build_path_points({lat1, lon1}, _), do: [{lat1, lon1}]
defp build_path_points(_, _), do: [] defp build_path_points(_, _), do: []
# Spherical-vector midpoint between two lat/lon points. The arithmetic
# mean of two longitudes folds across the anti-meridian (e.g.
# `(179 + -179) / 2 == 0` back at Greenwich), so paths that cross the
# date line otherwise get a bogus center.
defp great_circle_midpoint(lat1, lon1, lat2, lon2) do
rlat1 = deg_to_rad(lat1)
rlat2 = deg_to_rad(lat2)
rlon1 = deg_to_rad(lon1)
dlon = deg_to_rad(lon2 - lon1)
bx = :math.cos(rlat2) * :math.cos(dlon)
by = :math.cos(rlat2) * :math.sin(dlon)
mid_lat =
:math.atan2(
:math.sin(rlat1) + :math.sin(rlat2),
:math.sqrt((:math.cos(rlat1) + bx) ** 2 + by ** 2)
)
mid_lon = rlon1 + :math.atan2(by, :math.cos(rlat1) + bx)
{rad_to_deg(mid_lat), normalize_lon(rad_to_deg(mid_lon))}
end
defp rad_to_deg(rad), do: rad * 180.0 / :math.pi()
defp normalize_lon(lon) when lon > 180.0, do: lon - 360.0
defp normalize_lon(lon) when lon < -180.0, do: lon + 360.0
defp normalize_lon(lon), do: lon
@earth_radius_km 6371.0 @earth_radius_km 6371.0
@spec haversine_km(number(), number(), number(), number()) :: float() @spec haversine_km(number(), number(), number(), number()) :: float()
@ -389,7 +418,10 @@ defmodule Microwaveprop.Radio do
:math.sin(dlat / 2) ** 2 + :math.sin(dlat / 2) ** 2 +
:math.cos(rlat1) * :math.cos(rlat2) * :math.sin(dlon / 2) ** 2 :math.cos(rlat1) * :math.cos(rlat2) * :math.sin(dlon / 2) ** 2
2 * @earth_radius_km * :math.asin(:math.sqrt(a)) # `atan2(√a, √(1a))` is the numerically stable form: when `a`
# rounds to slightly above 1 the asin form returns NaN, while
# atan2 stays finite.
2 * @earth_radius_km * :math.atan2(:math.sqrt(a), :math.sqrt(1 - a))
end end
@spec backfill_distances([Contact.t()]) :: Postgrex.Result.t() | nil @spec backfill_distances([Contact.t()]) :: Postgrex.Result.t() | nil

View file

@ -53,8 +53,11 @@ defmodule Microwaveprop.Radio.Maidenhead do
""" """
@spec from_latlon(number(), number(), pos_integer()) :: String.t() @spec from_latlon(number(), number(), pos_integer()) :: String.t()
def from_latlon(lat, lon, precision \\ 6) when is_number(lat) and is_number(lon) and precision in [4, 6, 8] do def from_latlon(lat, lon, precision \\ 6) when is_number(lat) and is_number(lon) and precision in [4, 6, 8] do
lat = lat + 90.0 # Clamp into the open ranges the field encoder expects. Exactly
lon = lon + 180.0 # `90.0` / `180.0` would otherwise overflow the trunc into 'S'
# (one past 'R'), producing grids outside the A..R alphabet.
lat = min(max(lat + 90.0, 0.0), 180.0 - 1.0e-9)
lon = min(max(lon + 180.0, 0.0), 360.0 - 1.0e-9)
# Field (A-R), 20° lon / 10° lat # Field (A-R), 20° lon / 10° lat
f1 = trunc(lon / 20.0) f1 = trunc(lon / 20.0)

View file

@ -366,9 +366,12 @@ defmodule Microwaveprop.Weather do
even after their data had landed. This reconciler closes that loop even after their data had landed. This reconciler closes that loop
as a single SQL UPDATE, invoked from the hourly enqueuer cron. as a single SQL UPDATE, invoked from the hourly enqueuer cron.
Radius encoded as a conservative ±1.5° lat/lon bounding box the Radius encoded as a ±1.5° latitude band; the longitude band is
same rectangle `weather_for_contact/2` uses at `radius_km: 150` scaled by `1 / cos(lat)` so the box covers the same physical
around the equator/mid-latitudes. east-west distance (~150 km) at every latitude. A fixed 1.5° lon
box collapses to ~110 km at lat 49° and would silently skip
observations the per-contact `weather_for_contact/2` query would
match.
""" """
@spec reconcile_weather_statuses() :: {:ok, non_neg_integer()} @spec reconcile_weather_statuses() :: {:ok, non_neg_integer()}
def reconcile_weather_statuses do def reconcile_weather_statuses do
@ -386,8 +389,8 @@ defmodule Microwaveprop.Weather do
AND o.observed_at <= c.qso_timestamp + interval '2 hours' AND o.observed_at <= c.qso_timestamp + interval '2 hours'
AND s.lat BETWEEN ((c.pos1->>'lat')::float - 1.5) AND s.lat BETWEEN ((c.pos1->>'lat')::float - 1.5)
AND ((c.pos1->>'lat')::float + 1.5) AND ((c.pos1->>'lat')::float + 1.5)
AND s.lon BETWEEN ((c.pos1->>'lon')::float - 1.5) AND s.lon BETWEEN ((c.pos1->>'lon')::float - 1.5 / GREATEST(cos(radians((c.pos1->>'lat')::float)), 0.01))
AND ((c.pos1->>'lon')::float + 1.5) AND ((c.pos1->>'lon')::float + 1.5 / GREATEST(cos(radians((c.pos1->>'lat')::float)), 0.01))
) )
OR ( OR (
c.pos2 IS NOT NULL c.pos2 IS NOT NULL
@ -399,8 +402,8 @@ defmodule Microwaveprop.Weather do
AND o.observed_at <= c.qso_timestamp + interval '2 hours' AND o.observed_at <= c.qso_timestamp + interval '2 hours'
AND s.lat BETWEEN ((c.pos2->>'lat')::float - 1.5) AND s.lat BETWEEN ((c.pos2->>'lat')::float - 1.5)
AND ((c.pos2->>'lat')::float + 1.5) AND ((c.pos2->>'lat')::float + 1.5)
AND s.lon BETWEEN ((c.pos2->>'lon')::float - 1.5) AND s.lon BETWEEN ((c.pos2->>'lon')::float - 1.5 / GREATEST(cos(radians((c.pos2->>'lat')::float)), 0.01))
AND ((c.pos2->>'lon')::float + 1.5) AND ((c.pos2->>'lon')::float + 1.5 / GREATEST(cos(radians((c.pos2->>'lat')::float)), 0.01))
) )
) )
) )

View file

@ -652,14 +652,24 @@ defmodule Microwaveprop.Weather.Grib2.Wgrib2 do
defp extract_var_level(_), do: nil defp extract_var_level(_), do: nil
# Parse "lon=242.958,lat=32.938,val=306.5" @doc false
defp parse_lon_val_segment(segment) do # Parse "lon=242.958,lat=32.938,val=306.5". Exposed (with @doc false)
# so the integer-lon edge case (wgrib2 occasionally drops the trailing
# `.0`) can be asserted directly.
def parse_lon_val_segment(segment) do
case Regex.run(~r/lon=([\d.]+),lat=([\d.]+),val=([\d.eE+-]+)/, segment) do case Regex.run(~r/lon=([\d.]+),lat=([\d.]+),val=([\d.eE+-]+)/, segment) do
[_, lon_str, lat_str, val_str] -> [_, lon_str, lat_str, val_str] ->
lon = denormalize_lon(String.to_float(lon_str)) # `Float.parse/1` accepts both `"243"` and `"243.0"`, where
{lat, ""} = Float.parse(lat_str) # `String.to_float/1` only accepts the latter. wgrib2 normally
{val, ""} = Float.parse(val_str) # emits the trailing `.0` but a build or locale flip that drops
{Float.round(lat, 3), Float.round(lon, 3), val} # it would crash the entire chain step.
with {lon_val, ""} <- Float.parse(lon_str),
{lat, ""} <- Float.parse(lat_str),
{val, ""} <- Float.parse(val_str) do
{Float.round(lat, 3), Float.round(denormalize_lon(lon_val), 3), val}
else
_ -> nil
end
_ -> _ ->
nil nil

View file

@ -215,16 +215,18 @@ defmodule Microwaveprop.Weather.HrrrClient do
end end
end end
@doc """
Truncates a DateTime down to the nearest published HRRR cycle hour.
Always rounds DOWN rounding to the nearest hour used to point at
the next cycle once we passed `:30`, which produced 404/403 errors
when that cycle had not yet been published by NOAA. Callers that need
a later valid_time should reach it via `forecast_hour` instead.
"""
@spec nearest_hrrr_hour(DateTime.t()) :: DateTime.t() @spec nearest_hrrr_hour(DateTime.t()) :: DateTime.t()
def nearest_hrrr_hour(dt) do def nearest_hrrr_hour(dt) do
total_seconds = dt.minute * 60 + dt.second total_seconds = dt.minute * 60 + dt.second
rounded_dt = DateTime.add(dt, -total_seconds, :second) DateTime.add(dt, -total_seconds, :second)
if dt.minute >= 30 do
DateTime.add(rounded_dt, 3600, :second)
else
rounded_dt
end
end end
@spec hrrr_url(Date.t(), non_neg_integer(), :surface | :pressure, non_neg_integer()) :: String.t() @spec hrrr_url(Date.t(), non_neg_integer(), :surface | :pressure, non_neg_integer()) :: String.t()

View file

@ -140,7 +140,10 @@ defmodule Microwaveprop.Weather.SoundingParams do
|> Enum.flat_map(fn [prev, curr] -> |> Enum.flat_map(fn [prev, curr] ->
dh_km = (curr.h_agl - prev.h_agl) / 1000.0 dh_km = (curr.h_agl - prev.h_agl) / 1000.0
if abs(dh_km) < 0.01 do # Only zero-thickness layers are skipped — the previous 10m guard
# discarded the very thin near-surface levels that HRRR's native
# vertical grid uses to resolve sharp inversions and surface ducts.
if abs(dh_km) < 1.0e-6 do
[] []
else else
dn = curr.n - prev.n dn = curr.n - prev.n
@ -175,14 +178,17 @@ defmodule Microwaveprop.Weather.SoundingParams do
classify_duct_transition(prev, curr, ducts, duct_state) classify_duct_transition(prev, curr, ducts, duct_state)
end) end)
# Handle case where profile ends while still in a duct # Surface-based ducts that extend past the top of a shallow
# profile still strongly affect ground-to-ground propagation, so
# close the duct using the highest sampled level as the top
# rather than discarding it for "lack of a clear top".
ducts = ducts =
case in_duct_state do case {in_duct_state, List.last(refract_profile)} do
%{base: _base, base_m: _base_m} -> {%{base: base, base_m: base_m}, %{} = top} ->
# Duct didn't close — discard (no clear top) {finalized, _} = finalize_duct(ducts, base, base_m, top)
ducts finalized
nil -> _ ->
ducts ducts
end end

View file

@ -183,14 +183,15 @@ pub fn merge_ranges(mut ranges: Vec<ByteRange>) -> Vec<ByteRange> {
out out
} }
/// Truncates a DateTime down to the nearest published HRRR cycle hour.
///
/// Always rounds DOWN — rounding up past `:30` used to point at the next
/// cycle, which produced 404/403 errors when that cycle had not yet been
/// published by NOAA. Callers that need a later valid_time should reach
/// it via `forecast_hour` instead.
pub fn nearest_hrrr_hour(dt: DateTime<Utc>) -> DateTime<Utc> { pub fn nearest_hrrr_hour(dt: DateTime<Utc>) -> DateTime<Utc> {
let total_seconds = dt.minute() as i64 * 60 + dt.second() as i64; let total_seconds = dt.minute() as i64 * 60 + dt.second() as i64;
let rounded = dt - chrono::Duration::seconds(total_seconds); dt - chrono::Duration::seconds(total_seconds)
if dt.minute() >= 30 {
rounded + chrono::Duration::hours(1)
} else {
rounded
}
} }
pub fn pressure_messages_grid() -> Vec<(String, String)> { pub fn pressure_messages_grid() -> Vec<(String, String)> {
@ -722,11 +723,15 @@ mod tests {
} }
#[test] #[test]
fn nearest_hour_rounds_properly() { fn nearest_hour_always_rounds_down() {
let base = Utc.with_ymd_and_hms(2026, 4, 19, 15, 29, 0).unwrap(); // Round-down keeps us from probing a cycle that NOAA has not
assert_eq!(nearest_hrrr_hour(base).hour(), 15); // published yet. Callers reach later valid_times via forecast_hour.
let up = Utc.with_ymd_and_hms(2026, 4, 19, 15, 30, 0).unwrap(); let early = Utc.with_ymd_and_hms(2026, 4, 19, 15, 29, 0).unwrap();
assert_eq!(nearest_hrrr_hour(up).hour(), 16); assert_eq!(nearest_hrrr_hour(early).hour(), 15);
let mid = Utc.with_ymd_and_hms(2026, 4, 19, 15, 30, 0).unwrap();
assert_eq!(nearest_hrrr_hour(mid).hour(), 15);
let late = Utc.with_ymd_and_hms(2026, 4, 19, 15, 59, 0).unwrap();
assert_eq!(nearest_hrrr_hour(late).hour(), 15);
} }
#[test] #[test]

View file

@ -247,12 +247,28 @@ fn cell_to_conditions(
.get("PRES:surface") .get("PRES:surface")
.copied() .copied()
.map(|pa| pa as f64 / 100.0); .map(|pa| pa as f64 / 100.0);
let rain_rate_mmhr = cell // Mirror Elixir's merged_rain_rate: pick the heavier of HRRR's
// accumulation-derived rate and NEXRAD's reflectivity-derived rate
// so a fast convective cell that hasn't yet shown up in the hourly
// APCP still triggers the rain penalty.
let hrrr_rate = cell
.get("APCP:surface") .get("APCP:surface")
.copied() .copied()
.map(|mm| mm as f64) .map(|mm| mm as f64)
.filter(|v| *v > 0.0); .unwrap_or(0.0);
let nexrad_rate = scorer::dbz_to_rain_rate_mmhr(
cell.get("nexrad_max_reflectivity_dbz")
.copied()
.map(|v| v as f64),
);
let merged_rate = hrrr_rate.max(nexrad_rate);
let rain_rate_mmhr = if merged_rate > 0.0 {
Some(merged_rate)
} else {
None
};
let bl_depth_m = cell.get("HPBL:surface").copied().map(|v| v as f64); let bl_depth_m = cell.get("HPBL:surface").copied().map(|v| v as f64);
let best_duct_band_ghz = cell.get("best_duct_freq_ghz").copied().map(|v| v as f64);
let levels: Vec<Level> = fetcher::grid_level_keys() let levels: Vec<Level> = fetcher::grid_level_keys()
.iter() .iter()
@ -287,7 +303,7 @@ fn cell_to_conditions(
min_refractivity_gradient, min_refractivity_gradient,
bl_depth_m, bl_depth_m,
pwat_mm, pwat_mm,
best_duct_band_ghz: None, best_duct_band_ghz,
bulk_richardson: None, bulk_richardson: None,
}) })
} }
@ -343,6 +359,43 @@ mod tests {
assert!(cell_to_conditions(&cell, 32.0, -97.0, &vt).is_none()); assert!(cell_to_conditions(&cell, 32.0, -97.0, &vt).is_none());
} }
#[test]
fn nexrad_reflectivity_lifts_rain_rate_above_hrrr_precip() {
// Mirrors the Elixir scorer's rain merge: when HRRR's hourly
// accumulation hasn't caught a fast convective cell but NEXRAD
// sees 45 dBZ overhead, take whichever rate is higher. Without
// this the Rust pipeline silently underestimates rain attenuation.
use chrono::TimeZone;
let mut cell = CellValues::new();
cell.insert("TMP:2 m above ground".into(), 295.0);
cell.insert("DPT:2 m above ground".into(), 285.0);
cell.insert("APCP:surface".into(), 0.0);
cell.insert("nexrad_max_reflectivity_dbz".into(), 45.0);
let vt = Utc.with_ymd_and_hms(2026, 6, 15, 18, 0, 0).unwrap();
let c = cell_to_conditions(&cell, 32.0, -97.0, &vt).unwrap();
let rate = c.rain_rate_mmhr.expect("nexrad-derived rain rate present");
assert!(rate > 5.0, "expected >5 mm/hr from 45 dBZ, got {rate}");
}
#[test]
fn duct_freq_threads_into_best_duct_band_ghz() {
// After native_duct::merge_duct_grid runs, every cell with a
// detected trapping layer carries `best_duct_freq_ghz`. The
// scorer reads `best_duct_band_ghz` from Conditions for the 15%
// Native Duct Boost. cell_to_conditions must wire one to the
// other or the boost never fires from the Rust pipeline.
use chrono::TimeZone;
let mut cell = CellValues::new();
cell.insert("TMP:2 m above ground".into(), 295.0);
cell.insert("DPT:2 m above ground".into(), 285.0);
cell.insert("best_duct_freq_ghz".into(), 24.0);
let vt = Utc.with_ymd_and_hms(2026, 6, 15, 18, 0, 0).unwrap();
let c = cell_to_conditions(&cell, 32.0, -97.0, &vt).unwrap();
assert_eq!(c.best_duct_band_ghz, Some(24.0));
}
#[test] #[test]
fn rejects_forecast_hour_zero() { fn rejects_forecast_hour_zero() {
use chrono::TimeZone; use chrono::TimeZone;

View file

@ -136,13 +136,26 @@ pub fn write_atomic(
Ok(path) Ok(path)
} }
/// Round a lat/lon to 3 decimal places — matches Elixir's /// Snap a lat/lon to the propagation grid step. Mirrors Elixir's
/// `Float.round(_, 3)` snap before keying into the profiles map. /// `Microwaveprop.Propagation.ProfilesFile.snap/2`:
///
/// ```elixir
/// Float.round(Float.round(coord / step) * step, 3)
/// ```
///
/// Two implementations of "snap" — Rust's old plain 3-decimal round
/// and Elixir's step-aware round — disagreed on every off-grid input,
/// so a coordinate written by Elixir was unreachable from Rust and
/// vice versa.
const GRID_STEP: f64 = 0.125;
pub fn snap_coords(lat: f64, lon: f64) -> (f64, f64) { pub fn snap_coords(lat: f64, lon: f64) -> (f64, f64) {
( (snap_step(lat), snap_step(lon))
(lat * 1000.0).round() / 1000.0, }
(lon * 1000.0).round() / 1000.0,
) fn snap_step(coord: f64) -> f64 {
let snapped = (coord / GRID_STEP).round() * GRID_STEP;
(snapped * 1000.0).round() / 1000.0
} }
/// Convenience helper: build a `rmpv::Value` map from an /// Convenience helper: build a `rmpv::Value` map from an
@ -239,7 +252,14 @@ mod tests {
} }
#[test] #[test]
fn snap_coords_rounds_to_three_decimals() { fn snap_coords_matches_elixir_step_snap() {
assert_eq!(snap_coords(32.8998, -97.04031), (32.900, -97.040)); // The Elixir reader writes the profiles map with keys snapped via
// `Float.round(Float.round(coord / 0.125) * 0.125, 3)`. This Rust
// helper has to produce the same keys or `read_point` looks up a
// coordinate that was never written and the UI shows "no data".
assert_eq!(snap_coords(32.8998, -97.04031), (32.875, -97.0));
assert_eq!(snap_coords(32.07, -97.07), (32.125, -97.125));
// Round-trip — already-snapped coordinates stay put.
assert_eq!(snap_coords(32.875, -97.0), (32.875, -97.0));
} }
} }

View file

@ -14,7 +14,9 @@ pub enum Region {
Other, Other,
} }
/// Lat min (inclusive), lat max (exclusive), lon min (inclusive), lon max (exclusive). /// Lat min, lat max, lon min, lon max — all inclusive on both ends so
/// that points sitting exactly on a region's upper boundary still pick
/// up the regional multiplier instead of falling through to `Other`.
struct BBox(Region, f64, f64, f64, f64); struct BBox(Region, f64, f64, f64, f64);
// Ordered to match Elixir's `Enum.find_value` short-circuit. // Ordered to match Elixir's `Enum.find_value` short-circuit.
@ -32,7 +34,7 @@ const REGIONS: &[BBox] = &[
pub fn for_point(lat: f64, lon: f64) -> Region { pub fn for_point(lat: f64, lon: f64) -> Region {
for bbox in REGIONS { for bbox in REGIONS {
let BBox(name, lat_min, lat_max, lon_min, lon_max) = *bbox; let BBox(name, lat_min, lat_max, lon_min, lon_max) = *bbox;
if lat >= lat_min && lat < lat_max && lon >= lon_min && lon < lon_max { if lat >= lat_min && lat <= lat_max && lon >= lon_min && lon <= lon_max {
return name; return name;
} }
} }
@ -102,4 +104,13 @@ mod tests {
// southeast (30..37, -90..-75). Elixir hits gulf_coast first. // southeast (30..37, -90..-75). Elixir hits gulf_coast first.
assert_eq!(for_point(30.5, -90.0), Region::GulfCoast); assert_eq!(for_point(30.5, -90.0), Region::GulfCoast);
} }
#[test]
fn includes_upper_boundary_points() {
// CornBelt is 38.0..48.0 lat; PacificNorthwest is 42.0..50.0.
// Half-open bounds dropped points sitting exactly at the upper
// edge into Region::Other and stripped their seasonal multiplier.
assert_ne!(for_point(48.0, -93.0), Region::Other);
assert_ne!(for_point(50.0, -122.0), Region::Other);
}
} }

View file

@ -81,7 +81,10 @@ pub fn min_refractivity_gradient(mut levels: Vec<Level>) -> Option<f64> {
let (h0, n0) = pair[0]; let (h0, n0) = pair[0];
let (h1, n1) = pair[1]; let (h1, n1) = pair[1];
let dh_km = (h1 - h0) / 1000.0; let dh_km = (h1 - h0) / 1000.0;
if dh_km.abs() < 0.01 { // Only zero-thickness layers are skipped. The previous 10 m
// guard discarded the very thin native-level layers HRRR uses
// to resolve sharp surface-based inversions and ducts.
if dh_km.abs() < 1.0e-9 {
continue; continue;
} }
let g = (n1 - n0) / dh_km; let g = (n1 - n0) / dh_km;

View file

@ -40,6 +40,14 @@ defmodule Microwaveprop.Propagation.RegionTest do
test "returns :other for points outside CONUS" do test "returns :other for points outside CONUS" do
assert Region.for_point(19.0, -155.0) == :other assert Region.for_point(19.0, -155.0) == :other
end end
test "includes points sitting exactly on the upper latitude boundary" do
# corn_belt is defined as 38.0..48.0 lat; pacific_northwest as
# 42.0..50.0. The exclusive upper bound dropped points exactly at
# 48.0 / 50.0 into `:other` and stripped their seasonal multiplier.
assert Region.for_point(48.0, -93.0) != :other
assert Region.for_point(50.0, -122.0) != :other
end
end end
describe "seasonal_adjustment/2" do describe "seasonal_adjustment/2" do

View file

@ -162,6 +162,35 @@ defmodule Microwaveprop.PropagationTest do
assert rain_with_nx == rain_no_nx assert rain_with_nx == rain_no_nx
end end
test "threads latitude into the season factor so regional adjustments fire" do
# August on the Gulf coast (lat ~29) gets the +1.15 regional multiplier
# and Iowa corn belt (lat ~42) gets the 0.80 multiplier on the same
# band/month. If score_grid_point dropped latitude on the floor the two
# season factors would be identical because the underlying base score
# only depends on `month`.
base_profile = %{
surface_temp_c: 25.0,
surface_dewpoint_c: 18.0,
surface_pressure_mb: 1013.0,
hpbl_m: 500.0,
wind_u: 3.0,
wind_v: 2.0,
cloud_cover_pct: 15.0,
precip_mm: 0.0,
min_refractivity_gradient: -100.0
}
valid_time = ~U[2026-08-15 13:00:00Z]
gulf = Propagation.score_grid_point(base_profile, valid_time, 29.0, -95.0)
iowa = Propagation.score_grid_point(base_profile, valid_time, 42.0, -93.0)
gulf_season = Enum.find(gulf, &(&1.band_mhz == 10_000)).factors.season
iowa_season = Enum.find(iowa, &(&1.band_mhz == 10_000)).factors.season
assert gulf_season > iowa_season
end
end end
describe "replace_scores/2" do describe "replace_scores/2" do

View file

@ -131,5 +131,17 @@ defmodule Microwaveprop.Radio.MaidenheadTest do
grid = Maidenhead.from_latlon(32.897, -97.038, 6) grid = Maidenhead.from_latlon(32.897, -97.038, 6)
assert String.starts_with?(grid, "EM12") assert String.starts_with?(grid, "EM12")
end end
test "clamps the upper field boundary at 90 latitude / 180 longitude" do
# Without clamping, lat = 90 + 90 = 180 → trunc(180/10) = 18, which
# the field encoder turns into 'S' — outside the legal A-R range.
# The same overflow happens at lon = 180. Both must round into the
# last valid field instead of producing an invalid grid character.
grid = Maidenhead.from_latlon(90.0, 180.0, 4)
assert Maidenhead.valid?(grid)
<<f1, f2, _, _>> = grid
assert f1 in ?A..?R
assert f2 in ?A..?R
end
end end
end end

View file

@ -175,6 +175,15 @@ defmodule Microwaveprop.RadioTest do
test "returns 0 for same point" do test "returns 0 for same point" do
assert Radio.haversine_km(32.9, -97.0, 32.9, -97.0) == 0.0 assert Radio.haversine_km(32.9, -97.0, 32.9, -97.0) == 0.0
end end
test "uses the numerically stable atan2 form" do
# Antipodal points sit on a 6371-km-radius sphere at exactly
# π·r ≈ 20_015 km apart. Both haversine forms agree here, but
# the asin form is fragile when `a` rounds to slightly above 1
# (sqrt returns >1 and asin returns NaN). atan2 stays finite.
dist = Radio.haversine_km(0.0, 0.0, 0.0, 180.0)
assert_in_delta dist, 20_015, 1.0
end
end end
describe "backfill_distances/1" do describe "backfill_distances/1" do
@ -393,6 +402,22 @@ defmodule Microwaveprop.RadioTest do
points = Radio.contact_path_points(contact) points = Radio.contact_path_points(contact)
assert length(points) == 3 assert length(points) == 3
end end
test "computes a great-circle midpoint across the anti-meridian" do
# Two stations on opposite sides of the 180° line. The arithmetic
# mean of -179 and +179 lands at 0° (Greenwich) — clearly wrong;
# the geodesic midpoint sits near the date line at ±180°.
contact =
create_contact(%{
pos1: %{"lat" => 0.0, "lon" => 179.0},
pos2: %{"lat" => 0.0, "lon" => -179.0}
})
[_, {mid_lat, mid_lon}, _] = Radio.contact_path_points(contact)
assert_in_delta mid_lat, 0.0, 0.01
assert abs(mid_lon) > 179.0
end
end end
describe "change_contact/2" do describe "change_contact/2" do

View file

@ -350,4 +350,28 @@ defmodule Microwaveprop.Weather.Grib2.Wgrib2Test do
end end
end end
end end
describe "parse_lon_val_segment/1" do
test "parses a typical decimal-formatted segment" do
assert {32.938, lon, 306.5} =
Wgrib2.parse_lon_val_segment("lon=242.958,lat=32.938,val=306.5")
# `denormalize_lon` flips the 0..360 wgrib2 form into -180..180.
assert_in_delta lon, -117.042, 0.001
end
test "tolerates a longitude string with no decimal point" do
# wgrib2 typically formats `lon=242.000` but a build that drops the
# trailing `.0` would crash `String.to_float/1` and bring the entire
# propagation chain step down. The parser should accept either form.
assert {32.938, lon, 306.5} =
Wgrib2.parse_lon_val_segment("lon=243,lat=32.938,val=306.5")
assert_in_delta lon, -117.0, 0.001
end
test "returns nil for unparseable segments" do
assert Wgrib2.parse_lon_val_segment("d=2024092219") == nil
end
end
end end

View file

@ -70,21 +70,24 @@ defmodule Microwaveprop.Weather.GridSnapPropertyTest do
end end
end end
property "output is within 30 minutes of the input" do property "output is at or before the input by at most one hour" do
check all(dt <- datetime_generator()) do check all(dt <- datetime_generator()) do
result = HrrrClient.nearest_hrrr_hour(dt) result = HrrrClient.nearest_hrrr_hour(dt)
delta_sec = abs(DateTime.diff(result, dt, :second)) # Round-down only — the result is always the published cycle
# Input minute < 30 rounds down; >= 30 rounds up. Worst-case # at or before `dt`, never a future cycle that may not exist.
# gap is exactly 30:00. delta_sec = DateTime.diff(dt, result, :second)
assert delta_sec <= 30 * 60 assert delta_sec >= 0
assert delta_sec < 60 * 60
end end
end end
property "rounds minute >= 30 up to next hour" do property "always rounds down (never points at a future cycle)" do
check all(dt <- datetime_generator(minute_range: 30..59)) do check all(dt <- datetime_generator(minute_range: 30..59)) do
result = HrrrClient.nearest_hrrr_hour(dt) result = HrrrClient.nearest_hrrr_hour(dt)
# result hour is strictly greater than input hour (modulo day/month rollover). # result is strictly before the input — the input minute > 0
assert DateTime.after?(result, dt) # ensures a non-zero offset.
assert DateTime.before?(result, dt)
assert result.hour == dt.hour
end end
end end
end end

View file

@ -22,7 +22,7 @@ defmodule Microwaveprop.Weather.HrrrClientTest do
refute HrrrClient.cycle_available?(~U[2026-04-25 00:00:00Z]) refute HrrrClient.cycle_available?(~U[2026-04-25 00:00:00Z])
end end
test "passes the rounded run_time to the probe" do test "passes the rounded-down run_time to the probe" do
test_pid = self() test_pid = self()
Application.put_env(:microwaveprop, :hrrr_cycle_available_fn, fn dt -> Application.put_env(:microwaveprop, :hrrr_cycle_available_fn, fn dt ->
@ -31,7 +31,7 @@ defmodule Microwaveprop.Weather.HrrrClientTest do
end) end)
HrrrClient.cycle_available?(~U[2026-04-25 00:38:00Z]) HrrrClient.cycle_available?(~U[2026-04-25 00:38:00Z])
assert_receive {:probed, ~U[2026-04-25 01:00:00Z]} assert_receive {:probed, ~U[2026-04-25 00:00:00Z]}
end end
end end
@ -41,9 +41,12 @@ defmodule Microwaveprop.Weather.HrrrClientTest do
~U[2026-03-28 18:00:00Z] ~U[2026-03-28 18:00:00Z]
end end
test "rounds up when past 30 minutes" do test "rounds down when past 30 minutes (never points at a future cycle)" do
# Rounding up was the historical behaviour; in production it caused
# 404/403 errors when the next cycle had not yet been published. The
# caller can always reach a later valid_time via `forecast_hour`.
assert HrrrClient.nearest_hrrr_hour(~U[2026-03-28 18:45:00Z]) == assert HrrrClient.nearest_hrrr_hour(~U[2026-03-28 18:45:00Z]) ==
~U[2026-03-28 19:00:00Z] ~U[2026-03-28 18:00:00Z]
end end
test "stays same when exactly on the hour" do test "stays same when exactly on the hour" do
@ -51,9 +54,9 @@ defmodule Microwaveprop.Weather.HrrrClientTest do
~U[2026-03-28 18:00:00Z] ~U[2026-03-28 18:00:00Z]
end end
test "rounds at exactly 30 minutes" do test "rounds down at exactly 30 minutes" do
result = HrrrClient.nearest_hrrr_hour(~U[2026-03-28 18:30:00Z]) assert HrrrClient.nearest_hrrr_hour(~U[2026-03-28 18:30:00Z]) ==
assert result in [~U[2026-03-28 18:00:00Z], ~U[2026-03-28 19:00:00Z]] ~U[2026-03-28 18:00:00Z]
end end
end end

View file

@ -89,6 +89,25 @@ defmodule Microwaveprop.Weather.SoundingParamsTest do
assert result.min_refractivity_gradient < 0 assert result.min_refractivity_gradient < 0
end end
test "captures sharp sub-10m inversions instead of skipping the layer" do
# HRRR native-level profiles routinely report 0/10/30/50 m thin layers;
# an inversion confined to the lowest 5m is exactly the kind of
# surface-based duct the propagation map cares about, so dropping any
# gradient with `dh < 10 m` quietly hides the strongest signals.
thin_layer_profile = [
%{"pres" => 1013.0, "hght" => 0.0, "tmpc" => 25.0, "dwpc" => 22.0},
%{"pres" => 1012.5, "hght" => 5.0, "tmpc" => 24.5, "dwpc" => 5.0},
%{"pres" => 1010.0, "hght" => 100.0, "tmpc" => 23.0, "dwpc" => 4.0},
%{"pres" => 950.0, "hght" => 600.0, "tmpc" => 18.0, "dwpc" => 0.0}
]
result = SoundingParams.derive(thin_layer_profile)
# The 0→5 m drop in dewpoint creates a strong negative dM/dh and
# therefore a sharp negative refractivity gradient. Without it we
# only see the broad 5→100 m baseline, which is much shallower.
assert result.min_refractivity_gradient < -200
end
test "detects temperature inversions" do test "detects temperature inversions" do
result = SoundingParams.derive(@inversion_profile) result = SoundingParams.derive(@inversion_profile)
assert result.inversions != [] assert result.inversions != []
@ -134,6 +153,24 @@ defmodule Microwaveprop.Weather.SoundingParamsTest do
assert result.ducts == [] assert result.ducts == []
end end
test "keeps a surface duct that has not yet closed by the top of the profile" do
# Three sampled levels with monotonically decreasing modified
# refractivity (M). The profile starts inside a trapping layer
# and the duct never "tops out" before the data runs out. The
# previous code threw the duct away because it never observed
# the layer close — hiding exactly the surface ducts the
# propagation map needs to surface.
open_top_duct = [
%{"pres" => 1013.0, "hght" => 0.0, "tmpc" => 30.0, "dwpc" => 27.0},
%{"pres" => 1010.0, "hght" => 20.0, "tmpc" => 33.0, "dwpc" => 5.0},
%{"pres" => 1005.0, "hght" => 50.0, "tmpc" => 35.0, "dwpc" => -10.0}
]
result = SoundingParams.derive(open_top_duct)
assert result.ducting_detected
assert result.ducts != []
end
test "computes boundary layer depth" do test "computes boundary layer depth" do
result = SoundingParams.derive(@inversion_profile) result = SoundingParams.derive(@inversion_profile)
# BL depth should be positive and reasonable # BL depth should be positive and reasonable

View file

@ -129,6 +129,40 @@ defmodule Microwaveprop.WeatherTest do
assert %{weather_status: :queued} = Repo.get(Contact, contact.id) assert %{weather_status: :queued} = Repo.get(Contact, contact.id)
end end
test "scales the longitude band by latitude so high-lat windows still cover ~150 km east/west" do
# At lat 49° (US/Canada border) cos(lat) ≈ 0.656, so a fixed 1.5°
# longitude box collapses to ~110 km — narrower than the 150 km
# window the function documents. A station 1.8° east of the QSO
# is roughly 132 km away, well inside the intended radius, but a
# latitude-blind bounding box would exclude it.
{:ok, station} =
Weather.find_or_create_station(%{
@station_attrs
| station_code: "KHIGH",
lat: 49.0,
lon: -121.7
})
Weather.upsert_surface_observation(station, %{
observed_at: ~U[2023-09-17 17:45:00Z],
temp_f: 60.0
})
contact =
create_contact(%{
station1: "HILAT",
station2: "HILAT2",
pos1: %{"lat" => 49.0, "lon" => -123.5},
pos2: nil,
grid2: nil
})
Radio.set_enrichment_status!([contact.id], :weather_status, :queued)
{:ok, _n} = Weather.reconcile_weather_statuses()
assert %{weather_status: :complete} = Repo.get(Contact, contact.id)
end
test "leaves :complete and :pending contacts alone" do test "leaves :complete and :pending contacts alone" do
{:ok, station} = Weather.find_or_create_station(@station_attrs) {:ok, station} = Weather.find_or_create_station(@station_attrs)