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
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

extern crate std;
use std::iter;

/// A cross-correlation profile, which can be shifted by fast linear interpolation.
pub struct Profile {
    pub rv: Vec<f64>,
    pub ccf: Vec<f64>,
    derivative: Vec<f64>,
    stepsize: f64,
}

impl Profile {
    /// Creates a profile with the provided radial velocity and
    /// cross-correlation function, and computes the derivative to enable fast
    /// linear interpolation to simulate viewing the profile at different
    /// relative velocities.
    pub fn new(rv: Vec<f64>, ccf: Vec<f64>) -> Self {
        let ccf_diff = ccf.windows(2).map(|s| s[0] - s[1]);
        let rv_diff = rv.windows(2).map(|s| s[0] - s[1]);

        let der = ccf_diff
            .zip(rv_diff)
            .map(|(c, r)| c / r)
            .chain(std::iter::once(0.0))
            .collect();

        Profile {
            rv: rv.clone(),
            ccf: ccf.clone(),
            derivative: der,
            stepsize: (rv[0] - rv[1]).abs(),
        }
    }

    /// Returns the number of elements in the profile.
    pub fn len(&self) -> usize {
        self.rv.len()
    }

    /// Returns `true` if the profile has a length of 0.
    pub fn is_empty(&self) -> bool {
        self.len() == 0
    }

    /// Uses the pre-computed derivative to compute a shifted version of
    /// this profile's cross-correlation function by linear interpolation.
    /// The units of velocity must match those of the radial velocity used
    /// to construct this profile.
    pub fn shift(&self, velocity: f64) -> Vec<f64> {
        let quotient = (velocity / self.stepsize).round() as isize;
        let remainder = velocity - (quotient as f64) * self.stepsize;

        if velocity >= 0.0 {
            iter::repeat(self.ccf[0])
                .take(quotient as usize)
                .chain(
                    self.ccf
                        .iter()
                        .zip(self.derivative.iter())
                        .take(self.ccf.len() - quotient as usize)
                        .map(|(ccf, der)| ccf - remainder * der),
                )
                .collect()
        } else {
            self.ccf
                .iter()
                .zip(self.derivative.iter())
                .skip((-quotient) as usize)
                .map(|(ccf, der)| ccf - remainder * der)
                .chain(iter::repeat(self.ccf[0]).take((-quotient) as usize))
                .collect()
        }
    }
}

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

    #[test]
    fn derivative_len() {
        let test_profile = Profile::new(rv(), ccf_quiet());
        assert_eq!(
            test_profile.derivative.len(),
            ccf_quiet().len(),
            "derivative length is {} but should be {}",
            test_profile.derivative.len(),
            ccf_quiet().len()
        );
    }

    #[test]
    fn pos_zero_shift_len() {
        let test_profile = Profile::new(rv(), ccf_quiet());
        let shifted = test_profile.shift(0.0);
        assert_eq!(
            shifted.len(),
            ccf_quiet().len(),
            "zero-shifted length is {} but should be {}",
            shifted.len(),
            ccf_quiet().len()
        );
    }

    #[test]
    fn neg_zero_shift_len() {
        let test_profile = Profile::new(rv(), ccf_quiet());
        let shifted = test_profile.shift(0.0);
        assert_eq!(
            shifted.len(),
            ccf_quiet().len(),
            "negative zero-shifted length is {} but should be {}",
            shifted.len(),
            ccf_quiet().len()
        );
    }

    #[test]
    fn pos_shift_len() {
        let test_profile = Profile::new(rv(), ccf_quiet());
        let shifted = test_profile.shift(1.5);
        assert_eq!(
            shifted.len(),
            ccf_quiet().len(),
            "postive-shifted length is {} but should be {}",
            shifted.len(),
            ccf_quiet().len()
        );
    }

    #[test]
    fn neg_shift_len() {
        let test_profile = Profile::new(rv(), ccf_quiet());
        let shifted = test_profile.shift(-1.5);
        assert_eq!(
            shifted.len(),
            ccf_quiet().len(),
            "negative-shifted length is {} but should be {}",
            shifted.len(),
            ccf_quiet().len()
        );
    }

    #[test]
    fn pos_zero_shift() {
        let test_profile = Profile::new(rv(), ccf_quiet());
        let shifted = test_profile.shift(0.0);
        for (original, shifted) in ccf_quiet().iter().zip(shifted.iter()) {
            assert_eq!(
                original,
                shifted,
                "zero-shifted profile value is {} but should be {}",
                original,
                shifted
            );
        }
    }

    #[test]
    fn neg_zero_shift() {
        let test_profile = Profile::new(rv(), ccf_quiet());
        let shifted = test_profile.shift(-0.0);
        for (original, shifted) in ccf_quiet().iter().zip(shifted.iter()) {
            assert_eq!(
                original,
                shifted,
                "negative zero-shifted profile value is {} but should be {}",
                original,
                shifted
            );
        }
    }
}