rust_roche/
sphere_eclipse.rs

1use std::f64::consts::TAU;
2use crate::Vec3;
3
4
5///
6/// sphere_eclipse tells you whether or not a given sphere will eclipse a 
7/// given point or not. If the answer is yes, it will also return with four 
8/// parameters to define the phase range and the multiplier range delimiting the
9/// region within which the spheres surface is crossed.
10/// These can then be used as the starting point for later computation.
11/// (The line of sight is described as the point in question plus a scalar multiplier
12/// times a unit vector pointing towards Earth -- this is the "multiplier" referred to above
13/// and below). The multiplier must be positive: in other words the routine does not
14/// project backwards. If the point in inside the sphere, phi1 will be set = 0, phi2 = 1,
15/// lam1 = 0, and lam2 = the largest value of the multiplier lambda
16///
17/// \param cosi cosine of orbital inclination
18/// \param sini sine of orbital inclination
19/// \param r       the position vector of the point in question (units of binary separation)
20/// \param c       the centre of the sphere enclosing the star (units of binary separation)
21/// \param rsphere the radius defining the sphere enclosing the star (units of binary separation)
22/// \param phi1    if eclipse by the sphere, this is the start of the phase range (0 -1)
23/// \param phi2    if eclipse by the sphere, this is the end of the phase range (least amount > phi1)
24/// \param lam1    if eclipse by the sphere, this is the start of the multiplier range (>=0)
25/// \param lam2    if eclipse by the sphere, this is the end of the multiplier range (>lam1)
26/// \return false = not eclipsed; true = eclipsed.
27///
28pub fn sphere_eclipse(cosi: f64, sini: f64, r: &Vec3, c: &Vec3, rsphere: f64, phi1: &mut f64, phi2: &mut f64, lam1: &mut f64, lam2: &mut f64) -> bool {
29
30    let d: Vec3 = *r - *c;
31
32    let pdist: f64 = (d.x*d.x + d.y*d.y).sqrt();
33    let bquad: f64 = d.z*cosi - pdist*sini;
34    if bquad >= 0. {
35        return false;
36    }
37    let cquad: f64 = d.sqr() - rsphere*rsphere;
38
39    let mut fac: f64 = bquad*bquad - cquad;
40    if fac <= 0.0 {
41        return false
42    }
43    fac = fac.sqrt();
44
45    *lam2 = -bquad + fac;
46    *lam1 = 0_f64.max(cquad/(*lam2));
47
48    if cquad < 0. {
49        *phi1 = 0.;
50        *phi2 = 1.
51    }else {
52        let delta: f64 = ((cosi*d.z + cquad.sqrt())/(sini*pdist)).acos();
53        let phi: f64 = d.y.atan2(-d.x);
54        *phi1 = (phi - delta)/TAU;
55        *phi1 -= phi1.floor();
56        *phi2 = *phi1 + 2.*delta/TAU;
57    }
58    return true;
59
60}
61
62
63///
64/// This version of sphere_eclipse tells you whether or not a given sphere will eclipse 
65/// a given point at a particular phase or not. If the answer is yes,
66/// it will also return with the multiplier values giving the cut points. These can then
67/// be used as starting points for Roche lobe computations. These can then be used as the 
68/// starting point for later computation. Points inside the sphere are regarded as being
69/// eclipsed with the lower mulitplier set = 0
70///
71/// \param earth   vector towards Earth
72/// \param r       the position vector of the point in question (units of binary separation)
73/// \param c       the centre of the sphere (units of binary separation)
74/// \param rsphere the radius defining the sphere enclosing the star (units of binary separation)
75/// \param lam1    if eclipse by the sphere, this is the start of the multiplier range
76/// \param lam2    if eclipse by the sphere, this is the end of the multiplier range
77/// \return false = not eclipsed; true = eclipsed.
78///
79pub fn sphere_eclipse_vector(earth: &Vec3, r: &Vec3, c: &Vec3, rsphere: f64, lam1: &mut f64, lam2: &mut f64) -> bool {
80
81    let d: Vec3 = *r - *c;
82
83    let bquad: f64 = earth.dot(&d);
84    if bquad >= 0. {
85        return false;
86    }
87    let cquad: f64 = d.sqr() - rsphere*rsphere;
88
89    let mut fac: f64 = bquad*bquad - cquad;
90    if fac <= 0. {
91        return false
92    }
93    fac = fac.sqrt();
94
95    *lam2 = -bquad + fac;
96    *lam1 = 0_f64.max(cquad/(*lam2));
97    return true;
98}
rust_roche/sphere_eclipse.rs

rust_roche/
sphere_eclipse.rs
