scatrs 0.3.0

Fast, linearized, rigorously validated single-scattering calculations for spherical and spheroidal particles (Lorenz–Mie, layered spheres, T-matrix spheroids, size distributions)
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

//! Values with derivatives and error estimates, propagated to first order through the
//! ratios that turn integrals into ensemble outputs.

/// A value, its error estimate, its derivatives with respect to the requested parameters,
/// and their error estimates.
#[derive(Clone, Debug, PartialEq)]
pub(crate) struct Lin {
    pub(crate) v: f64,
    pub(crate) e: f64,
    pub(crate) d: Vec<f64>,
    pub(crate) de: Vec<f64>,
}

impl Lin {
    /// An exact value without derivatives (`n` parameters).
    pub(crate) fn constant(v: f64, n: usize) -> Self {
        Self {
            v,
            e: 0.0,
            d: vec![0.0; n],
            de: vec![0.0; n],
        }
    }

    /// `self`/`b`: the quotient rule, and the first-order error bounds
    /// e(a/b) ≈ (e_a + |a/b| e_b)/|b| and
    /// e(∂(a/b)) ≈ (e_∂a + |a/b| e_∂b + |∂(a/b)| e_b + e(a/b) |∂b|)/|b|.
    pub(crate) fn div(&self, b: &Self) -> Self {
        let v = self.v / b.v;
        let e = (self.e + v.abs() * b.e) / b.v.abs();
        let d: Vec<f64> = self
            .d
            .iter()
            .zip(&b.d)
            .map(|(da, db)| (da - v * db) / b.v)
            .collect();
        let de = (0..d.len())
            .map(|k| {
                (self.de[k] + v.abs() * b.de[k] + d[k].abs() * b.e + e * b.d[k].abs()) / b.v.abs()
            })
            .collect();
        Self { v, e, d, de }
    }

    /// `c`·`self` for an exact constant `c`.
    pub(crate) fn scale(&self, c: f64) -> Self {
        Self {
            v: c * self.v,
            e: c.abs() * self.e,
            d: self.d.iter().map(|x| c * x).collect(),
            de: self.de.iter().map(|x| c.abs() * x).collect(),
        }
    }
}

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

    #[test]
    fn the_quotient_rule_matches_finite_differences() {
        // a(θ) = 3 + 2θ, b(θ) = 5 − θ at θ = 0.
        let a = Lin {
            v: 3.0,
            e: 1e-3,
            d: vec![2.0],
            de: vec![1e-4],
        };
        let b = Lin {
            v: 5.0,
            e: 2e-3,
            d: vec![-1.0],
            de: vec![0.0],
        };
        let q = a.div(&b);
        let h = 1e-6;
        let fd = ((3.0 + 2.0 * h) / (5.0 - h) - (3.0 - 2.0 * h) / (5.0 + h)) / (2.0 * h);
        assert!((q.d[0] - fd).abs() < 1e-9);
        assert!((q.v - 0.6).abs() < 1e-15);
        // Perturbing a and b by their errors moves a/b by the bound, to first order.
        let worst = ((3.0f64 + 1e-3) / (5.0 - 2e-3) - 0.6).abs();
        assert!(worst <= q.e * (1.0 + 1e-3) && worst >= q.e * (1.0 - 1e-3));
        assert_eq!(a.scale(2.0).d, vec![4.0]);
    }
}