nyx_space/propagators/rk_methods/
rk.rs

1/*
2    Nyx, blazing fast astrodynamics
3    Copyright (C) 2018-onwards Christopher Rabotin <christopher.rabotin@gmail.com>
4
5    This program is free software: you can redistribute it and/or modify
6    it under the terms of the GNU Affero General Public License as published
7    by the Free Software Foundation, either version 3 of the License, or
8    (at your option) any later version.
9
10    This program is distributed in the hope that it will be useful,
11    but WITHOUT ANY WARRANTY; without even the implied warranty of
12    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
13    GNU Affero General Public License for more details.
14
15    You should have received a copy of the GNU Affero General Public License
16    along with this program.  If not, see <https://www.gnu.org/licenses/>.
17*/
18
19use super::RK;
20
21/// `CashKarp45` is a [Runge Kutta Cash Karp integrator](https://en.wikipedia.org/wiki/Cash%E2%80%93Karp_method).
22pub(crate) struct CashKarp45 {}
23
24impl RK for CashKarp45 {
25    const ORDER: u8 = 5;
26    const STAGES: usize = 6;
27    const A_COEFFS: &'static [f64] = &[
28        1.0 / 5.0,
29        3.0 / 40.0,
30        9.0 / 40.0,
31        3.0 / 10.0,
32        -9.0 / 10.0,
33        6.0 / 5.0,
34        -11.0 / 54.0,
35        5.0 / 2.0,
36        -70.0 / 27.0,
37        35.0 / 27.0,
38        1_631.0 / 55_296.0,
39        175.0 / 512.0,
40        575.0 / 13_824.0,
41        44_275.0 / 110_592.0,
42        253.0 / 4_096.0,
43    ];
44    const B_COEFFS: &'static [f64] = &[
45        37.0 / 378.0,
46        0.0,
47        250.0 / 621.0,
48        125.0 / 594.0,
49        0.0,
50        512.0 / 1_771.0,
51        2_825.0 / 27_648.0,
52        0.0,
53        18_575.0 / 48_384.0,
54        13_525.0 / 55_296.0,
55        277.0 / 14_336.0,
56        1.0 / 4.0,
57    ];
58}
59
60/// `RK4Fixed` is a fixed step RK4.
61///
62/// If initialized with an `PropOpts.with_adaptive_step`, the variable step will **not** be taken into consideration.
63#[allow(clippy::upper_case_acronyms)]
64pub(crate) struct RK4Fixed {}
65
66impl RK for RK4Fixed {
67    const ORDER: u8 = 4;
68    const STAGES: usize = 4;
69    const A_COEFFS: &'static [f64] = &[0.5, 0.0, 0.5, 0.0, 0.0, 1.0];
70    const B_COEFFS: &'static [f64] = &[
71        1.0 / 6.0,
72        1.0 / 3.0,
73        1.0 / 3.0,
74        1.0 / 6.0,
75        // NOTE: Duplicating the B coefficients for force the error to zero.
76        1.0 / 6.0,
77        1.0 / 3.0,
78        1.0 / 3.0,
79        1.0 / 6.0,
80    ];
81}
82
83const SQRT6: f64 = 2.449_489_742_783_178;
84
85/// `RK89` is a Runge Kutta 8-9 integrator.
86///
87/// Coefficients taken from GMAT `src/base/propagator/RungeKutta89.cpp`.
88#[allow(clippy::upper_case_acronyms)]
89pub(crate) struct RK89 {}
90
91impl RK for RK89 {
92    const ORDER: u8 = 9;
93    const STAGES: usize = 16;
94    const A_COEFFS: &'static [f64] = &[
95        1.0 / 12.0,
96        1.0 / 27.0,
97        2.0 / 27.0,
98        1.0 / 24.0,
99        0.0,
100        1.0 / 8.0,
101        (4.0 + 94.0 * SQRT6) / 375.0,
102        0.0,
103        (-94.0 - 84.0 * SQRT6) / 125.0,
104        (328.0 + 208.0 * SQRT6) / 375.0,
105        (9.0 - SQRT6) / 150.0,
106        0.0,
107        0.0,
108        (312.0 + 32.0 * SQRT6) / 1425.0,
109        (69.0 + 29.0 * SQRT6) / 570.0,
110        (927.0 - 347.0 * SQRT6) / 1250.0,
111        0.0,
112        0.0,
113        (-16248.0 + 7328.0 * SQRT6) / 9375.0,
114        (-489.0 + 179.0 * SQRT6) / 3750.0,
115        (14268.0 - 5798.0 * SQRT6) / 9375.0,
116        2.0 / 27.0,
117        0.0,
118        0.0,
119        0.0,
120        0.0,
121        (16.0 - SQRT6) / 54.0,
122        (16.0 + SQRT6) / 54.0,
123        19.0 / 256.0,
124        0.0,
125        0.0,
126        0.0,
127        0.0,
128        (118.0 - 23.0 * SQRT6) / 512.0,
129        (118.0 + 23.0 * SQRT6) / 512.0,
130        -9.0 / 256.0,
131        11.0 / 144.0,
132        0.0,
133        0.0,
134        0.0,
135        0.0,
136        (266.0 - SQRT6) / 864.0,
137        (266.0 + SQRT6) / 864.0,
138        -1.0 / 16.0,
139        -8.0 / 27.0,
140        (5034.0 - 271.0 * SQRT6) / 61440.0,
141        0.0,
142        0.0,
143        0.0,
144        0.0,
145        0.0,
146        (7859.0 - 1626.0 * SQRT6) / 10240.0,
147        (-2232.0 + 813.0 * SQRT6) / 20480.0,
148        (-594.0 + 271.0 * SQRT6) / 960.0,
149        (657.0 - 813.0 * SQRT6) / 5120.0,
150        (5996.0 - 3794.0 * SQRT6) / 405.0,
151        0.0,
152        0.0,
153        0.0,
154        0.0,
155        (-4342.0 - 338.0 * SQRT6) / 9.0,
156        (154_922.0 - 40458.0 * SQRT6) / 135.0,
157        (-4176.0 + 3794.0 * SQRT6) / 45.0,
158        (-340_864.0 + 242_816.0 * SQRT6) / 405.0,
159        (26304.0 - 15176.0 * SQRT6) / 45.0,
160        -26624.0 / 81.0,
161        (3793.0 + 2168.0 * SQRT6) / 103_680.0,
162        0.0,
163        0.0,
164        0.0,
165        0.0,
166        (4042.0 + 2263.0 * SQRT6) / 13824.0,
167        (-231_278.0 + 40717.0 * SQRT6) / 69120.0,
168        (7947.0 - 2168.0 * SQRT6) / 11520.0,
169        (1048.0 - 542.0 * SQRT6) / 405.0,
170        (-1383.0 + 542.0 * SQRT6) / 720.0,
171        2624.0 / 1053.0,
172        3.0 / 1664.0,
173        -137.0 / 1296.0,
174        0.0,
175        0.0,
176        0.0,
177        0.0,
178        (5642.0 - 337.0 * SQRT6) / 864.0,
179        (5642.0 + 337.0 * SQRT6) / 864.0,
180        -299.0 / 48.0,
181        184.0 / 81.0,
182        -44.0 / 9.0,
183        -5120.0 / 1053.0,
184        -11.0 / 468.0,
185        16.0 / 9.0,
186        (33617.0 - 2168.0 * SQRT6) / 518_400.0,
187        0.0,
188        0.0,
189        0.0,
190        0.0,
191        (-3846.0 + 31.0 * SQRT6) / 13824.0,
192        (155_338.0 - 52807.0 * SQRT6) / 345_600.0,
193        (-12537.0 + 2168.0 * SQRT6) / 57600.0,
194        (92.0 + 542.0 * SQRT6) / 2025.0,
195        (-1797.0 - 542.0 * SQRT6) / 3600.0,
196        320.0 / 567.0,
197        -1.0 / 1920.0,
198        4.0 / 105.0,
199        0.0,
200        (-36487.0 - 30352.0 * SQRT6) / 279_600.0,
201        0.0,
202        0.0,
203        0.0,
204        0.0,
205        (-29666.0 - 4499.0 * SQRT6) / 7456.0,
206        (2_779_182.0 - 615_973.0 * SQRT6) / 186_400.0,
207        (-94329.0 + 91056.0 * SQRT6) / 93200.0,
208        (-232_192.0 + 121_408.0 * SQRT6) / 17475.0,
209        (101_226.0 - 22764.0 * SQRT6) / 5825.0,
210        -169_984.0 / 9087.0,
211        -87.0 / 30290.0,
212        492.0 / 1165.0,
213        0.0,
214        1260.0 / 233.0,
215    ];
216    const B_COEFFS: &'static [f64] = &[
217        23.0 / 525.0,
218        0.0,
219        0.0,
220        0.0,
221        0.0,
222        0.0,
223        0.0,
224        171.0 / 1400.0,
225        86.0 / 525.0,
226        93.0 / 280.0,
227        -2048.0 / 6825.0,
228        -3.0 / 18200.0,
229        39.0 / 175.0,
230        0.0,
231        9.0 / 25.0,
232        233.0 / 4200.0,
233        // NOTE: The b_i_star are defined here as a subtraction since the coefficients come from GMAT,
234        // and GMAT actually hard codes the b_i (as `cj[]`) and the errors.
235        23.0 / 525.0 + 7.0 / 400.0,
236        0.0,
237        0.0,
238        0.0,
239        0.0,
240        0.0,
241        0.0,
242        171.0 / 1400.0 - 63.0 / 200.0,
243        86.0 / 525.0 + 14.0 / 25.0,
244        93.0 / 280.0 - 21.0 / 20.0,
245        -2048.0 / 6825.0 + 1024.0 / 975.0,
246        -3.0 / 18200.0 + 21.0 / 36400.0,
247        39.0 / 175.0 + 3.0 / 25.0,
248        9.0 / 280.0,
249        0.0,
250        0.0,
251    ];
252}
nyx_space/propagators/rk_methods/rk.rs

nyx_space/propagators/rk_methods/
rk.rs
