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

extern crate std;
extern crate itertools;
use std::f64::consts;
use std::sync::Arc;
use self::itertools::cons_tuples;

use star::Star;
use boundingshape::BoundingShape;
use linspace::floatrange;

/// A circular starspot
pub struct Spot {
    pub star: Arc<Star>,
    pub latitude: f64,
    pub longitude: f64,
    pub radius: f64,
    pub temperature: f64,
    pub plage: bool,
    pub mortal: bool,
    pub time_appear: f64,
    pub time_disappear: f64,
    pub intensity: f64,
}

impl Spot {
    /// Create a new spot on `star` and at the coordinates specified, where
    /// latitude is 0 at the equator and both latitude and longitude are
    /// measured in degrees.
    pub fn new(
        star: Arc<Star>,
        latitude: f64,
        longitude: f64,
        fillfactor: f64,
        plage: bool,
        mortal: bool,
    ) -> Self {
        let temperature = star.temperature - star.spot_temp_diff;
        Spot {
            star: star,
            latitude: latitude * consts::PI / 180.0,
            longitude: longitude * consts::PI / 180.0,
            radius: (2.0 * fillfactor).sqrt(),
            temperature: temperature,
            plage: plage,
            mortal: mortal,
            time_appear: 0.0,
            time_disappear: 15.0,
            intensity: 0.0,
        }
    }

    pub fn get_flux(&self, time: f64) -> f64 {
        let bounds = BoundingShape::new(self, time);
        if let Some(y_bounds) = bounds.y_bounds() {
            let limb_integral: f64 = floatrange(
                y_bounds.lower,
                y_bounds.upper,
                2.0 / self.star.grid_size as f64,
            ).map(|y| {
                // TODO I should be able to filter_map here
                if let Some(z_bounds) = bounds.z_bounds(y) {
                    self.star.limb_integral(&z_bounds, y)
                } else {
                    0.0
                }
            })
                .sum();
            (1.0 - self.intensity) * limb_integral
        } else {
            0.0
        }
    }

    pub fn get_ccf(&self, time: f64) -> Vec<f64> {
        let mut profile = vec![0.0; self.star.profile_active.len()];
        let bounds = BoundingShape::new(self, time);
        if let Some(y_bounds) = bounds.y_bounds() {
            for y in floatrange(
                y_bounds.lower,
                y_bounds.upper,
                2.0 / self.star.grid_size as f64,
            )
            {
                let quiet_shifted = self.star.profile_quiet.shift(
                    y * self.star.equatorial_velocity,
                );
                let active_shifted = self.star.profile_active.shift(
                    y * self.star.equatorial_velocity,
                );

                if let Some(z_bounds) = bounds.z_bounds(y) {
                    let limb_integral = self.star.limb_integral(&z_bounds, y);
                    for (tot, qshift, ashift) in
                        cons_tuples(profile.iter_mut().zip(quiet_shifted.iter()).zip(
                            active_shifted.iter(),
                        ))
                    {
                        *tot += (qshift - self.intensity * ashift) * limb_integral;
                    }
                }
            }
        };
        profile
    }

    /// Returns `true` if this spot currently exists.
    pub fn alive(&self, time: f64) -> bool {
        if !self.mortal {
            true
        } else {
            time >= self.time_appear && time <= self.time_disappear
        }
    }

    /// Returns whether a spot _ever_ collides with the `other`.
    pub fn collides_with(&self, other: &Spot) -> bool {
        let bounds = BoundingShape::new(self, 0.0);
        let other_bounds = BoundingShape::new(other, 0.0);
        bounds.collides_with(&other_bounds)
    }
}