rust-roche 0.2.3

Rust translation of Tom Marsh's cpp-roche package for modelling Roche distorted stars/binary systems.
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
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

use crate::errors::RocheError;
use crate::{Star, Vec3};
use crate::{fblink, pot_min, ref_sphere, set_earth, sphere_eclipse};
use pyo3::prelude::*;
use std::f64::consts::TAU;

///
/// ingress_egress tests for whether a given point is eclipsed by a Roche-distorted star. If
/// it is, it computes the ingress and egress phases using a binary chop. The accuracy on the
/// phase should be set to be below the expected uncertainties of the phases of your data.
///
/// Arguments:
///
/// * `q`:       the mass ratio = M2/M1.
/// * `star`:    which star, primary or secondary, is doing the eclipsing
/// * `spin`:    ratio of spin to orbit of eclipsing star
/// * `ffac`:    linear filling factor
/// * `iangle`:  inclination angle
/// * `delta`:  the accuracy in phase wanted.
/// * `r`:       position vector of point of interest.
/// * `ingress`: ingress phase (if eclipsed)
/// * `egress`:  egress phase
///
/// Returns:
///
/// * false = not eclipsed; true = eclipsed.
///
pub fn ingress_egress(
    q: f64,
    star: Star,
    spin: f64,
    ffac: f64,
    iangle: f64,
    delta: f64,
    r: &Vec3,
    ingress: &mut f64,
    egress: &mut f64,
) -> Result<bool, RocheError> {
    let rref: f64;
    let pref: f64;
    (rref, pref) = ref_sphere(q, star, spin, ffac)?;
    let ri: f64 = iangle.to_radians();
    let (sini, cosi) = ri.sin_cos();

    let cofm: Vec3 = match star {
        Star::Primary => Vec3::cofm1(),
        Star::Secondary => Vec3::cofm2(),
    };

    let mut phi1: f64 = 0.0;
    let mut phi2: f64 = 0.0;
    let mut lam1: f64 = 0.0;
    let mut lam2: f64 = 0.0;
    let mut phi: f64 = 0.0;
    let mut lam: f64 = 0.0;

    if sphere_eclipse(
        cosi, sini, r, &cofm, rref, &mut phi1, &mut phi2, &mut lam1, &mut lam2,
    ) {
        let acc: f64 = 2. * (2.0 * TAU * (lam2 - lam1) * delta).sqrt();

        if pot_min(
            q, star, spin, cosi, sini, r, phi1, phi2, lam1, lam2, rref, pref, acc, &mut phi,
            &mut lam,
        )? {
            let mut pin: f64 = phi;
            let mut pout: f64 = phi1;
            let mut pmid: f64;

            while (pin - pout).abs() > delta {
                pmid = (pin + pout) / 2.0;
                if fblink(q, star, spin, ffac, acc, &set_earth(cosi, sini, pmid), r).unwrap() {
                    pin = pmid;
                } else {
                    pout = pmid;
                }
            }
            *ingress = (pin + pout) / 2.0;
            *ingress = *ingress - ingress.floor();

            pin = phi;
            pout = phi2;
            while (pin - pout).abs() > delta {
                pmid = (pin + pout) / 2.;
                if fblink(q, star, spin, ffac, acc, &set_earth(cosi, sini, pmid), r).unwrap() {
                    pin = pmid;
                } else {
                    pout = pmid;
                }
            }
            *egress = (pin + pout) / 2.0;
            *egress = *egress - egress.floor();
            if *egress < *ingress {
                *egress += 1.0;
            }
            return Ok(true);
        } else {
            return Ok(false);
        }
    } else {
        return Ok(false);
    }
}

// wrapper for python library, avoiding mutable references

///
/// ingress_egress tests for whether a given point is eclipsed by a Roche-distorted star. If
/// it is, it computes the ingress and egress phases using a binary chop. The accuracy on the
/// phase should be set to be below the expected uncertainties of the phases of your data.
///
/// \param q       the mass ratio = M2/M1.
/// \param star    which star, primary or secondary, is doing the eclipsing
/// \param spin    ratio of spin to orbit of eclipsing star
/// \param ffac    linear filling factor
/// \param iangle  inclination angle
/// \param delta   the accuracy in phase wanted.
/// \param r       position vector of point of interest.
/// \return (eclipsed, ingress_phase, egress_phase)
///
#[pyfunction]
#[pyo3(name = "ingress_egress")]
pub fn ingress_egress_wrapper(
    q: f64,
    star: Star,
    spin: f64,
    ffac: f64,
    iangle: f64,
    delta: f64,
    r: &Vec3,
) -> Result<(bool, f64, f64), RocheError> {
    let mut ingress: f64 = 0.0;
    let mut egress: f64 = 0.0;
    let eclipsed: bool = ingress_egress(
        q,
        star,
        spin,
        ffac,
        iangle,
        delta,
        r,
        &mut ingress,
        &mut egress,
    )?;
    Ok((eclipsed, ingress, egress))
}