lather 0.0.0

A Rust crate and Python extension for modeling starspot effects on radial velocity and photometric observations, inspired by complaints about the SOAP project of similar goals.
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

extern crate std;
extern crate rulinalg;
use self::rulinalg::matrix::Matrix;
use self::rulinalg::matrix::BaseMatrix;

pub fn fit_rv(rv: &[f64], ccf: &[f64]) -> f64 {
    let (min_index, _) = ccf.iter().enumerate().fold(
        (0, std::f64::INFINITY),
        |(min_ind, min_val),
         (current_ind, current_val)| {
            if *current_val < min_val {
                (current_ind, *current_val)
            } else {
                (min_ind, min_val)
            }
        },
    );

    let peak_rv = rv.iter().skip(min_index - 3).take(7).cloned();
    let peak_ccf: Vec<f64> = ccf.iter().skip(min_index - 3).take(7).cloned().collect();
    let rv_matrix_values: Vec<f64> = std::iter::repeat(1.0)
        .take(7)
        .chain(peak_rv.clone())
        .chain(peak_rv.map(|x| x.powi(2)))
        .collect();

    let rv_matrix_transpose = Matrix::new(3, 7, rv_matrix_values.clone());
    let rv_matrix = rv_matrix_transpose.transpose();
    let ccf_matrix = Matrix::new(7, 1, peak_ccf);
    let coefficients = (rv_matrix_transpose.clone() * rv_matrix.clone())
        .inverse()
        .unwrap() * rv_matrix_transpose * ccf_matrix;

    -coefficients[[1, 0]] / (2.0 * coefficients[[2, 0]])
}

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

    fn is_close(actual: f64, expected: f64) -> bool {
        let precision = 5;
        let pow = 10.0f64.powi(precision + 1);
        let delta = (expected - actual).abs();
        let max_delta = 10.0f64.powi(-precision) / 2.0;
        return (delta * pow).round() / pow <= max_delta;
    }

    #[test]
    fn simple_gaussian() {
        use linspace::linspace;

        let test_len = 101;
        let rv: Vec<f64> = linspace(-1.0, 1.0, test_len).collect();
        let ccf: Vec<f64> = rv.iter().map(|x| -(-(x - 0.5).powi(2)).exp()).collect();

        let poly_rv = fit_rv(&rv, &ccf);
        assert!(is_close(poly_rv, 0.5));
    }
}