1use crate::field::{ProfileField, WindField};
10use crate::sounding::Sounding;
11
12const LAYER_TOPS_M: [f64; 4] = [1_000.0, 3_000.0, 6_000.0, 10_000.0];
15
16#[derive(Debug, Clone, Copy, PartialEq)]
18pub struct LayerError {
19 pub top_m: f64,
22 pub mean_abs_speed_error_mps: f64,
24 pub levels: usize,
26}
27
28#[derive(Debug, Clone, PartialEq)]
30pub struct ValidationReport {
31 pub levels: usize,
33 pub speed_rmse_mps: f64,
35 pub direction_rmse_deg: f64,
37 pub layers: Vec<LayerError>,
39}
40
41pub fn validate(profile: &ProfileField, sounding: &Sounding) -> ValidationReport {
45 let mut speed_sq = 0.0;
46 let mut direction_sq = 0.0;
47 let mut layer_abs = vec![0.0; LAYER_TOPS_M.len() + 1];
48 let mut layer_n = vec![0usize; LAYER_TOPS_M.len() + 1];
49 for level in &sounding.levels {
50 let wind = profile.wind(0.0, [0.0, level.altitude_m, 0.0]);
51 let speed_err = wind.speed_mps() - level.speed_mps;
52 speed_sq += speed_err * speed_err;
53 let mut dir_err = (wind.direction_from_deg() - level.direction_from_deg) % 360.0;
54 if dir_err > 180.0 {
55 dir_err -= 360.0;
56 } else if dir_err < -180.0 {
57 dir_err += 360.0;
58 }
59 direction_sq += dir_err * dir_err;
60 let layer = LAYER_TOPS_M
61 .iter()
62 .position(|top| level.altitude_m < *top)
63 .unwrap_or(LAYER_TOPS_M.len());
64 layer_abs[layer] += speed_err.abs();
65 layer_n[layer] += 1;
66 }
67 let n = sounding.levels.len().max(1) as f64;
68 ValidationReport {
69 levels: sounding.levels.len(),
70 speed_rmse_mps: (speed_sq / n).sqrt(),
71 direction_rmse_deg: (direction_sq / n).sqrt(),
72 layers: layer_abs
73 .into_iter()
74 .zip(layer_n)
75 .enumerate()
76 .filter(|(_, (_, count))| *count > 0)
77 .map(|(i, (abs, count))| LayerError {
78 top_m: LAYER_TOPS_M.get(i).copied().unwrap_or(f64::INFINITY),
79 mean_abs_speed_error_mps: abs / count as f64,
80 levels: count,
81 })
82 .collect(),
83 }
84}