rust_roche/
stream_physics.rs

1use bulirsch::{self, Integrator};
2use crate::Vec3;
3use crate::x_l1;
4
5
6/// 
7/// strinit sets a particle just inside the L1 point with the 
8/// correct velocity as given in Lubow and Shu.
9///
10/// \param q mass ratio = M2/M1
11/// \param r start position returned
12/// \param v start velocity returned
13///
14pub fn strinit(q: f64) -> (Vec3, Vec3) {
15    
16    const SMALL: f64 = 1.0e-5;
17    let rl1: f64 = x_l1(q);
18    let mu: f64 = q/(1.0+q);
19    let a: f64 = (1.0-mu)/rl1.powi(3)+mu/(1.0-rl1).powi(3);
20    let lambda1: f64 = (((a-2.0) + (a*(9.0*a-8.0)).sqrt())/2.0).sqrt();
21    let m1: f64 = (lambda1*lambda1-2.0*a-1.0)/2.0/lambda1;
22
23    let r: Vec3 = Vec3::new(rl1-SMALL, -m1*SMALL, 0.0);
24    let v: Vec3 = Vec3::new(-lambda1*SMALL, -lambda1*m1*SMALL, 0.0);
25
26    (r, v)
27
28}
29
30
31/// 
32/// stradv advances a particle of given position and velocity until
33/// it reaches a specified radius. It then returns with updated position and
34/// velocity. It is up to the user not to request a value that cannot be reached.
35///
36/// \param q    mass ratio = M2/M1
37/// \param r    Initial and final position
38/// \param v    Initial and final velocity
39/// \param rad  Radius to aim for
40/// \param acc  Accuracy with which to place output point at rad.
41/// \param smax Largest time step allowed. It is possible that the
42/// routine could take such a large step that it misses
43/// the point when the stream is inside the requested
44/// radius. This allows one to control this. Typical
45/// value = 1.e-3.
46///
47/// Returns time step taken
48///
49pub fn stradv(q: f64, r: &mut Vec3, v: &mut Vec3, rad: f64, acc: f64, smax: f64) -> f64 {
50
51    const TMAX: f64 = 10.0;
52    let t_next: f64 = 1.0e-2;
53
54    let mut time: f64 = 0.0;
55
56    // let to: f64;
57    let mut ro = *r;
58    let mut vo = *v;
59    
60    // Store initial radius
61    let rinit: f64 = r.length();
62    let mut rnow: f64 = rinit;
63
64    // set up Bulirsch-Stoer integrator
65    let system = OrbitalSystem{ q: q };
66    let mut integrator = Integrator::default().with_abs_tol(1.0e-8).with_rel_tol(1.0e-8).into_adaptive();
67    // Initialise arrays
68    let mut y = ndarray::array![r.x, r.y, r.z, v.x, v.y, v.z];
69    let mut y_next = ndarray::Array::zeros(y.raw_dim());
70    
71    let mut yo = y.clone();
72    let mut delta_t = t_next.min(smax);
73    // Step until radius crossed
74    while (rinit > rad && rnow > rad) || (rinit < rad && rnow < rad) {
75        ro = *r;
76        vo = *v;
77        yo = y.clone();
78        integrator
79            .step(&system, delta_t, y.view(), y_next.view_mut())
80            .unwrap();
81        y.assign(&y_next);
82        r.set(y[0], y[1], y[2]);
83        v.set(y[3], y[4], y[5]);
84        rnow = r.length();
85        time += delta_t;
86        
87        if time > TMAX {
88            panic!("roche::stradv taken too long without crossing given radius.")
89        }
90    }
91
92    // Now refine by reinitialising and binary chopping until
93    // close enough to requested radius.
94
95    let mut lo: f64 = 0.0;
96    let mut hi: f64 = delta_t;
97    let mut rlo: f64 = ro.length();
98    let mut rhi: f64 = rnow;
99    let to: f64 = time;
100
101    while (rhi-rlo).abs() > acc {
102        delta_t = (lo+hi)/2.0;
103        y = yo.clone();
104        *r = ro;
105        *v = vo;
106        time = to;
107
108        integrator
109            .step(&system, delta_t, y.view(), y_next.view_mut())
110            .unwrap();
111        y.assign(&y_next);
112
113        r.set(y[0], y[1], y[2]);
114        v.set(y[3], y[4], y[5]);
115        rnow = r.length();
116
117        if (rhi > rad && rnow > rad) || (rhi < rad && rnow < rad) {
118            rhi = rnow;
119            hi = delta_t;
120        } else {
121            rlo = rnow;
122            lo = delta_t;
123        }
124    }
125
126    time
127
128}
129
130
131///
132/// rocacc calculates and returns the acceleration (in the rotating frame)
133/// in a Roche potential of a particle of given position and velocity.
134///
135/// \param q mass ratio = M2/M1
136/// \param r position, scaled in units of separation.
137/// \param v velocity, scaled in units of separation
138///
139pub fn rocacc(q: f64, r: &Vec3, v: &Vec3) -> (f64, f64, f64) {
140
141
142    let f1: f64 = 1.0 / (1.0+q);
143    let f2: f64 = f1*q;
144
145    let yzsq: f64 = r.y*r.y + r.z*r.z;
146    let r1sq: f64 = r.x*r.x + yzsq;
147    let r2sq: f64 = (r.x-1.0)*(r.x-1.0) + yzsq;
148    let fm1: f64 = f1/(r1sq*(r1sq.sqrt()));
149    let fm2: f64 = f2/(r2sq*(r2sq.sqrt()));
150    let fm3 = fm1+fm2;
151
152    let x: f64 = -fm3*r.x + fm2 + 2.0*v.y + r.x - f2;
153    let y: f64 = -fm3*r.y       - 2.0*v.x + r.y;
154    let z: f64 = -fm3*r.z;
155    (x, y, z)
156}
157
158
159struct OrbitalSystem {
160    q: f64,
161}
162
163impl bulirsch::System for OrbitalSystem {
164    type Float = f64;
165    
166    fn system(&self, y: bulirsch::ArrayView1<Self::Float>, mut dydt: bulirsch::ArrayViewMut1<Self::Float>) {
167        dydt[[0]] = y[[3]];
168        dydt[[1]] = y[[4]];
169        dydt[[2]] = y[[5]];
170        let r = Vec3::new(y[[0]], y[[1]], y[[2]]);
171        let v = Vec3::new(y[[3]], y[[4]], y[[5]]);
172        (dydt[[3]], dydt[[4]], dydt[[5]]) = rocacc(self.q, &r, &v);
173    }
174}
rust_roche/stream_physics.rs

rust_roche/
stream_physics.rs
