roche/
ref_sphere.rs

1use crate::errors::RocheError;
2use crate::{Star, Vec3};
3use crate::{rpot1, rpot2, x_l1_1, x_l1_2};
4use pyo3::prelude::*;
5
6///
7/// ref_sphere computes the radius of a reference sphere just outside a Roche distorted star
8/// along the line of centres and centred upon its centre of mass. This sphere, which is guaranteed
9/// to enclose the equipotential in question can then be used to define regions for searching for
10/// equipotential crossing when computing eclipses. The Roche-distorted star is defined by the mass
11/// ratio and the (linear) filling factor defined as the distance from the centre of mass of the
12/// star to its surface in the direction of the L1 point divided by the distance to the L1 point.
13/// A filling factor = 1 is Roche filling. Note that the size of the Roche lobe is calculated as
14/// appropriate given any asynchronism.
15///
16/// Arguments:
17///
18/// * `q`: the mass ratio = M2/M1
19/// * `star`: specifies which star, primary or secondary is under consideration
20/// * `spin`: ratio spin/orbital frequencies to allow for asynchronism
21/// * `ffac`: linear filling factor.
22///
23/// Returns:
24///
25/// * the radius of the reference sphere. This will be 0.1% expanded above the minimum
26/// size to avoid round off bugs, if it remains within Roche lobe.
27/// * the reference potential. Roche potential on surface of distorted star.
28///
29#[pyfunction]
30pub fn ref_sphere(q: f64, star: Star, spin: f64, ffac: f64) -> Result<(f64, f64), RocheError> {
31    let tref: f64;
32    let rref: f64;
33    let pref: f64;
34
35    if star == Star::Primary {
36        tref = x_l1_1(q, spin)?;
37        rref = tref * 1.0_f64.min(1.001 * ffac);
38        pref = rpot1(
39            q,
40            spin,
41            &Vec3 {
42                x: ffac * tref,
43                y: 0.0,
44                z: 0.0,
45            },
46        )?;
47        Ok((rref, pref))
48    } else if star == Star::Secondary {
49        tref = 1.0 - x_l1_2(q, spin)?;
50        rref = tref * 1.0_f64.min(1.001 * ffac);
51        pref = rpot2(
52            q,
53            spin,
54            &Vec3 {
55                x: 1.0 - ffac * tref,
56                y: 0.0,
57                z: 0.0,
58            },
59        )?;
60        Ok((rref, pref))
61    } else {
62        let message = format!("{:?} is not and instance of Star.", star);
63        Err(RocheError::ParameterError(message))
64    }
65}
66
67#[cfg(test)]
68mod tests {
69    use super::*;
70
71    #[test]
72    fn ref_sphere_test() -> Result<(), RocheError> {
73        // Values from trm.roche.ref_sphere
74        assert_eq!(
75            ref_sphere(0.2, Star::Secondary, 1.0, 0.8)?,
76            (0.27342861229381593, -2.3996722470168605)
77        );
78        assert!(ref_sphere(-0.2, Star::Secondary, 1.0, 0.8).is_err());
79        Ok(())
80    }
81}
roche/ref_sphere.rs

roche/
ref_sphere.rs
