scicalc_rs/
measurement.rs

1use float_cmp::{ApproxEq, F64Margin};
2use std::fmt;
3use std::ops::{Add, Div, Mul, Neg, Sub};
4
5/**A Measurement 'x' is written as x = (mean +- sigma)
6
7where 'mean' is the mean value
8
9and 'sigma' is the uncertainty(also called error or standard deviation from the mean)*/
10#[derive(Debug, Clone, Copy)]
11pub struct Measurement {
12    pub mean: f64,  //mean value
13    pub sigma: f64, //std deviation, error or uncertainty
14}
15
16impl Measurement {
17    pub fn new(mean: f64, sigma: f64) -> Measurement {
18        Measurement { mean, sigma }
19    }
20}
21
22impl Neg for Measurement {
23    type Output = Self;
24
25    fn neg(self) -> Self {
26        Self {
27            mean: -self.mean,
28            sigma: self.sigma,
29        }
30    }
31}
32
33impl Add for Measurement {
34    type Output = Self;
35
36    fn add(self, other: Self) -> Self {
37        Self {
38            mean: self.mean + other.mean,
39            sigma: quadrature(self.sigma, other.sigma).sqrt(),
40        }
41    }
42}
43
44impl Add<f64> for Measurement {
45    type Output = Self;
46
47    fn add(self, other: f64) -> Self {
48        Self {
49            mean: self.mean + other,
50            sigma: self.sigma,
51        }
52    }
53}
54
55impl Sub for Measurement {
56    type Output = Self;
57
58    fn sub(self, other: Self) -> Self {
59        Self {
60            mean: self.mean - other.mean,
61            sigma: quadrature(self.sigma, other.sigma).sqrt(),
62        }
63    }
64}
65
66impl Sub<f64> for Measurement {
67    type Output = Self;
68
69    fn sub(self, other: f64) -> Self {
70        Self {
71            mean: self.mean - other,
72            sigma: self.sigma,
73        }
74    }
75}
76
77impl Mul for Measurement {
78    type Output = Self;
79
80    fn mul(self, other: Self) -> Self {
81        let new_mean = self.mean * other.mean;
82        Self {
83            mean: new_mean,
84            sigma: quadrature(self.sigma / self.mean, other.sigma / other.mean).sqrt() * new_mean,
85        }
86    }
87}
88
89impl Mul<f64> for Measurement {
90    type Output = Self;
91
92    fn mul(self, other: f64) -> Self {
93        Self {
94            mean: self.mean * other,
95            sigma: self.sigma * other,
96        }
97    }
98}
99
100impl Div for Measurement {
101    type Output = Self;
102
103    fn div(self, other: Self) -> Self {
104        let new_mean = self.mean / other.mean;
105        Self {
106            mean: new_mean,
107            sigma: quadrature(self.sigma / self.mean, other.sigma / other.mean).sqrt() * new_mean,
108        }
109    }
110}
111
112impl Div<f64> for Measurement {
113    type Output = Self;
114
115    fn div(self, other: f64) -> Self {
116        Self {
117            mean: self.mean / other,
118            sigma: self.sigma / other,
119        }
120    }
121}
122
123impl fmt::Display for Measurement {
124    fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result {
125        write!(f, "{} ± {}", self.mean, self.sigma)
126    }
127}
128
129impl PartialEq for Measurement {
130    fn eq(&self, other: &Self) -> bool {
131        self.mean == other.mean && self.sigma == other.sigma
132    }
133}
134
135impl ApproxEq for Measurement {
136    type Margin = F64Margin;
137
138    fn approx_eq<T: Into<Self::Margin>>(self, other: Self, margin: T) -> bool {
139        let margin = margin.into();
140        self.mean.approx_eq(other.mean, margin) && self.sigma.approx_eq(other.sigma, margin)
141    }
142}
143
144impl From<f64> for Measurement {
145    fn from(x: f64) -> Self {
146        Measurement {
147            mean: x,
148            sigma: 0.0,
149        }
150    }
151}
152
153pub fn pow(x: Measurement, y: Measurement) -> Measurement {
154    //x^y
155    let [xm, ym, sx, sy] = [x.mean, y.mean, x.sigma, y.sigma];
156    let dfdx = ym * xm.powf(ym - 1.0);
157    let dfdy = xm.ln() * xm.powf(ym);
158
159    Measurement {
160        mean: xm.powf(ym),
161        sigma: quadrature(dfdx * sx, dfdy * sy).sqrt(),
162    }
163}
164
165fn quadrature(x: f64, y: f64) -> f64 {
166    x * x + y * y
167}
168
169#[cfg(test)]
170mod tests {
171    use super::*;
172    #[test]
173    fn simple_operations() {
174        let x = Measurement::new(1.0, 0.01);
175        let y = Measurement::new(2.0, 0.01);
176
177        let added = Measurement::new(3.0, 0.0002f64.sqrt());
178        let subtracted = Measurement::new(-1.0, 0.0002f64.sqrt());
179        let multiplied = Measurement::new(2.0, 2.0 * (quadrature(0.01 / 1.0, 0.01 / 2.0)).sqrt());
180        let divided = Measurement::new(0.5, 0.5 * (quadrature(0.01 / 1.0, 0.01 / 2.0)).sqrt());
181
182        assert!(added.approx_eq(x + y, F64Margin::default()));
183        assert!(subtracted.approx_eq(x - y, F64Margin::default()));
184        assert!(multiplied.approx_eq(x * y, F64Margin::default()));
185        assert!(divided.approx_eq(x / y, F64Margin::default()));
186    }
187    #[test]
188    fn approximate_equality() {
189        /* Tests the approximate equality due to floating point errors */
190        let x = Measurement::new(1.0, 0.01);
191        let y = Measurement::new(1.0, 0.001);
192        let x_prime = Measurement::new(1.0, 2.0 * 0.005);
193
194        //A number is equal to itself. Therefore it's also approximately equal to itself
195        assert!(x.approx_eq(x, F64Margin::default()));
196        //x and y should NOT be approximately equal
197        assert_eq!(false, x.approx_eq(y, F64Margin::default()));
198        //x and x_prime should be equal
199        assert!(x.approx_eq(x_prime, F64Margin::default()));
200    }
201    #[test]
202    fn simple_pow() {
203        let x = Measurement::from(2.0);
204        let y = Measurement::from(3.0);
205        let r = Measurement {
206            mean: 8.0,
207            sigma: 0.0,
208        };
209        assert_eq!(pow(x, y), r);
210    }
211
212    #[test]
213    fn negative_pow() {
214        let x = Measurement::from(2.0);
215        let y = Measurement::from(-3.0);
216        let r = Measurement {
217            mean: 1.0 / 8.0,
218            sigma: 0.0,
219        };
220        assert_eq!(pow(x, y), r);
221    }
222}
scicalc_rs/measurement.rs

scicalc_rs/
measurement.rs
