cescobar
9d7d3bac45
Fork de SPLAT! (John A. Magliacane, KD2BD, 2002-2014), bajo GPLv2 heredada del original. El núcleo de cálculo de propagación (ITWOM v3.0, Sid Shumate) se mantiene sin modificar en cpp/itwom3.0.cpp, llamado vía FFI. El resto del pipeline (lectura de formatos SPLAT!, reportes, mapas, KML, gnuplot) se reescribió en Rust, con paralelización nativa vía rayon reemplazando el paralelismo por múltiples procesos del original. Validado contra el corpus golden de SPLAT! con drift < 0.02 dB en cálculos ITWOM.
418 lines
16 KiB
Rust
418 lines
16 KiB
Rust
//! Parallel azimuthal (radial) propagation sweep.
|
||
//!
|
||
//! The expensive step of a path-loss map is computing the signal attenuation
|
||
//! along every radial from the transmitter out to the coverage radius. SPLAT!'s
|
||
//! original C++ does this with a single-threaded sweep from 0° to 360°; the
|
||
//! work is embarrassingly parallel (each radial is independent of the others
|
||
//! once the DEM is loaded). Rayon + thread-local ITWOM state let us scale
|
||
//! nearly linearly with CPU count.
|
||
//!
|
||
//! Output strategy: each radial returns a `Vec<PixelWrite>` and the driver
|
||
//! merges them after the parallel sweep. Writes never contend. When two
|
||
//! radials hit the same pixel (near the TX), the merge keeps the *lower* loss
|
||
//! (strongest signal) — matching SPLAT!'s sequential overwrite semantics with
|
||
//! `PutSignal`.
|
||
|
||
use rayon::prelude::*;
|
||
|
||
use crate::geo::{self, EARTH_RADIUS_MI, FEET_PER_METER, METERS_PER_MILE};
|
||
use crate::itwom::{self, LrParams, Model};
|
||
use crate::types::{Config, Site};
|
||
|
||
/// A single pixel write produced by a radial: geographic coordinate plus
|
||
/// path loss in dB. Aggregated by the map builder.
|
||
#[derive(Debug, Clone, Copy)]
|
||
pub struct PixelWrite {
|
||
pub lat: f64,
|
||
pub lon: f64,
|
||
pub loss_db: f64,
|
||
}
|
||
|
||
/// Per-radial context passed to the compute function. Holds *read-only*
|
||
/// references so Rayon can share it freely across threads.
|
||
pub struct RadialCtx<'a> {
|
||
pub config: &'a Config,
|
||
pub source: &'a Site,
|
||
/// RX antenna height in meters AGL — the "probe" used to compute loss
|
||
/// at every radial sample point, matching SPLAT!'s `-L altitude` flag.
|
||
pub rx_alt_m: f64,
|
||
pub max_range_km: f64,
|
||
pub params: &'a LrParams,
|
||
pub model: Model,
|
||
pub dem: &'a dyn DemLookup,
|
||
}
|
||
|
||
/// Minimal interface the radial sweep needs from the DEM. Concrete impl
|
||
/// will back this with an in-memory SDF tile set. Defining the trait here
|
||
/// keeps `compute_radial` testable with synthetic terrain.
|
||
pub trait DemLookup: Sync {
|
||
/// Returns elevation in meters, or some large negative sentinel for
|
||
/// no-data. Callers should treat z < -1000 as "outside DEM coverage".
|
||
fn elevation(&self, lat: f64, lon: f64) -> f64;
|
||
}
|
||
|
||
/// Compute one radial from the transmitter outward along `azimuth` degrees.
|
||
///
|
||
/// Pure function — no shared mutable state, safe to call concurrently. This
|
||
/// is the unit of work `plot_lr_map_parallel` distributes across cores.
|
||
///
|
||
/// Algorithm (port of SPLAT!'s `PlotLRPath` driver):
|
||
/// 1. Walk outward along the great-circle bearing in steps of
|
||
/// `1 / samples_per_radian` radians, sampling the DEM at each step.
|
||
/// 2. For every sample past the 3rd, call ITWOM with the *prefix* profile
|
||
/// (samples 0..=i). That gives path loss from TX to the point at step `i`.
|
||
/// 3. Emit a `PixelWrite` per outer point with the loss in dB.
|
||
///
|
||
/// Complexity per radial is O(N²) in ITWOM calls (one per step, with the
|
||
/// call size proportional to the step index). Matches SPLAT!'s cost model.
|
||
pub fn compute_radial(ctx: &RadialCtx<'_>, azimuth_deg: f64) -> Vec<PixelWrite> {
|
||
let max_d_mi = ctx.max_range_km / 1.609344;
|
||
let spr = geo::samples_per_radian(ctx.config.ppd());
|
||
let path_length = (max_d_mi / EARTH_RADIUS_MI) * spr;
|
||
let n = path_length.max(2.0) as usize;
|
||
let miles_per_sample = max_d_mi / path_length;
|
||
let delta_m = miles_per_sample * METERS_PER_MILE;
|
||
|
||
// Sample the terrain profile once along the full radial.
|
||
let mut lats = Vec::with_capacity(n + 1);
|
||
let mut lons = Vec::with_capacity(n + 1);
|
||
let mut elev_m = Vec::with_capacity(n + 1);
|
||
for i in 0..=n {
|
||
let d_mi = miles_per_sample * i as f64;
|
||
let (lat, lon) =
|
||
geo::step_great_circle(ctx.source.lat, ctx.source.lon, azimuth_deg, d_mi);
|
||
let z = ctx.dem.elevation(lat, lon);
|
||
if z < -1000.0 { break; }
|
||
lats.push(lat);
|
||
lons.push(lon);
|
||
elev_m.push(z);
|
||
}
|
||
|
||
if elev_m.len() < 4 {
|
||
return Vec::new(); // too short to do any meaningful path loss
|
||
}
|
||
|
||
let source_alt_m = ctx.source.alt;
|
||
// Reused elev buffer for ITWOM — first two slots are header (npts-1, delta).
|
||
let mut elev_buf = Vec::with_capacity(elev_m.len() + 2);
|
||
let mut out = Vec::with_capacity(elev_m.len().saturating_sub(2));
|
||
|
||
for i in 2..elev_m.len() {
|
||
elev_buf.clear();
|
||
elev_buf.push((i - 1) as f64); // number_of_points - 1
|
||
elev_buf.push(delta_m);
|
||
elev_buf.extend_from_slice(&elev_m[0..=i]);
|
||
|
||
// Hot path: we only need the dB number, not the mode label.
|
||
let loss_db = itwom::path_loss_db(
|
||
ctx.model,
|
||
&mut elev_buf,
|
||
source_alt_m,
|
||
ctx.rx_alt_m,
|
||
ctx.params,
|
||
);
|
||
out.push(PixelWrite {
|
||
lat: lats[i],
|
||
lon: lons[i],
|
||
loss_db,
|
||
});
|
||
}
|
||
let _ = FEET_PER_METER; // silence unused import warning if feet conv unused
|
||
out
|
||
}
|
||
|
||
/// LOS-only radial — no ITWOM calls. Checks whether each point along the
|
||
/// radial has line of sight to the TX, taking Earth curvature (4/3 R) and
|
||
/// antenna heights into account. Emits `loss_db = 0.0` for visible points
|
||
/// and `loss_db = f64::INFINITY` for obstructed ones, so downstream
|
||
/// `CoverageMap::ingest` (which keeps the *min* loss per pixel) produces
|
||
/// a binary visibility raster.
|
||
///
|
||
/// Dramatically cheaper than the full path-loss sweep — use it as a quick
|
||
/// "where would this transmitter even be seen from" pass before paying
|
||
/// for the ITWOM computation.
|
||
pub fn compute_los_radial(ctx: &RadialCtx<'_>, azimuth_deg: f64) -> Vec<PixelWrite> {
|
||
let max_d_mi = ctx.max_range_km / 1.609344;
|
||
let spr = geo::samples_per_radian(ctx.config.ppd());
|
||
let path_length = (max_d_mi / EARTH_RADIUS_MI) * spr;
|
||
let n = path_length.max(2.0) as usize;
|
||
let miles_per_sample = max_d_mi / path_length;
|
||
|
||
// 4/3-Earth effective radius in meters, for curvature correction.
|
||
let four_thirds_earth_m: f64 = (4.0 / 3.0) * 20_902_230.97 * 0.3048;
|
||
let tx_ground_m = ctx.dem.elevation(ctx.source.lat, ctx.source.lon);
|
||
let tx_top_m = tx_ground_m + ctx.source.alt;
|
||
|
||
// Compare *slopes* (rise/run) rather than cosines — slopes at different
|
||
// distances compose correctly. The curvature-corrected apparent height
|
||
// of a point at distance d is `z - d²/(2R)`, so the slope from the TX
|
||
// antenna top is `(z - tx_top)/d - d/(2R)`.
|
||
let slope = |z: f64, d_m: f64| (z - tx_top_m) / d_m - d_m / (2.0 * four_thirds_earth_m);
|
||
|
||
let mut out = Vec::with_capacity(n);
|
||
let mut max_terrain_slope = f64::NEG_INFINITY;
|
||
for i in 1..n {
|
||
let d_mi = miles_per_sample * i as f64;
|
||
let (lat, lon) =
|
||
geo::step_great_circle(ctx.source.lat, ctx.source.lon, azimuth_deg, d_mi);
|
||
let z = ctx.dem.elevation(lat, lon);
|
||
if z < -1000.0 { break; }
|
||
|
||
let d_m = d_mi * METERS_PER_MILE;
|
||
let s_terrain = slope(z, d_m);
|
||
if s_terrain > max_terrain_slope {
|
||
max_terrain_slope = s_terrain;
|
||
}
|
||
let s_rx = slope(z + ctx.rx_alt_m, d_m);
|
||
|
||
let visible = s_rx >= max_terrain_slope;
|
||
out.push(PixelWrite {
|
||
lat,
|
||
lon,
|
||
loss_db: if visible { 0.0 } else { f64::INFINITY },
|
||
});
|
||
}
|
||
out
|
||
}
|
||
|
||
/// Parallel 360° LOS coverage sweep. Same shape as `plot_lr_map_parallel`
|
||
/// but uses `compute_los_radial` — typically 100-500× faster because no
|
||
/// ITWOM invocations are needed.
|
||
pub fn plot_los_coverage_parallel(ctx: &RadialCtx<'_>) -> Vec<PixelWrite> {
|
||
let ppd = ctx.config.ppd() as f64;
|
||
let n_steps = (360.0 * ppd) as u32;
|
||
(0..n_steps)
|
||
.into_par_iter()
|
||
.map(|i| i as f64 / ppd)
|
||
.flat_map_iter(|az| compute_los_radial(ctx, az))
|
||
.filter(|p| p.loss_db.is_finite())
|
||
.collect()
|
||
}
|
||
|
||
/// Sweep the full 360° azimuth range in parallel, using one radial per
|
||
/// `1.0 / ppd` degree step (matching SPLAT!'s native resolution).
|
||
///
|
||
/// The call blocks until all radials complete. Returns flat pixel writes
|
||
/// that the caller folds into a coverage map.
|
||
pub fn plot_lr_map_parallel(ctx: &RadialCtx<'_>) -> Vec<PixelWrite> {
|
||
let ppd = ctx.config.ppd() as f64;
|
||
let n_steps = (360.0 * ppd) as u32;
|
||
|
||
(0..n_steps)
|
||
.into_par_iter()
|
||
.map(|i| i as f64 / ppd)
|
||
.flat_map_iter(|az| compute_radial(ctx, az))
|
||
.collect()
|
||
}
|
||
|
||
/// Same sweep as `plot_lr_map_parallel`, but keeps the per-radial grouping
|
||
/// so the caller can rasterize consecutive samples as line segments
|
||
/// (eliminating diagonal-radial aliasing artifacts in the output map).
|
||
///
|
||
/// Use `CoverageMap::from_radials(&result, ppd, pad)` to render without
|
||
/// streak gaps.
|
||
pub fn plot_lr_map_parallel_radials(ctx: &RadialCtx<'_>) -> Vec<Vec<PixelWrite>> {
|
||
let ppd = ctx.config.ppd() as f64;
|
||
let n_steps = (360.0 * ppd) as u32;
|
||
|
||
(0..n_steps)
|
||
.into_par_iter()
|
||
.map(|i| i as f64 / ppd)
|
||
.map(|az| compute_radial(ctx, az))
|
||
.collect()
|
||
}
|
||
|
||
/// Sequential baseline of the same sweep, used for correctness diffing and
|
||
/// benchmark comparison against the parallel version.
|
||
pub fn plot_lr_map_sequential(ctx: &RadialCtx<'_>) -> Vec<PixelWrite> {
|
||
let ppd = ctx.config.ppd() as f64;
|
||
let n_steps = (360.0 * ppd) as u32;
|
||
|
||
(0..n_steps)
|
||
.map(|i| i as f64 / ppd)
|
||
.flat_map(|az| compute_radial(ctx, az))
|
||
.collect()
|
||
}
|
||
|
||
#[cfg(test)]
|
||
mod tests {
|
||
use super::*;
|
||
use crate::sdf::DemTileSet;
|
||
|
||
struct FlatDem(f64);
|
||
impl DemLookup for FlatDem {
|
||
fn elevation(&self, _lat: f64, _lon: f64) -> f64 { self.0 }
|
||
}
|
||
|
||
#[test]
|
||
fn single_radial_returns_monotonic_growing_loss() {
|
||
// Flat Earth: loss should increase (mostly monotonically) with distance.
|
||
let cfg = Config::default();
|
||
let src = Site { name: "X".into(), lat: -1.0, lon: 79.0, alt: 30.0 };
|
||
let p = LrParams {
|
||
eps_dielect: 15.0,
|
||
sgm_conductivity: 0.005,
|
||
eno_ns_surfref: 301.0,
|
||
frq_mhz: 426.0,
|
||
radio_climate: 5,
|
||
pol: 1,
|
||
conf: 0.5,
|
||
rel: 0.5,
|
||
};
|
||
let dem = FlatDem(50.0);
|
||
let ctx = RadialCtx {
|
||
config: &cfg,
|
||
source: &src,
|
||
rx_alt_m: 10.0,
|
||
max_range_km: 20.0,
|
||
params: &p,
|
||
model: Model::Itwom,
|
||
dem: &dem,
|
||
};
|
||
let r = compute_radial(&ctx, 0.0);
|
||
assert!(!r.is_empty());
|
||
// Spot-check: the last point should have more loss than the first.
|
||
let first = r.first().unwrap().loss_db;
|
||
let last = r.last().unwrap().loss_db;
|
||
assert!(last > first, "expected loss to grow: first={first} last={last}");
|
||
for p in &r {
|
||
assert!(p.loss_db.is_finite(), "non-finite: {:?}", p);
|
||
}
|
||
}
|
||
|
||
/// Sweeps a small fraction of the 360° range with real Ecuador DEM,
|
||
/// comparing sequential vs parallel. Small radius so both variants
|
||
/// finish quickly — the full map is proved out by
|
||
/// `full_rvt_queve_parallel_only` below.
|
||
#[test]
|
||
#[ignore = "expensive; run with --ignored"]
|
||
fn rvt_queve_quadrant_ser_vs_par() {
|
||
let tiles_dir = "/home/cescobar/Laboratorio/splat-1.4.2-Charles/sdf";
|
||
if !std::path::Path::new(tiles_dir).exists() {
|
||
eprintln!("skipping: {tiles_dir} not present");
|
||
return;
|
||
}
|
||
let mut dem = DemTileSet::new();
|
||
let loaded = dem.load_dir(tiles_dir).unwrap();
|
||
println!("loaded {loaded} tiles");
|
||
|
||
let src = Site { name: "Q".into(), lat: -1.023206, lon: 79.458669, alt: 14.0 };
|
||
let p = LrParams {
|
||
eps_dielect: 15.0, sgm_conductivity: 0.005, eno_ns_surfref: 301.0,
|
||
frq_mhz: 426.0, radio_climate: 5, pol: 1, conf: 0.5, rel: 0.5,
|
||
};
|
||
let cfg = Config::default();
|
||
|
||
// 10 km radius, first 90° only — ~108 k radials × ~118 samples each.
|
||
let ppd = cfg.ppd() as f64;
|
||
let n = (90.0 * ppd) as u32;
|
||
let build_ctx = |range| RadialCtx {
|
||
config: &cfg, source: &src, rx_alt_m: 10.0, max_range_km: range,
|
||
params: &p, model: Model::Itwom, dem: &dem,
|
||
};
|
||
let ctx = build_ctx(10.0);
|
||
|
||
let quadrant_par = |ctx: &RadialCtx<'_>| -> Vec<PixelWrite> {
|
||
use rayon::prelude::*;
|
||
(0..n).into_par_iter()
|
||
.map(|i| i as f64 / ppd)
|
||
.flat_map_iter(|az| compute_radial(ctx, az))
|
||
.collect()
|
||
};
|
||
let quadrant_seq = |ctx: &RadialCtx<'_>| -> Vec<PixelWrite> {
|
||
(0..n)
|
||
.map(|i| i as f64 / ppd)
|
||
.flat_map(|az| compute_radial(ctx, az))
|
||
.collect()
|
||
};
|
||
|
||
let t0 = std::time::Instant::now();
|
||
let px_seq = quadrant_seq(&ctx);
|
||
let d_seq = t0.elapsed();
|
||
|
||
let t1 = std::time::Instant::now();
|
||
let px_par = quadrant_par(&ctx);
|
||
let d_par = t1.elapsed();
|
||
|
||
println!("pixels: seq={} par={}", px_seq.len(), px_par.len());
|
||
println!("sequential: {:?}", d_seq);
|
||
println!("parallel: {:?}", d_par);
|
||
println!("speedup: {:.2}x", d_seq.as_secs_f64() / d_par.as_secs_f64());
|
||
|
||
assert_eq!(px_seq.len(), px_par.len(),
|
||
"quadrant sweep produced different pixel counts");
|
||
}
|
||
|
||
/// Compute a single radial toward Cerro Cochabamba and check the path
|
||
/// loss at that distance matches SPLAT!'s reported 145.5 dB. This is our
|
||
/// primary correctness gate for the pipeline.
|
||
#[test]
|
||
#[ignore = "needs real DEM"]
|
||
fn rvt_queve_cerro_cochabamba_matches_splat() {
|
||
let tiles_dir = "/home/cescobar/Laboratorio/splat-1.4.2-Charles/sdf";
|
||
if !std::path::Path::new(tiles_dir).exists() { return; }
|
||
let mut dem = DemTileSet::new();
|
||
dem.load_dir(tiles_dir).unwrap();
|
||
|
||
let src = Site { name: "Q".into(), lat: -1.023206, lon: 79.458669, alt: 14.0 };
|
||
let p = LrParams {
|
||
eps_dielect: 15.0, sgm_conductivity: 0.005, eno_ns_surfref: 301.0,
|
||
frq_mhz: 426.0, radio_climate: 5, pol: 1, conf: 0.5, rel: 0.5,
|
||
};
|
||
let cfg = Config::default();
|
||
let ctx = RadialCtx {
|
||
config: &cfg, source: &src, rx_alt_m: 20.0, max_range_km: 90.0,
|
||
params: &p, model: Model::Itwom, dem: &dem,
|
||
};
|
||
let radial = compute_radial(&ctx, 152.52);
|
||
// SPLAT reports 84.66 km to Cerro Cochabamba. Find the sample closest
|
||
// to that distance along this radial.
|
||
let d = |p: &PixelWrite| crate::geo::distance(
|
||
src.lat, src.lon, p.lat, p.lon) * 1.609344;
|
||
let target_km = 84.66;
|
||
let best = radial.iter()
|
||
.min_by(|a, b| (d(a) - target_km).abs().total_cmp(&(d(b) - target_km).abs()))
|
||
.expect("non-empty radial");
|
||
println!(
|
||
"Cerro Cochabamba sample: {:.2} km, loss {:.2} dB (SPLAT golden: 145.50 dB)",
|
||
d(best), best.loss_db
|
||
);
|
||
// The SPLAT C++ report uses RX AGL=20 m at Cerro Cochabamba; our
|
||
// radial uses a uniform rx_alt_m=20. Matching within 5 dB is a
|
||
// generous tolerance that still catches gross errors — bit-for-bit
|
||
// equivalence would need same SDF no-data handling and interpolation.
|
||
assert!(
|
||
(best.loss_db - 145.5).abs() < 5.0,
|
||
"loss drifted too far from SPLAT golden: {:.2} dB", best.loss_db
|
||
);
|
||
}
|
||
|
||
/// Full 360°/40 km sweep, parallel only — validates the end-to-end
|
||
/// pipeline under realistic workload. Aim: <3 min on a 16-core box.
|
||
#[test]
|
||
#[ignore = "expensive; run with --ignored"]
|
||
fn full_rvt_queve_parallel_only() {
|
||
let tiles_dir = "/home/cescobar/Laboratorio/splat-1.4.2-Charles/sdf";
|
||
if !std::path::Path::new(tiles_dir).exists() {
|
||
eprintln!("skipping: {tiles_dir} not present");
|
||
return;
|
||
}
|
||
let mut dem = DemTileSet::new();
|
||
dem.load_dir(tiles_dir).unwrap();
|
||
let src = Site { name: "Q".into(), lat: -1.023206, lon: 79.458669, alt: 14.0 };
|
||
let p = LrParams {
|
||
eps_dielect: 15.0, sgm_conductivity: 0.005, eno_ns_surfref: 301.0,
|
||
frq_mhz: 426.0, radio_climate: 5, pol: 1, conf: 0.5, rel: 0.5,
|
||
};
|
||
let cfg = Config::default();
|
||
let ctx = RadialCtx {
|
||
config: &cfg, source: &src, rx_alt_m: 10.0, max_range_km: 40.73,
|
||
params: &p, model: Model::Itwom, dem: &dem,
|
||
};
|
||
let t = std::time::Instant::now();
|
||
let px = plot_lr_map_parallel(&ctx);
|
||
println!("full sweep: {} pixels in {:?}", px.len(), t.elapsed());
|
||
assert!(px.len() > 100_000_000);
|
||
}
|
||
}
|