roche/
roche_shadow.rs

1use crate::errors::RocheError;
2use crate::{Vec3, Star, fblink, rpot, set_earth_iangle, x_l1};
3use crate::lobes::{rtsafe, LineRoche};
4use std::f64::consts::PI;
5use pyo3::prelude::*;
6use numpy::{IntoPyArray, PyArray1};
7
8
9pub fn roche_shadow(q: f64, iangle: f64, phi: f64, n: i32, dist: f64, acc: f64) -> Result<(Vec<f64>, Vec<f64>, Vec<bool>), RocheError> {
10
11    if q <= 0.0 {
12        return Err(RocheError::ParameterError("q must be positive.".to_string()));
13    }
14
15    if n < 1 {
16        return Err(RocheError::ParameterError("n must be > 1.".to_string()));
17    }
18
19    let mut x: Vec<f64> = Vec::with_capacity(n as usize);
20    let mut y: Vec<f64> = Vec::with_capacity(n as usize);
21    let mut shade: Vec<bool> = Vec::with_capacity(n as usize);
22
23    // Compute L1 point and critical potential there.
24    let rl1: f64 = x_l1(q)?;
25    let earth: crate::Vec3 = set_earth_iangle(iangle, phi);
26    let cofm: Vec3 = Vec3::cofm2();
27    let mut p: Vec3 = Vec3::new(rl1, 0.0, 0.0);
28    let cpot: f64 = rpot(q, &p)?;
29    let mut dirn = Vec3::new(0.0, 0.0, 0.0);
30
31    // Limits for Roche lobe computation.
32    let upper: f64 = 1.0 - rl1;
33    let lower: f64 = upper / 4.0;
34
35    // Now compute Roche lobe in regular steps of angle looking
36    // from centre of Roche lobe, then search fro the shadow in between
37    // the Roche lobe and the maximum distance
38
39    let mut r1: f64;
40    let mut r2: f64;
41
42    for i in 0..n {
43
44        // L1 point is a special case because derivative becomes zero there.
45        // lambda is set so that after i=0, there is a decent starting 
46        // multiplier
47    
48        let theta: f64 = 2.0 * PI * (i as f64) / (n as f64 - 1.0);
49        let (sin_theta, cos_theta) = theta.sin_cos();
50        let dx: f64 = -cos_theta;
51        let dy: f64 = sin_theta;
52        if i == 0 || i == n - 1 {
53            r1 = 1.0 - rl1 + acc;
54        } else {
55            // Locate critical surface using rtsafe.
56            // Based on assuming that rl1 is maximum distance
57            // from centre of mass and that at no point is the
58            // surface closer than 1/4 of this.
59            let line: LineRoche = LineRoche::new(q, Star::Secondary, dx, dy, cpot);
60            r1 = rtsafe(lower, upper, |lam| line.cost(lam), acc)? + acc;
61        }
62        r2 = dist;
63        
64        // First check status of end points
65        dirn.set(dx, dy, 0.0);
66        let mut p1: Vec3 = cofm + r1 * dirn;
67        let mut p2: Vec3 = cofm + r2 * dirn;
68
69        if !fblink(q, crate::Star::Secondary, 1.0, 1.0, acc, &earth, &p1)? {
70            x.push(p1.x);
71            y.push(p1.y);
72            shade.push(false)
73        } else if fblink(q, crate::Star::Secondary, 1.0, 1.0, acc, &earth, &p2)? {
74            x.push(p2.x);
75            y.push(p2.y);
76            shade.push(true)
77        } else {
78            while r2 - r1 > acc {
79                p = (p1 + p2) / 2.0;
80                if fblink(q, Star::Secondary, 1.0, 1.0, acc, &earth, &p)? {
81                    p1 = p;
82                    r1 = (r1 + r2) / 2.0;
83                } else {
84                    p2 = p;
85                    r2 = (r1 + r2) / 2.0;
86                }
87            }
88            p = (p1 + p2) / 2.0;
89            x.push(p.x);
90            y.push(p.y);
91            shade.push(true);
92        }
93    }
94
95    Ok((x, y, shade))
96
97}
98
99#[pyfunction]
100#[pyo3(name = "shadow", signature = (q, iangle, phi, n=200, dist=5.0, acc=1.0e-4))]
101///
102///    Compute roche shadow region in equatorial plane
103///    x,y,s = shadow(q, iangle, phi, n=200, dist=5., acc=1.e-4),
104///
105///    Arguments:
106///        q (float):  M2/M1.
107///        iangle (float): inclination
108///        phi (float): orbital phase
109///        n (int): number of points in output arrays
110///        dist (float): maximum distance to search for shadow measured
111///        acc (float): numerical accuracy parameter
112///    Returns:
113///        (np.ndarray[float]): x array of roche lobe shadow
114///        (np.ndarray[float]): y array of roche lobe shadow
115///        (np.ndarray[bool]): true/false if genuine shade or not.
116///                            The array goes all the way round and when not in shade it
117///                            will be glued to the red star. This array allows you to see if this is the case or not.
118/// 
119pub fn roche_shadow_py(py: Python, q: f64, iangle: f64, phi: f64, n: i32, dist: f64, acc: f64) -> PyResult<(Py<PyArray1<f64>>, Py<PyArray1<f64>>, Py<PyArray1<bool>>)> {
120    let (xarr, yarr, shade) = roche_shadow(q, iangle, phi, n, dist, acc)?;
121    Ok((xarr.into_pyarray(py).unbind(), yarr.into_pyarray(py).unbind(), shade.into_pyarray(py).unbind()))
122}
roche/roche_shadow.rs

roche/
roche_shadow.rs
