rao 0.2.2

Robust and scalable Adaptive Optics tools
Documentation 
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300

use std::ops;

use serde::{Deserialize, Serialize};

/// 2D geometric vector, associated with transverse plane of optical
/// path.
///
/// This struct is used for defining displacements, coordinates, and directions
/// along the transverse plane (in a slice of the propogation axis).
/// # Examples
/// ```
/// let a = rao::Vec2D::new(3.0, 4.0);
/// assert_eq!(a.norm(),5.0);
/// let b = &a * 2.0;
/// assert_eq!(b.norm(),10.0);
/// let e_x = rao::Vec2D::x_unit();
/// let e_y = rao::Vec2D::y_unit();
/// assert_eq!(e_x.dot(&e_y), 0.0);
/// assert_eq!(e_x.dot(&e_x), 1.0);
/// assert_eq!(e_x.dot(&a), 3.0);
/// assert_eq!((e_x / 2.0).x, 0.5);
/// ```
#[derive(Debug, Clone, PartialEq, Serialize, Deserialize)]
pub struct Vec2D {
    pub x: f64,
    pub y: f64,
}

impl Vec2D {
    /// Create a new [`Vec2D`] from the coordinates in the transverse plane.
    #[must_use]
    pub fn new(x: f64, y: f64) -> Self {
        Self { x, y }
    }
    /// create a unit vector along x-axis.
    #[must_use]
    pub fn x_unit() -> Self {
        Self { x: 1.0, y: 0.0 }
    }
    /// create a unit vector along y-axis.
    #[must_use]
    pub fn y_unit() -> Self {
        Self { x: 0.0, y: 1.0 }
    }
    /// return the (Euclidean) norm of the [`Vec2D`].
    #[must_use]
    pub fn norm(&self) -> f64 {
        self.norm2().powf(0.5)
    }
    /// return the (Euclidean) squared-norm of the [`Vec2D`].
    #[must_use]
    pub fn norm2(&self) -> f64 {
        self.x.powf(2.0) + self.y.powf(2.0)
    }
    /// return the dot (inner) product of the [`Vec2D`] with another (borrowed) [`Vec2D`].
    #[must_use]
    pub fn dot(&self, other: &Self) -> f64 {
        self.x * other.x + self.y * other.y
    }
    /// Calculate a uniformly spaced set of points between two [`Vec2D`]s,
    /// including the specified endpoints
    #[must_use]
    pub fn linspace(a: &Self, b: &Self, npoints: u32) -> Vec<Self> {
        match npoints {
            0 => vec![],
            1 => vec![a.clone()],
            _ => (0..npoints)
                .map(|u| f64::from(u) / (f64::from(npoints) - 1.0))
                .map(|t| (1.0 - t) * a + t * b)
                .collect(),
        }
    }
    /// Calculate the centre of equal sized intervals between two [`Vec2D`]s
    #[must_use]
    pub fn linspread(a: &Self, b: &Self, npoints: u32) -> Vec<Self> {
        (0..npoints)
            .map(|u| (f64::from(u) / f64::from(npoints)) + 1.0 / (2.0 * f64::from(npoints)))
            .map(|t| (1.0 - t) * a + t * b)
            .collect()
    }
    /// Calculates a [`Vec2D`] that is rotated by +90 degrees, such that it is
    /// orthogonal to the input [`Vec2D`].
    #[must_use]
    pub fn ortho(&self) -> Self {
        Self {
            x: -self.y,
            y: self.x,
        }
    }
}

impl_op!(-|a: Vec2D| -> Vec2D { Vec2D::new(-a.x, -a.y) });

impl_op_ex_commutative!(/ |a:&Vec2D,b:&f64| -> Vec2D {
    Vec2D {
        x: a.x / b,
        y: a.y / b,
    }
});
impl_op_ex_commutative!(*|a: &Vec2D, b: &f64| -> Vec2D {
    Vec2D {
        x: a.x * b,
        y: a.y * b,
    }
});
impl_op_ex!(+ |a:&Vec2D,b:&Vec2D| -> Vec2D {
    Vec2D {
        x: a.x + b.x,
        y: a.y + b.y,
    }
});
impl_op_ex!(-|a: &Vec2D, b: &Vec2D| -> Vec2D {
    Vec2D {
        x: a.x - b.x,
        y: a.y - b.y,
    }
});

impl_op_ex_commutative!(+
    |a:&Vec3D,b:&Vec2D| -> Vec3D
    {
        Vec3D {x:a.x+b.x,y:a.y+b.y,z:a.z}
    }
);
impl_op_ex_commutative!(+
    |a:&Line,b:&Vec2D| -> Line
    {
        Line::new(a.x0+b.x, a.xz, a.y0+b.y, a.yz)
    }
);

/// 3D geometric vector, associated with both transverse (x,y) and
/// propagation (z) dimensions.
///
/// This struct is used for defining coordinates of optical components that
/// have a meaninful position along the optical propagation axis (e.g., an actuator
/// of a deformable mirror conjugated to some altitude). In particular this
/// struct is useful for finding the intersection of an optical ray (see [Line])
/// with the influence function of an actuator where that [Line] intersects the
/// altitude plane of the actuator.
/// # Examples
/// ```
/// // `point` at x=6.0 m, y=3.0 m, z=10.0 m
/// let point = rao::Vec3D::new(6.0, 3.0, 10.0);
/// // `line` with z-component of x = 1.0 metres/metre.
/// let line = rao::Line::new(0.0, 1.0, 0.0, 0.0);
/// // when `line` propagates to z=10.0 m, it will have
/// // x=10.0 m, y=0.0 m, which is +4.0 metres away from
/// // the `point` in x, and -3.0 metres away in y.
/// assert_eq!(
///     point.displacement_at_altitude(&line),
///     rao::Vec2D::new(4.0, -3.0)
/// );
/// assert_eq!(
///     point.distance_at_altitude(&line),
///     5.0
/// );
/// ```
#[derive(Debug, Clone, PartialEq, Serialize, Deserialize)]
pub struct Vec3D {
    /// Position in x relative to the optical axis
    pub x: f64,
    /// Position in y relative to the optical axis
    pub y: f64,
    /// Position in z relative to the optical axis
    /// e.g., 0km => pupil-plane, +infinity => object-plane
    pub z: f64,
}

impl Vec3D {
    /// Create a new [`Vec3D`] from a 3D coordinate.
    #[must_use]
    pub fn new(x: f64, y: f64, z: f64) -> Self {
        Self { x, y, z }
    }
    /// Convenience function for creating a [`Vec3D`] at the origin.
    #[must_use]
    pub fn origin() -> Self {
        Self {
            x: 0.0,
            y: 0.0,
            z: 0.0,
        }
    }
    /// Return the Euclidean norm of the displacement vector between self
    /// and the intersection of a [Line] at the altitude of self.
    #[must_use]
    pub fn distance_at_altitude(&self, line: &Line) -> f64 {
        self.displacement_at_altitude(line).norm()
    }
    /// Return the displacement vector between self
    /// and the intersection of a [Line] at the altitude of self.
    #[must_use]
    pub fn displacement_at_altitude(&self, line: &Line) -> Vec2D {
        let line_intersection = line.position_at_altitude(self.z);
        Vec2D {
            x: (line_intersection.x - self.x),
            y: (line_intersection.y - self.y),
        }
    }
}

/// Line that propagates along the optical axis.
///
/// The parametric form of the line as defined here is convenient for
/// finding the intersection of the line with a plane parallel with the
/// x-y plane, that is, points of constant altitude.
///
/// The `(x, y, z)` tuple is parameterised by `t`:
/// ```math
/// x = x0 + xz*t
/// y = y0 + yz*t
/// z = t
/// ```
#[derive(Debug, Clone, PartialEq, Serialize, Deserialize)]
pub struct Line {
    /// The x-position at z=0.
    pub x0: f64,
    /// The component of the x-position which depends linearly on z.
    pub xz: f64,
    /// The y-position at z=0.
    pub y0: f64,
    /// The component of the y-position which depends linearly on z.
    pub yz: f64,
}

impl Line {
    /// Constructor for a new [Line].
    #[must_use]
    pub fn new(x0: f64, xz: f64, y0: f64, yz: f64) -> Line {
        Line { x0, xz, y0, yz }
    }
    /// Convenience function for [Line]s that do not depend on z
    #[must_use]
    pub fn new_on_axis(x0: f64, y0: f64) -> Line {
        Line::new(x0, 0.0, y0, 0.0)
    }
    /// Convenience function for defining a [Line] given two points
    /// in 3D. Results in singularities if a.z == b.z, otherwise is
    /// gauaranteed to be stable.
    #[must_use]
    pub fn new_from_two_points(a: &Vec3D, b: &Vec3D) -> Line {
        let xz = (b.x - a.x) / (b.z - a.z);
        let yz = (b.y - a.y) / (b.z - a.z);
        let x0 = a.x - (a.z) / (b.z - a.z) * (b.x - a.x);
        let y0 = a.y - (a.z) / (b.z - a.z) * (b.y - a.y);
        Line::new(x0, xz, y0, yz)
    }
    /// Calculate the [`Vec2D`] coordinates in the transverse plane at
    /// a specified altitude.
    #[must_use]
    pub fn position_at_altitude(&self, alt: f64) -> Vec2D {
        Vec2D::new(alt * self.xz + self.x0, alt * self.yz + self.y0)
    }

    #[must_use]
    pub fn distance_at_ground(&self, other: &Line) -> f64 {
        Vec2D::new(self.x0 - other.x0, self.y0 - other.y0).norm()
    }
}

#[cfg(test)]
mod tests {
    use super::*;
    use approx::assert_abs_diff_eq;

    #[test]
    fn linspread() {
        let a = Vec2D::new(0.0, 3.0);
        let b = Vec2D::new(3.0, 0.0);
        let ls = Vec2D::linspread(&a, &b, 3);
        assert_abs_diff_eq!(ls[0].x, 0.5);
        assert_abs_diff_eq!(ls[0].y, 2.5);
        assert_abs_diff_eq!(ls[1].x, 1.5);
        assert_abs_diff_eq!(ls[1].y, 1.5);
        assert_abs_diff_eq!(ls[2].x, 2.5);
        assert_abs_diff_eq!(ls[2].y, 0.5);
    }

    #[test]
    fn linspace() {
        let a = Vec2D::new(1.0, 2.0);
        let b = Vec2D::new(4.0, 0.0);
        let ls = Vec2D::linspace(&a, &b, 3);
        assert_abs_diff_eq!(ls[1].x, 2.5);
        assert_abs_diff_eq!(ls[1].y, 1.0);
        assert_eq!(ls[0], a);
        assert_eq!(ls[2], b);
    }

    #[test]
    fn adding() {
        let a = Vec2D::new(1.0, 2.0);
        let b = Vec3D::new(10.0, 20.0, 30.0);
        let c = a + b;
        assert_abs_diff_eq!(c.x, 11.0);
        assert_abs_diff_eq!(c.y, 22.0);
        assert_abs_diff_eq!(c.z, 30.0);
    }
}