rust_roche/
potential.rs

1use std::f64::consts::PI;
2use crate::{Vec3, Star};
3
4
5pub fn rpot_val(q: f64, star: Star, spin: f64, earth: &Vec3, p: &Vec3, lam: f64) -> f64 {
6
7    let r: Vec3 = *p + lam* *earth;
8    match star {
9        Star::Primary => rpot1(q, spin, &r),
10        Star::Secondary => rpot2(q, spin, &r),
11        }
12}
13
14
15pub fn rpot_val_grad(q: f64, star: Star, spin: f64, earth: &Vec3, p: &Vec3, lam: f64) -> (f64, f64, f64) {
16
17        let r: Vec3 = *p + lam* *earth;
18        let d: Vec3 = match star {
19            Star::Primary => drpot1(q, spin, &r),
20            Star::Secondary => drpot2(q, spin, &r),
21        };
22        let rpot: f64 = match star {
23            Star::Primary => rpot1(q, spin, &r),
24            Star::Secondary => rpot2(q, spin, &r),
25        };
26        let ed: Vec3 = Vec3::new(earth.y, -earth.x, 0.);
27        let dphi: f64 = 2.*PI*lam*d.dot(&ed);
28        let dlam: f64 = d.dot(earth);
29
30        (rpot, dphi, dlam)
31    }
32
33
34pub fn rpot_grad(
35    q: f64,
36    star: Star,
37    spin: f64,
38    earth: &Vec3,
39    p: &Vec3,
40    lam: f64
41) -> (f64, f64) {
42
43        let r: Vec3 = *p + lam * *earth;
44
45        let d: Vec3 = match star {
46            Star::Primary => drpot1(q, spin, &r),
47            Star::Secondary => drpot2(q, spin, &r),
48        };
49
50        let ed: Vec3 = Vec3::new(earth.y, -earth.x, 0.);
51
52        let dphi: f64 = 2.*PI * lam*d.dot(&ed);
53        let dlam: f64 = d.dot(earth);
54
55        (dphi, dlam)
56}
57
58
59///
60/// rpot computes the Roche potential at a given point. This is for the standard synchronised Roche geometry
61/// \param q mass ratio = M2/M1
62/// \param p the point in question (units scaled by separation)
63/// \return the Roche potential.
64///
65pub fn rpot(q: f64, p: &Vec3) -> f64 {
66    if q <= 0. {
67        panic!("q = {} <= 0", q);
68    }
69    let mu: f64 = q/(1. + q);
70    let comp: f64 = 1. - mu;
71    let x2y2: f64 = p.x*p.x + p.y*p.y;
72    let z2: f64 = p.z*p.z;
73    let r1sq: f64 = x2y2+z2;
74    let r1: f64 = r1sq.sqrt();
75    let r2: f64 = (r1sq + 1. - 2.*p.x).sqrt();
76    -comp/r1 - mu/r2 - (x2y2 + mu*(mu - 2.*p.x))/2.
77}
78
79
80///
81/// rpot1 computes the Roche potential at a given point allowing for non-synchronous rotation of the primary.
82/// \param q mass ratio = M2/M1
83/// \paran spin ratio spin/orbital frequencies
84/// \param p the point in question (units scaled by separation)
85/// \return the Roche potential.
86///
87pub fn rpot1(q: f64, spin: f64, p: &Vec3) -> f64 {
88    if q <= 0. {
89        panic!("q = {} <= 0", q);
90    }
91    let mu: f64 = q/(1. + q);
92    let comp: f64 = 1. - mu;
93    let x2y2: f64 = p.x*p.x + p.y*p.y;
94    let z2: f64 = p.z*p.z;
95    let r1sq: f64 = x2y2+z2;
96    let r1: f64 = r1sq.sqrt();
97    let r2: f64 = (r1sq + 1. - 2.*p.x).sqrt();
98    -comp/r1 - mu/r2 - spin*spin*x2y2/2. + mu*p.x
99}
100
101
102///
103/// rpot2 computes the Roche potential at a given point allowing for non-synchronous rotation of the secondary.
104/// \param q mass ratio = M2/M1
105/// \paran spin ratio spin/orbital frequencies
106/// \param p the point in question (units scaled by separation)
107/// \return the Roche potential.
108///
109pub fn rpot2(q: f64, spin: f64, p: &Vec3) -> f64 {
110    if q <= 0. {
111        panic!("q = {} <= 0", q);
112    }
113    let mu: f64 = q/(1. + q);
114    let comp: f64 = 1. - mu;
115    let x2y2: f64 = p.x*p.x + p.y*p.y;
116    let z2: f64 = p.z*p.z;
117    let r1sq: f64 = x2y2+z2;
118    let r1: f64 = r1sq.sqrt();
119    let r2: f64 = (r1sq + 1. - 2.*p.x).sqrt();
120    -comp/r1 - mu/r2 - spin*spin*(0.5 + 0.5*x2y2 - p.x) - comp*p.x
121}
122
123
124///
125/// drpot computes partial derivatives of Roche potential with respect
126/// to position at p for mass ratio q.
127/// \param q mass ratio  = M2/M1
128/// \param p position (units of separation)
129///
130pub fn drpot(q: f64, p: &Vec3) -> Vec3 {
131    if q <= 0. {
132        panic!("q = {} <= 0", q);
133    }
134    let r1sq: f64 = p.sqr();
135    let r1: f64 = r1sq.sqrt();
136    let r2sq: f64 = r1sq + 1. - 2.*p.x;
137    let r2: f64 = r2sq.sqrt();
138    let mu: f64 = q/(1. + q);
139    let mu1: f64 = mu/r2/r2sq;
140    let comp: f64 = (1. - mu)/r1/r1sq;
141    Vec3::new(comp*p.x + mu1*(p.x - 1.) - p.x + mu,
142              comp*p.y + mu1*p.y - p.y,
143              comp*p.z + mu1*p.z)
144}
145
146
147///
148/// drpot computes partial derivatives of Roche potential with respect
149/// to position at p for mass ratio q.
150/// \param q mass ratio  = M2/M1
151/// \param p position (units of separation)
152///
153pub fn drpot1(q: f64, spin: f64, p: &Vec3) -> Vec3 {
154    if q <= 0. {
155        panic!("q = {} <= 0", q);
156    }
157    let r1sq: f64 = p.sqr();
158    let r1: f64 = r1sq.sqrt();
159    let r2sq: f64 = r1sq + 1. - 2.*p.x;
160    let r2: f64 = r2sq.sqrt();
161    let mu: f64 = q/(1. + q);
162    let mu1: f64 = mu/r2/r2sq;
163    let comp: f64 = (1. - mu)/r1/r1sq;
164    let ssq: f64 = spin*spin;
165    Vec3::new(comp*p.x + mu1*(p.x - 1.) - ssq*p.x + mu,
166              comp*p.y + mu1*p.y - ssq*p.y,
167              comp*p.z + mu1*p.z)
168}
169
170
171///
172/// drpot2 computes partial derivatives of asynchronous Roche potential with respect
173/// to position at p for mass ratio q, star 2.
174/// \param q mass ratio  = M2/M1
175/// \param spin ratio of spin to orbital frequency.
176/// \param p position (units of separation)
177///
178pub fn drpot2(q: f64, spin: f64, p: &Vec3) -> Vec3 {
179    if q <= 0. {
180        panic!("q = {} <= 0", q);
181    }
182    let r1sq: f64 = p.sqr();
183    let r1: f64 = r1sq.sqrt();
184    let r2sq: f64 = r1sq + 1. - 2.*p.x;
185    let r2: f64 = r2sq.sqrt();
186    let mu: f64 = q/(1. + q);
187    let mu1: f64 = mu/r2/r2sq;
188    let comp: f64 = (1. - mu)/r1/r1sq;
189    let ssq: f64 = spin*spin;
190    Vec3::new(comp*p.x + mu1*(p.x - 1.) - ssq*(p.x - 1.) + mu - 1.,
191              comp*p.y + mu1*p.y - ssq*p.y,
192              comp*p.z + mu1*p.z)
193}
rust_roche/potential.rs

rust_roche/
potential.rs
