parse_monitors/
lib.rs

1//! GMT Computation Fluid Dynamic Post Processing API
2//!
3//! A library to query and to process the GMT CFD Baseline databases.
4//!
5//! The path to the CFD database is set with the env variable: `CFD_REPO`
6
7use nalgebra as na;
8
9pub mod error;
10mod vector;
11pub use vector::Vector;
12mod monitors;
13pub use monitors::{Exertion, Mirror, Monitors, MonitorsError, MonitorsLoader};
14pub mod cfd;
15pub mod domeseeing;
16pub use domeseeing::{Band, Data, DomeSeeing};
17pub mod pressure;
18pub mod report;
19pub mod temperature;
20
21#[cfg(all(feature = "2020", not(feature = "default")))]
22pub const CFD_YEAR: u32 = 2020;
23#[cfg(all(feature = "2021", not(feature = "default")))]
24pub const CFD_YEAR: u32 = 2021;
25#[cfg(feature = "default")]
26pub const CFD_YEAR: u32 = 2025;
27
28pub const FORCE_SAMPLING_FREQUENCY: f64 = 20_f64; // Hz
29pub const FORCE_SAMPLING: f64 = 1. / FORCE_SAMPLING_FREQUENCY; // Hz
30pub const TEMPERATURE_SAMPLING_FREQUENCY: f64 = 5_f64; // Hz
31
32pub fn polyfit<T: na::RealField + Copy>(
33    x_values: &[T],
34    y_values: &[T],
35    polynomial_degree: usize,
36) -> Result<Vec<T>, &'static str> {
37    let number_of_columns = polynomial_degree + 1;
38    let number_of_rows = x_values.len();
39    let mut a = na::DMatrix::zeros(number_of_rows, number_of_columns);
40
41    for (row, &x) in x_values.iter().enumerate() {
42        // First column is always 1
43        a[(row, 0)] = T::one();
44
45        for col in 1..number_of_columns {
46            a[(row, col)] = x.powf(na::convert(col as f64));
47        }
48    }
49
50    let b = na::DVector::from_row_slice(y_values);
51
52    let decomp = na::SVD::new(a, true, true);
53
54    match decomp.solve(&b, na::convert(1e-18f64)) {
55        Ok(mat) => Ok(mat.data.into()),
56        Err(error) => Err(error),
57    }
58}
59pub fn detrend<T: na::RealField + Copy>(
60    x_values: &[T],
61    y_values: &[T],
62    polynomial_degree: usize,
63) -> Result<Vec<T>, &'static str> {
64    let number_of_columns = polynomial_degree + 1;
65    let number_of_rows = x_values.len();
66    let mut a = na::DMatrix::zeros(number_of_rows, number_of_columns);
67
68    for (row, &x) in x_values.iter().enumerate() {
69        // First column is always 1
70        a[(row, 0)] = T::one();
71
72        for col in 1..number_of_columns {
73            a[(row, col)] = x.powf(na::convert(col as f64));
74        }
75    }
76
77    let b = na::DVector::from_row_slice(y_values);
78
79    let decomp = na::SVD::new(a.clone(), true, true);
80
81    match decomp.solve(&b, na::convert(1e-18f64)) {
82        Ok(mat) => {
83            let y_detrend = b - a * &mat;
84            Ok(y_detrend.data.into())
85        }
86        Err(error) => Err(error),
87    }
88}
89pub fn detrend_mut<T: na::RealField + Copy>(
90    x_values: &[T],
91    y_values: &mut [T],
92    polynomial_degree: usize,
93) -> Result<(), &'static str> {
94    let number_of_columns = polynomial_degree + 1;
95    let number_of_rows = x_values.len();
96    let mut a = na::DMatrix::zeros(number_of_rows, number_of_columns);
97
98    for (row, &x) in x_values.iter().enumerate() {
99        // First column is always 1
100        a[(row, 0)] = T::one();
101
102        for col in 1..number_of_columns {
103            a[(row, col)] = x.powf(na::convert(col as f64));
104        }
105    }
106
107    let b = na::DVector::from_row_slice(y_values);
108
109    let decomp = na::SVD::new(a.clone(), true, true);
110
111    match decomp.solve(&b, na::convert(1e-18f64)) {
112        Ok(mat) => {
113            let mut y_detrend = b - a * &mat;
114            y_values.swap_with_slice(&mut y_detrend.as_mut_slice());
115            Ok(())
116        }
117        Err(error) => Err(error),
118    }
119}
120
121#[cfg(feature = "plot")]
122pub fn plot_monitor<S: AsRef<std::path::Path> + std::convert::AsRef<std::ffi::OsStr>>(
123    time: &[f64],
124    monitor: &[Exertion],
125    key: &str,
126    path: S,
127) {
128    use plotters::prelude::*;
129    let max_value = |x: &[f64]| -> f64 {
130        x.iter()
131            .cloned()
132            .rev()
133            .take(400 * 20)
134            .fold(std::f64::NEG_INFINITY, f64::max)
135    };
136    let min_value = |x: &[f64]| -> f64 {
137        x.iter()
138            .cloned()
139            .rev()
140            .take(400 * 20)
141            .fold(std::f64::INFINITY, f64::min)
142    };
143
144    let file_path = std::path::Path::new(&path).join("TOTAL_FORCES.png");
145    let filename = if let Some(filename) = file_path.to_str() {
146        filename.to_string()
147    } else {
148        eprintln!("Invalid file path: {:?}", file_path);
149        return;
150    };
151    let plot = BitMapBackend::new(&filename, (768, 512)).into_drawing_area();
152    plot.fill(&WHITE).unwrap();
153
154    let (min_value, max_value) = {
155        let force_magnitude: Option<Vec<f64>> =
156            monitor.iter().map(|e| e.force.magnitude()).collect();
157        (
158            min_value(force_magnitude.as_ref().unwrap()),
159            max_value(force_magnitude.as_ref().unwrap()),
160        )
161    };
162    let xrange = time.last().unwrap() - time[0];
163    let minmax_padding = 0.1;
164    let mut chart = ChartBuilder::on(&plot)
165        .set_label_area_size(LabelAreaPosition::Left, 60)
166        .set_label_area_size(LabelAreaPosition::Bottom, 40)
167        .margin(10)
168        .build_cartesian_2d(
169            -xrange * 1e-2..xrange * (1. + 1e-2),
170            min_value * (1. - minmax_padding)..max_value * (1. + minmax_padding),
171        )
172        .unwrap();
173    chart
174        .configure_mesh()
175        .x_desc("Time [s]")
176        .y_desc(format!("{} Force [N]", key))
177        .draw()
178        .unwrap();
179
180    let mut colors = colorous::TABLEAU10.iter().cycle();
181
182    let color = colors.next().unwrap();
183    let rgb = RGBColor(color.r, color.g, color.b);
184    chart
185        .draw_series(LineSeries::new(
186            time.iter()
187                .zip(monitor.iter())
188                .map(|(&x, y)| (x - time[0], y.force.magnitude().unwrap())),
189            &rgb,
190        ))
191        .unwrap()
192        .label(key)
193        .legend(move |(x, y)| PathElement::new(vec![(x, y), (x + 20, y)], &rgb));
194
195    chart
196        .configure_series_labels()
197        .border_style(&BLACK)
198        .background_style(&WHITE.mix(0.8))
199        .position(SeriesLabelPosition::UpperRight)
200        .draw()
201        .unwrap();
202}
203
204#[cfg(test)]
205mod tests {
206    use std::path::PathBuf;
207
208    use crate::cfd::{Baseline, BaselineTrait, CfdCase};
209
210    use super::*;
211    /*
212       use nalgebra as na;
213       #[test]
214       fn test_arm() {
215           let force = [100f64, -33f64, 250f64];
216           let force_v = na::Vector3::from_column_slice(&force);
217           //let arm = na::Vector3::<f64>::new_random() * 2f64 - na::Vector3::repeat(1f64);
218           let arm = na::Vector3::<f64>::from_column_slice(&[1., 1., 1.]);
219           println!("ARM: {:?}", arm);
220           let moment = arm.cross(&force_v);
221           println!("Moment: {:?}", moment);
222           let amat = na::Matrix3::new(
223               0., force[2], -force[1], -force[2], 0., force[0], force[1], -force[0], 0.,
224           );
225           println!("A: {:#?}", amat);
226           println!("Moment: {:?}", amat * arm);
227           let qr = amat.svd(true, true);
228           let x = qr.solve(&moment, 1e-3).unwrap();
229           println!("ARM: {:?}", x);
230           println!("Moment: {:?}", x.cross(&force_v));
231       }
232    */
233    #[test]
234    fn monitors() -> anyhow::Result<()> {
235        let cfd_case = CfdCase::<2025>::colloquial(30, 0, "os", 7)?;
236        let path = Baseline::<2025>::path()?.join(cfd_case.to_string());
237        println!("{path:?}");
238        let mon = Monitors::loader::<PathBuf, 2025>(path).load()?;
239        println!(
240            "Monitors F&M key #: {:}",
241            mon.forces_and_moments.keys().collect::<Vec<_>>().len()
242        );
243        Ok(())
244    }
245    #[test]
246    fn cfd_2020() {
247        let monitors = MonitorsLoader::<2020>::default()
248            .data_path("/fsx/Baseline2020/b2019_30z_0az_os_7ms/")
249            .header_filter("Total".to_string())
250            .load()
251            .unwrap();
252        println!(
253            "Time: {:.3?}s",
254            (monitors.time[0], monitors.time.last().unwrap())
255        );
256        println!("Force entries #: {}", monitors.forces_and_moments.len());
257        monitors
258            .forces_and_moments
259            .keys()
260            .for_each(|k| println!("Key: {}", k));
261        println!(
262            "Total force entries #: {}",
263            monitors.total_forces_and_moments.len()
264        );
265    }
266    /*
267        #[test]
268        fn load_mirror_table() {
269            let mut m1 = Mirror::m1();
270            m1.load(
271                "/fsx/Baseline2021/Baseline2021/Baseline2021/CASES/zen00az180_OS2",
272                false,
273            )
274            .unwrap();
275            let t = m1.time().front().unwrap();
276            let f: Vec<_> = m1
277                .forces_and_moments()
278                .filter_map(|f| f.front().map(|v| v.force.clone()))
279                .collect();
280            println!("{}: {:?}", t, f);
281            let t = m1.time().back().unwrap();
282            let f: Vec<_> = m1
283                .forces_and_moments()
284                .filter_map(|f| f.back().map(|v| v.force.clone()))
285                .collect();
286            println!("{}: {:?}", t, f);
287        }
288    */
289    #[test]
290    fn test_polyfit() {
291        let (a, b) = (-1.5f64, 5f64);
292        let (x, y): (Vec<_>, Vec<_>) = (0..10).map(|k| (k as f64, a * k as f64 + b)).unzip();
293        let ba = polyfit(&x, &y, 1).unwrap();
294        println!("ab: {:#?}", ba);
295        assert!((ba[0] - b).abs() < 1e-6 && (ba[1] - a).abs() < 1e-6)
296    }
297    #[test]
298    fn test_detrend() {
299        let (a, b) = (-1.5f64, 5f64);
300        let (x, y): (Vec<_>, Vec<_>) = (0..10).map(|k| (k as f64, a * k as f64 + b)).unzip();
301        let ydtd = detrend(&x, &y, 1).unwrap();
302        let ba = polyfit(&x, &ydtd, 1).unwrap();
303        println!("ab: {:#?}", ba);
304        //assert!((ba[0] - b).abs() < 1e-6 && (ba[1] - a).abs() < 1e-6)
305    }
306    #[test]
307    fn test_detrend_mut() {
308        let (a, b) = (-1.5f64, 5f64);
309        let (x, mut y): (Vec<_>, Vec<_>) = (0..10).map(|k| (k as f64, a * k as f64 + b)).unzip();
310        detrend_mut(&x, &mut y, 1).unwrap();
311        let ba = polyfit(&x, &y, 1).unwrap();
312        println!("ab: {:#?}", ba);
313        //assert!((ba[0] - b).abs() < 1e-6 && (ba[1] - a).abs() < 1e-6)
314    }
315}
parse_monitors/lib.rs

parse_monitors/
lib.rs
