refraction/
lib.rs

1// Copyright (C) The Strand-Braid Authors
2// SPDX-License-Identifier: MIT OR Apache-2.0
3
4use bisection_search::{BisectionSearch, Interval};
5use simba::scalar::RealField;
6
7/// Define the parameters required to solve a refraction problem.
8///
9/// In the following drawing, we want to find `x` when we know `w`, `h`, `d`,
10/// and the ratio of refractive indices `n` where `n = n2/n1`.
11///
12/// `a` and `b` are angles, `n1` and `n2` are refractive indices, `d`, `h` and
13/// `w` are lengths.
14///
15/// ```text
16///                                                      ---
17///                                                 ----/ a|
18///               n1                           ----/       |
19///                                       ----/            |
20///                                 -----/                 |
21///                            ----/                       | w
22///                       ----/                            |
23///                  ----/                                 |
24///           x   --/                    d-x               |
25///      |--------|----------------------------------------|
26///      |       /
27///      | h     |
28///      |      /
29///      |     /
30///      |    /                   n2
31///      |    |
32///      |   /
33///      |  /
34///      | /
35///      |b|
36///      |/
37///
38///    sin(a)   n2
39///    ------ = -- = n
40///    sin(b)   n1
41/// ```
42///
43/// See [this](https:///math.stackexchange.com/questions/150769) for an
44/// explanation of the problem and this solution. Ascii drawing done with
45/// [textik.com](https://textik.com).
46#[derive(Clone)]
47pub struct RefractionEq<T> {
48    /// Distance along the refractive boundary
49    pub d: T,
50    /// Height below the refractive boundary
51    pub h: T,
52    /// Height above the refractive boundary
53    pub w: T,
54    /// Ratio of refractive indices
55    pub n: T,
56}
57
58impl<T: RealField + Copy> RefractionEq<T> {
59    /// Evaluate the refraction equation at location `x`.
60    pub fn f(&self, x: T) -> T {
61        let RefractionEq { d, h, w, n } = self.clone();
62
63        let d_minus_x = d - x;
64        d_minus_x * (x * x + h * h).sqrt() / (x * (d_minus_x * d_minus_x + w * w).sqrt()) - n
65    }
66}
67
68pub fn find_root<T>(a: T, b: T, eq: RefractionEq<T>, tolerance: T) -> Option<T>
69where
70    T: RealField + Copy,
71{
72    let interval = Interval::new(a, b)?;
73
74    let mut bisect = BisectionSearch::new(interval, |x| eq.f(*x));
75    loop {
76        bisect = bisect.step();
77        if bisect.interval.size() < tolerance {
78            break;
79        }
80    }
81    Some(*bisect.interval.a())
82}
83
84#[cfg(test)]
85mod tests {
86    #[test]
87    fn example1() {
88        // https://math.stackexchange.com/questions/150769/refraction-equation-quartic-equation/151249#comment347929_151050
89        // https://www.wolframalpha.com/input/?i=plot+%28%28d-x%29%2Fsqrt%28%28d-x%29%5E2%2Bw%5E2%29%29%2F%28x%2Fsqrt%28x%5E2%2Bh%5E2%29%29+-+1.33+where+w%3D1%2C+h%3D5%2C+d%3D30+for+x%3D0..30
90        let f = crate::RefractionEq {
91            d: 30.0,
92            h: 5.0,
93            w: 1.0,
94            n: 1.33,
95        };
96        let eps = 0.01;
97        assert!(f64::abs(f.f(5.7)) < eps);
98    }
99
100    #[test]
101    fn test_derivative() {
102        let refraction = crate::RefractionEq {
103            d: 30.0,
104            h: 5.0,
105            w: 1.0,
106            n: 1.33,
107        };
108
109        // Evaluate until a given tolerance.
110        let val = crate::find_root(
111            0.0, 30.0, // Initial guess
112            refraction, 1e-10,
113        )
114        .unwrap();
115        let eps = 0.01;
116        assert!(f64::abs(val - 5.7) < eps);
117    }
118}
refraction/lib.rs

refraction/
lib.rs
