satkit 0.3.14

Satellite Toolkit
Documentation 
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175

//! Tsitouras Order 5(4) Runge-Kutta integrator
//!
//! See:
//! <https://dblp.org/rec/journals/cma/Tsitouras11>
//!
//! Note: in paper, sign is reversed on Bhat[7] ...
//! should be -1.0/66.0, not 1.0/66.0
//!
//! Note also: bhat values are actually error values
//! The nomenclature is confusing
//!

use super::rk_adaptive::RKAdaptive;

const A32: f64 = 0.3354806554923570;
const A42: f64 = -6.359448489975075;
const A52: f64 = -11.74888356406283;
const A43: f64 = 4.362295432869581;
const A53: f64 = 7.495539342889836;
const A54: f64 = -0.09249506636175525;
const A62: f64 = -12.92096931784711;
const A63: f64 = 8.159367898576159;
const A64: f64 = -0.07158497328140100;
const A65: f64 = -0.02826905039406838;

const BI11: f64 = -1.0530884977290216;
const BI12: f64 = -1.3299890189751412;
const BI13: f64 = -1.4364028541716351;
const BI14: f64 = 0.7139816917074209;

const BI21: f64 = 0.1017;
const BI22: f64 = -2.1966568338249754;
const BI23: f64 = 1.2949852507374631;

const BI31: f64 = 2.490627285651252793;
const BI32: f64 = -2.38535645472061657;
const BI33: f64 = 1.57803468208092486;

const BI41: f64 = -16.54810288924490272;
const BI42: f64 = -1.21712927295533244;
const BI43: f64 = -0.61620406037800089;

const BI51: f64 = 47.37952196281928122;
const BI52: f64 = -1.203071208372362603;
const BI53: f64 = -0.658047292653547382;

const BI61: f64 = -34.87065786149660974;
const BI62: f64 = -1.2;
const BI63: f64 = -2.0 / 3.0;

const BI71: f64 = 2.5;
const BI72: f64 = -1.0;
const BI73: f64 = -0.6;

pub struct RKTS54 {}

impl RKTS54 {}

impl RKAdaptive<7, 4> for RKTS54 {
    const C: [f64; 7] = [0.0, 0.161, 0.327, 0.9, 0.9800255409045097, 1.0, 1.0];

    const B: [f64; 7] = [
        0.09646076681806523,
        0.01,
        0.4798896504144996,
        1.379008574103742,
        -3.290069515436081,
        2.324710524099774,
        0.0,
    ];

    const BERR: [f64; 7] = [
        0.001780011052226,
        0.000816434459657,
        -0.007880878010262,
        0.144711007173263,
        -0.582357165452555,
        0.458082105929187,
        -1.0 / 66.0,
    ];

    const A: [[f64; 7]; 7] = [
        [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0],
        [Self::C[1], 0.0, 0.0, 0.0, 0.0, 0.0, 0.0],
        [Self::C[2] - A32, A32, 0.0, 0.0, 0.0, 0.0, 0.0],
        [Self::C[3] - A42 - A43, A42, A43, 0.0, 0.0, 0.0, 0.0],
        [Self::C[4] - A52 - A53 - A54, A52, A53, A54, 0.0, 0.0, 0.0],
        [
            Self::C[5] - A62 - A63 - A64 - A65,
            A62,
            A63,
            A64,
            A65,
            0.0,
            0.0,
        ],
        Self::B,
    ];

    const ORDER: usize = 5;

    const FSAL: bool = false;

    // From expanding expressions in Tsitorous paper...
    const BI: [[f64; 4]; 7] = [
        [
            BI11 * BI12 * BI14,
            BI11 * (BI14 + BI12 * BI13),
            BI11 * (BI13 + BI12),
            BI11,
        ],
        [0.0, BI21 * BI23, BI21 * BI22, BI21],
        [0.0, BI31 * BI33, BI31 * BI32, BI31],
        [0.0, BI41 * BI42 * BI43, BI41 * (BI42 + BI43), BI41],
        [0.0, BI51 * BI52 * BI53, BI51 * (BI52 + BI53), BI51],
        [0.0, BI61 * BI62 * BI63, BI61 * (BI62 + BI63), BI61],
        [0.0, BI71 * BI72 * BI73, BI71 * (BI72 + BI73), BI71],
    ];
}

#[cfg(test)]
mod tests {
    use super::super::types::*;
    use super::super::RKAdaptive;
    use super::super::RKAdaptiveSettings;

    use super::*;
    type State = nalgebra::Vector2<f64>;

    #[test]
    fn test_nointerp() -> ODEResult<()> {
        let y0 = State::new(1.0, 0.0);

        use std::f64::consts::PI;

        let mut settings = RKAdaptiveSettings::default();
        settings.dense_output = false;
        settings.abserror = 1e-8;
        settings.relerror = 1e-8;

        let _res = RKTS54::integrate(
            0.0,
            2.0 * PI,
            &y0,
            |_t, y| Ok([y[1], -y[0]].into()),
            &settings,
        );

        Ok(())
    }

    #[test]
    fn testit() -> ODEResult<()> {
        let y0 = State::new(1.0, 0.0);

        use std::f64::consts::PI;

        let mut settings = RKAdaptiveSettings::default();
        settings.abserror = 1e-14;
        settings.relerror = 1e-14;
        settings.dense_output = true;

        let sol = RKTS54::integrate(0.0, PI, &y0, |_t, &y| Ok([y[1], -y[0]].into()), &settings)?;

        let testcount = 100;
        (0..100).for_each(|idx| {
            let x = idx as f64 * PI / testcount as f64;
            let interp = RKTS54::interpolate(x, &sol).unwrap();
            assert!((interp[0] - x.cos()).abs() < 1e-11);
            assert!((interp[1] + x.sin()).abs() < 1e-11);
        });

        Ok(())
    }
}