crseo/
atmosphere.rs

1use serde::{Deserialize, Serialize};
2use std::{
3    f32,
4    ffi::CString,
5    fmt::Display,
6    ops::{Div, Mul},
7    path::Path,
8};
9
10use crate::builders::atmosphere::{AtmosphereBuilder, AtmosphereBuilderError};
11
12use super::{Cu, FromBuilder, Propagation, Single, Source};
13use ffi::atmosphere;
14
15#[derive(Debug, thiserror::Error)]
16pub enum AtmosphereError {
17    #[error("cannot create `::crseo::AtmosphereBuilder`")]
18    Builder(#[from] AtmosphereBuilderError),
19}
20pub type Result<T> = std::result::Result<T, AtmosphereError>;
21
22#[allow(dead_code)]
23#[derive(Deserialize, Debug)]
24#[serde(rename = "lower_case")]
25struct GmtAtmosphere {
26    r0: f32,
27    #[serde(rename = "L0")]
28    l_not: f32,
29    #[serde(rename = "L")]
30    length: f32,
31    nxy_pupil: i32,
32    fov: f32,
33    duration: f32,
34    n_duration: i32,
35    filename: String,
36    seed: i32,
37}
38
39#[derive(Debug, Clone, PartialEq, Serialize, Deserialize)]
40#[doc(hidden)]
41pub struct TurbulenceProfile {
42    pub n_layer: usize,
43    pub altitude: Vec<f32>,
44    pub xi0: Vec<f32>,
45    pub wind_speed: Vec<f32>,
46    pub wind_direction: Vec<f32>,
47}
48impl Default for TurbulenceProfile {
49    fn default() -> Self {
50        TurbulenceProfile {
51            n_layer: 7,
52            altitude: [25.0, 275.0, 425.0, 1_250.0, 4_000.0, 8_000.0, 13_000.0].to_vec(),
53            xi0: [0.1257, 0.0874, 0.0666, 0.3498, 0.2273, 0.0681, 0.0751].to_vec(),
54            wind_speed: [5.6540, 5.7964, 5.8942, 6.6370, 13.2925, 34.8250, 29.4187].to_vec(),
55            wind_direction: [0.0136, 0.1441, 0.2177, 0.5672, 1.2584, 1.6266, 1.7462].to_vec(),
56        }
57    }
58}
59#[derive(Debug, Clone, PartialEq, Serialize, Deserialize)]
60#[doc(hidden)]
61pub struct RayTracing {
62    pub width: f32,
63    pub n_width_px: i32,
64    pub field_size: f32,
65    pub duration: f32,
66    pub filepath: Option<String>,
67    pub n_duration: Option<i32>,
68}
69/// Default properties:
70///  * width        : 25.5m
71///  * width_px     : 512px
72///  * field_size   : 0rd
73///  * duration     : 1s
74///  * filepath     : None
75///  * n_duration   : None
76impl Default for RayTracing {
77    fn default() -> Self {
78        Self {
79            width: 25.5,
80            n_width_px: 512,
81            field_size: 0.0,
82            duration: 1.0,
83            filepath: None,
84            n_duration: None,
85        }
86    }
87}
88impl RayTracing {
89    /// Size in meters of the phase screen at altitude 0m
90    pub fn width(mut self, width: f64) -> Self {
91        self.width = width as f32;
92        self
93    }
94    /// Size in pixels of the phase screen at altitude 0m
95    pub fn n_width_px(mut self, n_width_px: usize) -> Self {
96        self.n_width_px = n_width_px as i32;
97        self
98    }
99    /// Field-of-view in radians
100    pub fn field_size(mut self, field_size: f64) -> Self {
101        self.field_size = field_size as f32;
102        self
103    }
104    /// Phase screen minimum time length in seconds
105    ///
106    /// Phase screens of that time length must fit with the GPU memory
107    pub fn duration(mut self, duration: f64) -> Self {
108        self.duration = duration as f32;
109        self
110    }
111    /// Path where to write the phase screens data file
112    pub fn filepath<P: AsRef<Path>>(mut self, filepath: P) -> Self {
113        let path = filepath.as_ref();
114        self.filepath = Some(path.to_str().unwrap().to_string());
115        self
116    }
117    /// Total number of durations
118    ///
119    /// The total time length of the phase screens is `n_duration X duration` seconds
120    pub fn n_duration(mut self, n_duration: u64) -> Self {
121        self.n_duration = Some(n_duration as i32);
122        self
123    }
124}
125
126pub struct Atmosphere {
127    pub(crate) _c_: atmosphere,
128    pub r0_at_zenith: f64,
129    pub oscale: f64,
130    pub zenith_angle: f64,
131    pub secs: f64,
132    //filename: String,
133    //k_duration: i32,
134    pub(crate) propagate_ptr: fn(&mut Atmosphere, &mut Source, f32),
135}
136impl Display for Atmosphere {
137    fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result {
138        writeln!(
139            f,
140            "Atmosphere: r0@{:.0}z={:.2}cm, outer scale={:.2}m",
141            self.zenith_angle.to_degrees(),
142            self.r0() * 1e2,
143            self.oscale,
144        )
145    }
146}
147impl FromBuilder for Atmosphere {
148    type ComponentBuilder = AtmosphereBuilder;
149}
150impl Atmosphere {
151    /*
152        pub fn new() -> Atmosphere {
153            Atmosphere {
154                _c_: unsafe { mem::zeroed() },
155                r0_at_zenith: 0.16,
156                oscale: 25.5,
157                zenith_angle: 0.0,
158                secs: 0.0,
159                //filename: String::new(),
160                //k_duration: 0,
161                propagate_ptr: |_, _, _| (),
162            }
163        }
164        pub fn build(
165            &mut self,
166            r_not: f32,
167            l_not: f32,
168            n_layer: i32,
169            mut altitude: Vec<f32>,
170            mut xi0: Vec<f32>,
171            mut wind_speed: Vec<f32>,
172            mut wind_direction: Vec<f32>,
173        ) -> &mut Self {
174            unsafe {
175                self._c_.setup(
176                    r_not,
177                    l_not,
178                    n_layer,
179                    altitude.as_mut_ptr(),
180                    xi0.as_mut_ptr(),
181                    wind_speed.as_mut_ptr(),
182                    wind_direction.as_mut_ptr(),
183                );
184            }
185            self.propagate_ptr = |a, s, t| unsafe {
186                let n_xy = s.pupil_sampling;
187                let d_xy = (s.pupil_size / (n_xy - 1) as f64) as f32;
188                a._c_
189                    .get_phase_screen4(s.as_raw_mut_ptr(), d_xy, n_xy, d_xy, n_xy, t);
190            };
191            self
192        }
193    */
194    pub fn as_raw_mut_ptr(&mut self) -> &mut atmosphere {
195        &mut self._c_
196    }
197    pub fn raytrace_build(
198        &mut self,
199        r_not: f32,
200        l_not: f32,
201        n_layer: i32,
202        mut altitude: Vec<f32>,
203        mut xi0: Vec<f32>,
204        mut wind_speed: Vec<f32>,
205        mut wind_direction: Vec<f32>,
206        width: f32,
207        n_width_px: i32,
208        field_size: f32,
209        duration: f32,
210        filepath: Option<&str>,
211        n_duration: Option<i32>,
212    ) -> &mut Self {
213        match filepath {
214            Some(file) => unsafe {
215                self._c_.setup2(
216                    r_not,
217                    l_not,
218                    n_layer,
219                    altitude.as_mut_ptr(),
220                    xi0.as_mut_ptr(),
221                    wind_speed.as_mut_ptr(),
222                    wind_direction.as_mut_ptr(),
223                    width,
224                    n_width_px,
225                    field_size,
226                    duration,
227                    CString::new(file.to_owned().into_bytes())
228                        .unwrap()
229                        .into_raw(),
230                    n_duration.unwrap_or(1),
231                );
232            },
233            None => unsafe {
234                self._c_.setup1(
235                    r_not,
236                    l_not,
237                    n_layer,
238                    altitude.as_mut_ptr(),
239                    xi0.as_mut_ptr(),
240                    wind_speed.as_mut_ptr(),
241                    wind_direction.as_mut_ptr(),
242                    width,
243                    n_width_px,
244                    field_size,
245                    duration,
246                );
247            },
248        }
249        self.propagate_ptr = |a, s, t| unsafe {
250            let n_xy = s.pupil_sampling;
251            let d_xy = (s.pupil_size / (n_xy - 1) as f64) as f32;
252            a._c_
253                .rayTracing1(s.as_raw_mut_ptr(), d_xy, n_xy, d_xy, n_xy, t);
254        };
255        self
256    }
257    pub fn gmt_build(&mut self, r_not: f32, l_not: f32) -> &mut Self {
258        unsafe {
259            self._c_.gmt_setup4(r_not, l_not, 2020);
260        }
261        self
262    }
263    pub fn get_phase_values<'a, T>(
264        &mut self,
265        src: &mut Source,
266        x: &'a [T],
267        y: &'a [T],
268        t: f64,
269    ) -> Vec<T>
270    where
271        &'a [T]: Into<Cu<Single>>,
272        Cu<Single>: Into<Vec<T>>,
273        T: 'a,
274    {
275        let n = x.len();
276        let mut gx: Cu<Single> = x.into();
277        let mut gy: Cu<Single> = y.into();
278        let mut ps = Cu::<Single>::vector(n);
279        ps.malloc();
280        unsafe {
281            self._c_.get_phase_screen(
282                ps.as_mut_ptr(),
283                gx.as_mut_ptr(),
284                gy.as_mut_ptr(),
285                n as i32,
286                src.as_raw_mut_ptr(),
287                t as f32,
288            )
289        }
290        ps.into()
291    }
292    pub fn get_phase_screen<'a, T>(
293        &mut self,
294        src: &mut Source,
295        t: f64,
296        (s_x, n_x): (T, usize),
297        other_side: Option<(T, usize)>,
298    ) -> Vec<T>
299    where
300        T: 'a + Copy + From<u32> + Div<Output = T> + Mul<Output = T>,
301        Vec<T>: Into<Cu<Single>>,
302        Cu<Single>: Into<Vec<T>>,
303    {
304        let (s_y, n_y) = other_side.unwrap_or((s_x, n_x));
305        let n = n_x * n_y;
306        let mut x: Vec<T> = Vec::with_capacity(n);
307        let mut y: Vec<T> = Vec::with_capacity(n);
308        let delta_x = s_x / T::try_from(n_x as u32 - 1).unwrap();
309        let delta_y = s_y / T::try_from(n_x as u32 - 1).unwrap();
310        for i in 0..n_x {
311            for j in 0..n_y {
312                x.push(delta_x * T::try_from(i as u32).unwrap());
313                y.push(delta_y * T::try_from(j as u32).unwrap());
314            }
315        }
316        let mut gx: Cu<Single> = x.into();
317        let mut gy: Cu<Single> = y.into();
318        let mut ps = Cu::<Single>::vector(n);
319        ps.malloc();
320        unsafe {
321            self._c_.get_phase_screen(
322                ps.as_mut_ptr(),
323                gx.as_mut_ptr(),
324                gy.as_mut_ptr(),
325                n as i32,
326                src.as_raw_mut_ptr(),
327                t as f32,
328            )
329        }
330        ps.into()
331    }
332    pub fn update_r0(&mut self, new_r0: f64) {
333        self._c_.r0 = new_r0 as f32;
334    }
335    pub fn r0(&self) -> f64 {
336        let secz = 1f64 / self.zenith_angle.cos();
337        (self.r0_at_zenith.powf(-5.0 / 3.0) * secz).powf(-3.0 / 5.0)
338    }
339    pub fn reset(&mut self) {
340        unsafe {
341            self._c_.reset();
342        }
343    }
344}
345impl Drop for Atmosphere {
346    fn drop(&mut self) {
347        unsafe {
348            self._c_.cleanup();
349        }
350    }
351}
352impl Propagation for Atmosphere {
353    fn time_propagate(&mut self, secs: f64, src: &mut Source) {
354        (self.propagate_ptr)(self, src, secs as f32);
355    }
356    fn propagate(&mut self, src: &mut Source) {
357        self.time_propagate(self.secs, src)
358    }
359}
360
361#[cfg(test)]
362mod tests {
363    use super::*;
364
365    // cargo test --release --package crseo --lib  -- atmosphere::tests::atmosphere_new --exact --nocapture
366    #[test]
367    fn atmosphere_new() {
368        crate::ceo!(Atmosphere);
369    }
370
371    // cargo test --release --package crseo --lib  -- atmosphere::tests::dump_toml --exact --nocapture
372    #[test]
373    fn dump_toml() -> anyhow::Result<()> {
374        let builder = AtmosphereBuilder::default().ray_tracing(Default::default());
375        builder.save("atm_builder.toml")?;
376        Ok(())
377    }
378
379    // cargo test --release --package crseo --lib  -- atmosphere::tests::load_toml --exact --nocapture
380    #[test]
381    fn load_toml() -> anyhow::Result<()> {
382        let builder = AtmosphereBuilder::load("atm_builder.toml")?;
383        dbg!(&builder);
384        Ok(())
385    }
386}
crseo/atmosphere.rs

crseo/
atmosphere.rs
