roche/
roche_context.rs

1use std::f64::consts::TAU;
2use crate::errors::RocheError;
3use crate::{Vec3, Star, Etype};
4use crate::{rpot_val, rpot_grad, rpot1, rpot2, drpot1, drpot2};
5use crate::{pot_min, dbrent, x_l1, x_l1_1, x_l1_2, x_l2, x_l3};
6use crate::{sphere_eclipse, sphere_eclipse_vector, set_earth};
7
8
9///
10/// RocheContext is a basic struct packaging the mass ratio, q=M2/M1,
11/// the star (i.e. the Primary or Secondary), the spin of that star,
12/// and the x-coordinate of the L1 point together with the functions
13/// that accept these values as inputs. x_l1 is calculated on initialisation
14/// with RocheContext::new(q, star, spin) and saves it being recalculated
15/// over and over in loops.
16/// 
17pub struct RocheContext {
18    pub q: f64,
19    pub star: Star,
20    pub spin: f64,
21    pub x_l1: f64,
22}
23
24impl RocheContext {
25    pub fn new(q: f64, star: Star, spin: f64) -> Result<Self, RocheError> {
26
27        let x_l1: f64 = match star {
28            Star::Primary => x_l1_1(q, spin)?,
29            Star::Secondary => x_l1_2(q, spin)?,
30        };
31        Ok(Self { q, star, spin, x_l1 })
32    }
33
34    pub fn potential(&self, earth: &Vec3, p:&Vec3, lam: f64) -> Result<f64, RocheError> {
35        Ok(rpot_val(self.q, self.star, self.spin, earth, p, lam)?)
36    }
37
38
39    pub fn gradient(&self, earth: &Vec3, p: &Vec3, lam: f64) -> Result<(f64, f64), RocheError> {
40        let (dp, dl) = rpot_grad(self.q, self.star, self.spin, earth, p, lam)?;
41        Ok((dp, dl))
42    }
43
44
45    pub fn potential_grad(&self, earth: &Vec3, p: &Vec3, lam: f64) -> Result<(f64, f64, f64), RocheError> {
46        let f = self.potential(earth, p, lam)?;
47        let (dp, dl) = rpot_grad(self.q, self.star, self.spin, earth, p, lam)?;
48        Ok((f, dp, dl))
49    }
50
51
52    pub fn ref_sphere(&self, ffac: f64) -> Result<(f64, f64), RocheError> {
53        let tref: f64;
54        let rref: f64;
55        let pref: f64;
56        if self.star == Star::Primary {
57            tref = self.x_l1;
58            rref = tref * 1.0_f64.min(1.001*ffac);
59            pref = rpot1(self.q, self.spin, &Vec3 { x: ffac*tref, y: 0.0, z: 0.0 })?;
60            Ok((rref, pref))
61        } else if self.star == Star::Secondary {
62            tref = 1.0 - self.x_l1;
63            rref = tref * 1.0_f64.min(1.001*ffac);
64            pref = rpot2(self.q, self.spin, &Vec3 { x: 1.0 - ffac*tref, y: 0.0, z: 0.0 })?;
65            Ok((rref, pref))
66        } else {
67            let message = format!("{:?} is not and instance of Star.", self.star);
68            return Err(RocheError::ParameterError(message));
69        }
70    }
71
72
73    pub fn fblink(&self, ffac: f64, acc: f64, earth: &Vec3, p: &Vec3) -> Result<bool, RocheError> {
74
75        let (rref, pref) = self.ref_sphere(ffac)?;
76
77        let cofm: Vec3 = match self.star {
78            Star::Primary => Vec3::cofm1(),
79            Star::Secondary => Vec3::cofm2(),
80        };
81        
82        // First compute the multipliers cutting the reference sphere (if any)
83        let mut lam1 = 0.0;
84        let mut lam2 = 0.0;
85        if !sphere_eclipse_vector(earth, p, &cofm, rref, &mut lam1, &mut lam2) {
86            return Ok(false);
87        }
88        if lam1 == 0.0 {
89            return Ok(true);
90        }
91
92        // Create function objects for 1D minimisation in lambda direction
93        let func = |lam: f64| {
94            Ok(self.potential(earth, p, lam)?)
95        };
96
97        // Now try to bracket a minimum. We just crudely compute function at regularly spaced intervals filling in the
98        // gaps until the step size between the points drops below the threshold. Take every opportunity to jump out early
99        // either if the potential is below the threshold or if we have bracketed a minimum.
100        let mut nstep: i32 = 1;
101        let mut step: f64 = lam2 - lam1;
102
103        let mut f1: f64 = 0.0;
104        let mut f2: f64 = 0.0;
105        let mut flam: f64 = 1.0;
106        let mut lam: f64 = lam1;
107
108        while step > acc {
109
110            lam = lam1 + step/2.0;
111
112            for _ in 0..nstep {
113
114                flam = func(lam)?;
115                if flam <= pref {
116                    return Ok(true);
117                }
118
119                // Calculate these as late as possible because they may often not be needed
120                if nstep == 1 {
121                    f1 = func(lam1)?;
122                    f2 = func(lam2)?;
123                }
124
125                if flam < f1 && flam < f2 {
126                    break;
127                }
128
129                lam += step;
130            }
131            if flam < f1 && flam < f2 {
132                break;
133            }
134            step /= 2.0;
135            nstep *= 2;
136        }
137
138        if flam < f1 && flam < f2 {
139
140            // OK, minimum bracketted, so finally pin it down accurately
141            // Possible that multiple minima could cause problems but I have
142            // never seen this in practice.
143            let dfunc = |lam: f64| {
144                let (_dp, dl) = self.gradient(earth, p, lam)?;
145                Ok(dl)
146            };
147
148            let (_xmin, flam) = dbrent(lam1, lam, lam2, |x| func(x), |x| dfunc(x), acc, true, pref)?;
149
150            Ok(flam < pref)
151        } else {
152            // Not bracketted even after a detailed search, and we have not jumped 
153	        // out either, so assume no eclipse
154            Ok(false)
155        }
156
157
158    }
159
160
161    pub fn face(&self, direction: Vec3, rref: f64, pref: f64, acc: f64) -> Result<(Vec3, Vec3, f64, f64), RocheError> {
162
163        let mut pvec: Vec3;
164        let mut r: f64;
165
166        let cofm: Vec3 = match self.star {
167            Star::Primary => Vec3::cofm1(),
168            Star::Secondary => Vec3::cofm2(),
169        };
170
171        let rp: fn(f64, f64, &Vec3) -> Result<f64, RocheError> = match self.star {
172            Star::Primary => rpot1,
173            Star::Secondary => rpot2,
174        };
175
176        let drp: fn(f64, f64, &Vec3) -> Result<Vec3, RocheError> = match self.star {
177            Star::Primary => drpot1,
178            Star::Secondary => drpot2,
179        };
180
181        let mut tref: f64 = rp(self.q, self.spin, &(cofm + rref*direction))?;
182        if tref < pref {
183            let message = format!(
184            "point at reference radius {} appears to be at lower potential {} than the reference potential {}", rref, tref, pref
185            );
186            return Err(RocheError::FaceError(message));
187        }
188
189        let mut r1: f64 = rref/2.;
190        let mut r2: f64 = rref;
191        tref = pref + 1.;
192
193        const MAXSEARCH: i32 = 30;
194        let mut i: i32 = 0;
195        while i < MAXSEARCH && tref > pref {
196            r1 = r2/2.;
197            tref = rp(self.q, self.spin, &(cofm + r1*direction))?;
198            if tref > pref {
199                r2 = r1;
200            }
201            i+=1;
202        }
203        if tref > pref {
204            let message = "could not find a radius with a potential below the reference potential; probably bad inputs.";
205            return Err(RocheError::FaceError(message.to_string()));
206        }
207
208        const MAXCHOP: i32 = 100;
209        let mut nchop: i32 = 0;
210        while r2 - r1 > acc && nchop < MAXCHOP {
211            r = (r1 + r2)/2.;
212            pvec = cofm + r*direction;
213            if rp(self.q, self.spin, &pvec)? < pref {
214                r1 = r;
215            }else {
216                r2 = r;
217            }
218            nchop += 1;
219        }
220        if nchop == MAXCHOP {
221            return Err(RocheError::FaceError("reached maximum number of binary chops".to_string()));
222        }
223        r = (r1 + r2)/2.;
224        pvec = cofm + r*direction;
225        let mut dvec: Vec3 = drp(self.q, self.spin, &pvec)?;
226        let g = dvec.length();
227        dvec /= g;
228        return Ok((pvec, dvec, r, g))
229    }
230
231
232    pub fn ingress_egress(&self, ffac: f64, iangle: f64, delta: f64, r: &Vec3, ingress: &mut f64, egress: &mut f64) -> Result<bool, RocheError> {
233        let rref: f64;
234        let pref: f64;
235        (rref, pref) = self.ref_sphere(ffac)?;
236        let ri: f64 = iangle.to_radians();
237        let (sini, cosi) = ri.sin_cos();
238
239        let cofm: Vec3 = match self.star {
240            Star::Primary => Vec3::cofm1(),
241            Star::Secondary => Vec3::cofm2(),
242        };
243
244        let mut phi1: f64 = 0.0;
245        let mut phi2: f64 = 0.0;
246        let mut lam1: f64 = 0.0;
247        let mut lam2: f64 = 0.0;
248        let mut phi: f64 = 0.0;
249        let mut lam: f64 = 0.0;
250
251        if sphere_eclipse(cosi, sini, r, &cofm, rref, &mut phi1, &mut phi2, &mut lam1, &mut lam2) {
252            
253            let acc: f64 = 2.*(2.*TAU*(lam2 - lam1)*delta).sqrt();
254
255            if self.pot_min(cosi, sini, r, phi1, phi2, lam1, lam2, rref, pref, acc, &mut phi, &mut lam)? {
256
257                let mut pin: f64 = phi;
258                let mut pout: f64 = phi1;
259                let mut pmid: f64;
260
261                while (pin - pout).abs() > delta {
262                    pmid = (pin + pout)/2.0;
263                    if self.fblink(ffac, acc, &set_earth(cosi, sini, pmid), r).unwrap() {
264                        pin = pmid;
265                    } else {
266                        pout = pmid;
267                    }
268                }
269                *ingress = (pin+pout)/2.0;
270                *ingress = *ingress - ingress.floor();
271
272                pin = phi;
273                pout = phi2;
274                while (pin-pout).abs() > delta {
275                    pmid = (pin+pout)/2.;
276                    if self.fblink(ffac, acc, &set_earth(cosi, sini, pmid), r).unwrap() {
277                        pin = pmid;
278                    } else {
279                        pout = pmid;
280                    }
281                }
282                *egress = (pin+pout)/2.0;
283                *egress = *egress - egress.floor();
284                if *egress < *ingress {
285                    *egress += 1.0;
286                }
287                return Ok(true);
288            } else {
289                return Ok(false);
290            }
291        } else {
292            return Ok(false);
293        }
294
295    }
296
297
298    pub fn star_eclipse(&self, r: f64, ffac: f64, iangle: f64, posn: &Vec3, delta: f64, roche: bool, star: Star, eclipses: &mut Etype) -> Result<(), RocheError> {
299        let ri = iangle.to_radians();
300        let (sini, cosi) = ri.sin_cos();
301        let cofm = match star {
302            Star::Primary => Vec3::cofm1(),
303            Star::Secondary => Vec3::cofm2(),
304        };
305        let mut lam1: f64 = 0.0;
306        let mut lam2: f64 = 0.0;
307        let mut ingress: f64 = 0.0;
308        let mut egress: f64 = 0.0;
309        // let mut eclipses = Etype::new();
310        if (roche && self.ingress_egress(ffac, iangle, delta, &posn, &mut ingress, &mut egress)?) ||
311            (!roche && sphere_eclipse(cosi, sini, &posn, &cofm, r, &mut ingress, &mut egress, &mut lam1, &mut lam2)) {
312            eclipses.push((ingress, egress));
313        }
314        Ok(())
315    }
316
317
318    pub fn pot_min(&self, cosi: f64, sini: f64, p: &Vec3, phi1: f64, phi2: f64, lam1: f64, lam2: f64, rref: f64, pref: f64, acc: f64, phi: &mut f64, lam: &mut f64) -> Result<bool, RocheError> {
319        Ok(pot_min(self.q, self.star, self.spin, cosi, sini, p, phi1, phi2, lam1, lam2, rref, pref, acc, phi, lam)?)
320    }
321
322
323    pub fn x_l1(&self) -> Result<f64, RocheError> {
324        Ok(x_l1(self.q)?)
325    }
326
327
328    pub fn x_l1_asyncronous(&self) -> Result<f64, RocheError> {
329        match self.star {
330            Star::Primary => Ok(x_l1_1(self.q, self.spin)?),
331            Star::Secondary => Ok(x_l1_2(self.q, self.spin)?),
332        }
333    }
334
335
336    pub fn x_l2(&self) -> Result<f64, RocheError> {
337        Ok(x_l2(self.q)?)
338    }
339
340
341    pub fn x_l3(&self) -> Result<f64, RocheError> {
342        Ok(x_l3(self.q)?)
343    }
344
345}
346
roche/roche_context.rs

roche/
roche_context.rs
