rust_roche/ref_sphere.rs
1use crate::{Vec3, Star};
2use crate::{x_l1_1, x_l1_2, rpot1, rpot2};
3
4
5///
6/// ref_sphere computes the radius of a reference sphere just outside a Roche distorted star
7/// along the line of centres and centred upon its centre of mass. This sphere, which is guaranteed
8/// to enclose the equipotential in question can then be used to define regions for searching for
9/// equipotential crossing when computing eclipses. The Roche-distorted star is defined by the mass
10/// ratio and the (linear) filling factor defined as the distance from the centre of mass of the
11/// star to its surface in the direction of the L1 point divided by the distance to the L1 point.
12/// A filling factor = 1 is Roche filling. Note that the size of the Roche lobe is calculated as
13/// appropriate given any asynchronism.
14///
15/// \param q the mass ratio = M2/M1
16/// \param star specifies which star, primary or secondary is under consideration
17/// \param spin ratio spin/orbital frequencies to allow for asynchronism
18/// \param ffac linear filling factor.
19/// \param rref radius of the reference sphere. This will be 0.1% expanded above the minimum
20/// size to avoid round off bugs, if it remains within Roche lobe.
21/// \param pref reference potential. Roche potential on surface of distorted star.
22///
23pub fn ref_sphere(q: f64, star: Star, spin: f64, ffac: f64) -> (f64, f64) {
24 let tref: f64;
25 let rref: f64;
26 let pref: f64;
27 if star == Star::Primary {
28 tref = x_l1_1(q, spin);
29 rref = tref * 1.0_f64.min(1.001*ffac);
30 pref = rpot1(q, spin, &Vec3 { x: ffac*tref, y: 0.0, z: 0.0 });
31 (rref, pref)
32 } else if star == Star::Secondary {
33 tref = 1.0 - x_l1_2(q, spin);
34 rref = tref * 1.0_f64.min(1.001*ffac);
35 pref = rpot2(q, spin, &Vec3 { x: 1.0 - ffac*tref, y: 0.0, z: 0.0 });
36 (rref, pref)
37 } else {
38 panic!("star is not an instance of Star")
39 }
40}