rust_roche/
face.rs

1use crate::{Vec3, Star};
2use crate::{rpot1, rpot2, drpot1, drpot2};
3
4
5/// 
6/// 'face' computes the position and orientation of a face on either star in a binary assuming Roche geometry given
7/// a direction, a reference radius and a potential.
8///
9/// \param q    the mass ratio = M2/M1.
10/// \param star specifies which star, primary or secondary is under consideration.
11/// \param spin ratio of star in questions spin to the orbital frequency
12/// \param dirn the direction (unit) vector from the centre of mass of the secondary to the face in question.
13/// \param rref reference radius. This is a radius large enough to guarantee crossing of the reference potential. See ref_sphere
14/// \param pref reference potential. This defines the precise location of the face.
15/// \param acc  location accuracy (units of separation)
16/// \param pvec position vector of centre of face (position vector in standard binary coordinates), returned
17/// \param dvec orientation vector perpendicular to face, returned
18/// \param r    distance from centre of mass of star, returned
19/// \param g    magnitude of gravity at face, returned
20/// \exception The routine throws exceptions if it cannot bracket the reference potential. This can occur if the reference radius fails to enclose
21/// 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
22/// 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.
23/// 
24pub fn face(q: f64, star: Star, spin: f64, direction: Vec3, rref: f64, pref: f64, acc: f64) -> (Vec3, Vec3, f64, f64) {
25
26    let mut pvec: Vec3;
27    let mut r: f64;
28
29    let cofm: Vec3 = match star {
30        Star::Primary => Vec3::cofm1(),
31        Star::Secondary => Vec3::cofm2(),
32    };
33
34    let rp: fn(f64, f64, &Vec3) -> f64 = match star {
35        Star::Primary => rpot1,
36        Star::Secondary => rpot2,
37    };
38
39    let drp: fn(f64, f64, &Vec3) -> Vec3 = match star {
40        Star::Primary => drpot1,
41        Star::Secondary => drpot2,
42    };
43
44    let mut tref: f64 = rp(q, spin, &(cofm + rref*direction));
45    if tref < pref {
46        panic!("stuff")
47    }
48
49    let mut r1: f64 = rref/2.;
50    let mut r2: f64 = rref;
51    tref = pref + 1.;
52
53    const MAXSEARCH: i32 = 30;
54    let mut i: i32 = 0;
55    while i < MAXSEARCH && tref > pref {
56        r1 = r2/2.;
57        tref = rp(q, spin, &(cofm + r1*direction));
58        if tref > pref {
59            r2 = r1;
60        }
61        i+=1;
62    }
63    if tref > pref {
64        panic!("other stuff");
65    }
66
67    const MAXCHOP: i32 = 100;
68    let mut nchop: i32 = 0;
69    while r2 - r1 > acc && nchop < MAXCHOP {
70        r = (r1 + r2)/2.;
71        pvec = cofm + r*direction;
72        if rp(q, spin, &pvec) < pref {
73            r1 = r;
74        }else {
75            r2 = r;
76        }
77        nchop += 1;
78    }
79    if nchop == MAXCHOP {
80        panic!("even more stuff");
81    }
82    r = (r1 + r2)/2.;
83    pvec = cofm + r*direction;
84    let mut dvec: Vec3 = drp(q, spin, &pvec);
85    let g = dvec.length();
86    dvec /= g;
87    return (pvec, dvec, r, g)
88}
rust_roche/face.rs

rust_roche/
face.rs
