roche/
face.rs

1use crate::errors::RocheError;
2use crate::{Star, Vec3};
3use crate::{drpot1, drpot2, rpot1, rpot2};
4use pyo3::prelude::*;
5
6///
7/// 'face' computes the position and orientation of a face on either star in a binary assuming Roche geometry given
8/// a direction, a reference radius and a potential.
9///
10/// Arguments:
11///
12/// * `q`:    the mass ratio = M2/M1.
13/// * `star`: specifies which star, primary or secondary is under consideration.
14/// * `spin`: ratio of star in questions spin to the orbital frequency
15/// * `dirn`: the direction (unit) vector from the centre of mass of the secondary to the face in question.
16/// * `rref`: reference radius. This is a radius large enough to guarantee crossing of the reference potential. See ref_sphere
17/// * `pref`: reference potential. This defines the precise location of the face.
18/// * `acc`:  location accuracy (units of separation)
19///
20/// Returns:
21///
22/// * pvec position vector of centre of face (position vector in standard binary coordinates), returned
23/// * dvec orientation vector perpendicular to face, returned
24/// * r    distance from centre of mass of star, returned
25/// * g    magnitude of gravity at face, returned
26/// \exception The routine throws exceptions if it cannot bracket the reference potential. This can occur if the reference radius fails to enclose
27/// the face in question, or if the face is so deep in the potential that the initial search fails to reach it. Finally if acc is set too low an
28/// exception may be thrown if too many binary chops occur. The behaviour at the L1 point is undefined so do not try to call it there.
29///
30#[pyfunction]
31pub fn face(
32    q: f64,
33    star: Star,
34    spin: f64,
35    direction: Vec3,
36    rref: f64,
37    pref: f64,
38    acc: f64,
39) -> Result<(Vec3, Vec3, f64, f64), RocheError> {
40    let mut pvec: Vec3;
41    let mut r: f64;
42
43    let cofm: Vec3 = match star {
44        Star::Primary => Vec3::cofm1(),
45        Star::Secondary => Vec3::cofm2(),
46    };
47
48    let rp: fn(f64, f64, &Vec3) -> Result<f64, RocheError> = match star {
49        Star::Primary => rpot1,
50        Star::Secondary => rpot2,
51    };
52
53    let drp: fn(f64, f64, &Vec3) -> Result<Vec3, RocheError> = match star {
54        Star::Primary => drpot1,
55        Star::Secondary => drpot2,
56    };
57
58    let mut tref: f64 = rp(q, spin, &(cofm + rref * direction))?;
59    if tref < pref {
60        let message = format!(
61            "point at reference radius {} appears to be at lower potential {} than the reference potential {}",
62            rref, tref, pref
63        );
64        return Err(RocheError::FaceError(message));
65    }
66
67    let mut r1: f64 = rref / 2.;
68    let mut r2: f64 = rref;
69    tref = pref + 1.;
70
71    const MAXSEARCH: i32 = 30;
72    let mut i: i32 = 0;
73    while i < MAXSEARCH && tref > pref {
74        r1 = r2 / 2.;
75        tref = rp(q, spin, &(cofm + r1 * direction))?;
76        if tref > pref {
77            r2 = r1;
78        }
79        i += 1;
80    }
81    if tref > pref {
82        let message = "could not find a radius with a potential below the reference potential; probably bad inputs.";
83        return Err(RocheError::FaceError(message.to_string()));
84    }
85
86    const MAXCHOP: i32 = 100;
87    let mut nchop: i32 = 0;
88    while r2 - r1 > acc && nchop < MAXCHOP {
89        r = (r1 + r2) / 2.;
90        pvec = cofm + r * direction;
91        if rp(q, spin, &pvec)? < pref {
92            r1 = r;
93        } else {
94            r2 = r;
95        }
96        nchop += 1;
97    }
98    if nchop == MAXCHOP {
99        return Err(RocheError::FaceError(
100            "reached maximum number of binary chops".to_string(),
101        ));
102    }
103    r = (r1 + r2) / 2.;
104    pvec = cofm + r * direction;
105    let mut dvec: Vec3 = drp(q, spin, &pvec)?;
106    let g = dvec.length();
107    dvec /= g;
108    Ok((pvec, dvec, r, g))
109}
110
111#[cfg(test)]
112mod tests {
113    use crate::ref_sphere::ref_sphere;
114
115    use super::*;
116
117    #[test]
118    fn face_test() -> Result<(), RocheError> {
119        // Values from trm.roche.face
120        let dirn = Vec3::new(1.0, 1.0, 1.0).norm();
121        let (rref, pref) = ref_sphere(0.2, Star::Secondary, 1.0, 0.8)?;
122        let (pvec, dvec, r, g) = face(0.2, Star::Secondary, 1.0, dirn, rref, pref, 1.0e-5)?;
123        assert_eq!(
124            pvec,
125            Vec3::new(1.1352429036726859, 0.13524290367268593, 0.13524290367268593)
126        );
127        assert_eq!(
128            dvec,
129            Vec3::new(0.49134297246922065, 0.5917653290265037, 0.6390585878988438)
130        );
131        assert_eq!(r, 0.2342475805242355);
132        assert_eq!(g, 2.8596655611690993);
133        Ok(())
134    }
135}
roche/face.rs

roche/
face.rs
