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