gmt_dos_clients_lom/
lib.rs

1//! # M1 & M2 Rigid Body Motions to Linear Optical Model
2//!
3//! Transforms M1 and M2 rigid body motions to optical metrics
4//! (tip-tilt, segment piston and  segment tip-tilt) using
5//! linear optical sensitivity matrices
6
7use std::{io::Read, path::Path, sync::Arc};
8
9use flate2::bufread::GzDecoder;
10use gmt_dos_clients_io::{
11    gmt_m1::M1RigidBodyMotions,
12    gmt_m2::{asm::M2ASMReferenceBodyNodes, M2RigidBodyMotions},
13    optics::{
14        MaskedWavefront, SegmentD21PistonRSS, SegmentPiston, SegmentTipTilt, SegmentWfeRms,
15        TipTilt, Wavefront, WfeRms,
16    },
17};
18use gmt_lom::{LinearOpticalModelError, Loader, LOM};
19use interface::{self, Data, Size, Units, Update, Write};
20
21mod optical_sensitivity;
22pub use optical_sensitivity::OpticalSensitivities;
23
24/// M1 & M2 Rigid Body Motions to Linear Optical Model
25#[derive(Debug, Clone)]
26pub struct LinearOpticalModel {
27    lom: LOM,
28    m1_rbm: Arc<Vec<f64>>,
29    m2_rbm: Arc<Vec<f64>>,
30}
31impl LinearOpticalModel {
32    pub fn new() -> std::result::Result<Self, LinearOpticalModelError> {
33        let sens = include_bytes!("optical_sensitivities.rs.bin.gz");
34        let mut gz = GzDecoder::new(sens.as_slice());
35        let mut bytes = vec![];
36        gz.read_to_end(&mut bytes)?;
37        Ok(Self {
38            lom: LOM::try_from(bytes.as_slice())?,
39            m1_rbm: Arc::new(vec![0f64; 42]),
40            m2_rbm: Arc::new(vec![0f64; 42]),
41        })
42    }
43    pub fn new_with_path(
44        path: impl AsRef<Path>,
45    ) -> std::result::Result<Self, LinearOpticalModelError> {
46        let path = path.as_ref();
47        let filename = path.file_name().unwrap();
48        let sens_loader = Loader::<gmt_lom::OpticalSensitivities>::default()
49            .path(path.parent().unwrap())
50            .filename(filename.to_str().unwrap());
51        Ok(Self {
52            lom: LOM::builder()
53                .load_optical_sensitivities(sens_loader)?
54                .build()?,
55            m1_rbm: Arc::new(vec![0f64; 42]),
56            m2_rbm: Arc::new(vec![0f64; 42]),
57        })
58    }
59}
60
61impl Units for LinearOpticalModel {}
62
63impl Update for LinearOpticalModel {
64    fn update(&mut self) {
65        self.lom.rbm = vec![(self.m1_rbm.as_slice(), self.m2_rbm.as_slice())]
66            .into_iter()
67            .collect();
68    }
69}
70impl interface::Read<M1RigidBodyMotions> for LinearOpticalModel {
71    fn read(&mut self, data: Data<M1RigidBodyMotions>) {
72        self.m1_rbm = data.into_arc();
73    }
74}
75impl interface::Read<M2RigidBodyMotions> for LinearOpticalModel {
76    fn read(&mut self, data: Data<M2RigidBodyMotions>) {
77        self.m2_rbm = data.into_arc();
78    }
79}
80impl interface::Read<M2ASMReferenceBodyNodes> for LinearOpticalModel {
81    fn read(&mut self, data: Data<M2ASMReferenceBodyNodes>) {
82        self.m2_rbm = data.into_arc();
83    }
84}
85
86impl Write<TipTilt> for LinearOpticalModel {
87    fn write(&mut self) -> Option<Data<TipTilt>> {
88        Some(Data::new(self.lom.tiptilt().into()))
89    }
90}
91impl Write<SegmentTipTilt> for LinearOpticalModel {
92    fn write(&mut self) -> Option<Data<SegmentTipTilt>> {
93        Some(Data::new(self.lom.segment_tiptilt().into()))
94    }
95}
96impl<const E: i32> Write<SegmentPiston<E>> for LinearOpticalModel {
97    fn write(&mut self) -> Option<Data<SegmentPiston<E>>> {
98        let mut piston: Vec<f64> = self.lom.segment_piston().into();
99        if E != 0 {
100            piston.iter_mut().for_each(|p| {
101                *p *= 10f64.powi(-E);
102            });
103        }
104        Some(piston.into())
105    }
106}
107impl<const E: i32> Write<SegmentWfeRms<E>> for LinearOpticalModel {
108    fn write(&mut self) -> Option<Data<SegmentWfeRms<E>>> {
109        Some(self.lom.segment_wfe_rms::<E>().into())
110    }
111}
112impl<const E: i32> Write<SegmentD21PistonRSS<E>> for LinearOpticalModel {
113    fn write(&mut self) -> Option<Data<SegmentD21PistonRSS<E>>> {
114        let piston: Vec<f64> = self.lom.segment_piston().into();
115        let mut sum_squared = 0f64;
116        for p_i in piston.iter().take(6) {
117            for p_j in piston.iter().skip(1) {
118                sum_squared += (p_i - p_j).powi(2);
119            }
120        }
121        let mut rss = sum_squared.sqrt();
122        if E != 0 {
123            rss *= 10f64.powi(-E);
124        }
125        Some(vec![rss].into())
126    }
127}
128
129impl Write<Wavefront> for LinearOpticalModel {
130    fn write(&mut self) -> Option<Data<Wavefront>> {
131        Some(Data::new(self.lom.wavefront().into()))
132    }
133}
134
135impl Write<MaskedWavefront> for LinearOpticalModel {
136    fn write(&mut self) -> Option<Data<MaskedWavefront>> {
137        Some(Data::new(self.lom.masked_wavefront().into()))
138    }
139}
140
141impl<const E: i32> Write<WfeRms<E>> for LinearOpticalModel {
142    fn write(&mut self) -> Option<Data<WfeRms<E>>> {
143        let wavefront = self.lom.masked_wavefront();
144        let n = wavefront.len() as f64;
145        let (s, m) = wavefront
146            .into_iter()
147            .fold((0f64, 0f64), |(mut s, mut m), w| {
148                s += w * w;
149                m += w;
150                (s, m)
151            });
152        let mut wfe_rms = ((s - m * m / n) / n).sqrt();
153        if E != 0 {
154            wfe_rms *= 10f64.powi(-E);
155        }
156        Some(Data::new(vec![wfe_rms]))
157    }
158}
159
160impl Size<TipTilt> for LinearOpticalModel {
161    fn len(&self) -> usize {
162        2
163    }
164}
165impl Size<SegmentTipTilt> for LinearOpticalModel {
166    fn len(&self) -> usize {
167        14
168    }
169}
170impl Size<SegmentPiston> for LinearOpticalModel {
171    fn len(&self) -> usize {
172        7
173    }
174}
175impl Size<Wavefront> for LinearOpticalModel {
176    fn len(&self) -> usize {
177        512 * 512
178    }
179}
180impl Size<WfeRms> for LinearOpticalModel {
181    fn len(&self) -> usize {
182        1
183    }
184}
185
186#[cfg(test)]
187mod tests {
188    use super::*;
189
190    #[test]
191    fn m1_segment_rxy() {
192        let mut rbm2lom = LinearOpticalModel::new().unwrap();
193        let s = 1e-6;
194        let m1_rbm: Vec<_> = (0..7)
195            .flat_map(|_| {
196                let mut v = vec![0f64; 6];
197                v[3] = s;
198                v
199            })
200            .collect();
201        <LinearOpticalModel as interface::Read<M1RigidBodyMotions>>::read(
202            &mut rbm2lom,
203            m1_rbm.into(),
204        );
205        rbm2lom.update();
206        let mut stt: Vec<f64> = <LinearOpticalModel as Write<SegmentTipTilt>>::write(&mut rbm2lom)
207            .unwrap()
208            .into();
209        stt.iter_mut().for_each(|x| *x /= s);
210        let stt_mag: Vec<_> = stt[..7]
211            .iter()
212            .zip(&stt[7..])
213            .map(|(x, y)| x.hypot(*y))
214            .collect();
215        println!("STT: {:.2?}", stt_mag);
216    }
217
218    #[test]
219    fn m2_segment_rxy() {
220        let mut rbm2lom = LinearOpticalModel::new().unwrap();
221        let s = 1e-6;
222        let m1_rbm: Vec<_> = (0..7)
223            .flat_map(|_| {
224                let mut v = vec![0f64; 6];
225                v[4] = s;
226                v
227            })
228            .collect();
229        <LinearOpticalModel as interface::Read<M2RigidBodyMotions>>::read(
230            &mut rbm2lom,
231            m1_rbm.into(),
232        );
233        rbm2lom.update();
234        let mut stt: Vec<f64> = <LinearOpticalModel as Write<SegmentTipTilt>>::write(&mut rbm2lom)
235            .unwrap()
236            .into();
237        stt.iter_mut().for_each(|x| *x /= s);
238        let stt_mag: Vec<_> = stt[..7]
239            .iter()
240            .zip(&stt[7..])
241            .map(|(x, y)| x.hypot(*y))
242            .collect();
243        println!("STT: {:.2?}", stt_mag);
244    }
245
246    #[test]
247    fn m1_segment_txy() {
248        let mut rbm2lom = LinearOpticalModel::new().unwrap();
249        let s = 1e-6;
250        let m1_rbm: Vec<_> = (0..7)
251            .flat_map(|_| {
252                let mut v = vec![0f64; 6];
253                v[0] = s;
254                v
255            })
256            .collect();
257        <LinearOpticalModel as interface::Read<M1RigidBodyMotions>>::read(
258            &mut rbm2lom,
259            m1_rbm.into(),
260        );
261        rbm2lom.update();
262        let mut stt: Vec<f64> = <LinearOpticalModel as Write<SegmentTipTilt>>::write(&mut rbm2lom)
263            .unwrap()
264            .into();
265        stt.iter_mut().for_each(|x| *x /= s);
266        let stt_mag: Vec<_> = stt[..7]
267            .iter()
268            .zip(&stt[7..])
269            .map(|(x, y)| x.hypot(*y))
270            .collect();
271        println!("STT: {:.3?}", stt_mag);
272    }
273}
gmt_dos_clients_lom/lib.rs

gmt_dos_clients_lom/
lib.rs
