gistools/space/propagation/dspace.rs
1use core::f64::consts::TAU;
2use libm::{cos, fabs, sin};
3
4/// Options for Dspace
5#[derive(Debug, Default, Clone, Copy, PartialEq)]
6pub struct DspaceOptions {
7 /// irez
8 pub irez: f64,
9 /// d2201
10 pub d2201: f64,
11 /// d2211
12 pub d2211: f64,
13 /// d3210
14 pub d3210: f64,
15 /// d3222
16 pub d3222: f64,
17 /// d4410
18 pub d4410: f64,
19 /// d4422
20 pub d4422: f64,
21 /// d5220
22 pub d5220: f64,
23 /// d5232
24 pub d5232: f64,
25 /// d5421
26 pub d5421: f64,
27 /// d5433
28 pub d5433: f64,
29 /// dedt
30 pub dedt: f64,
31 /// del1
32 pub del1: f64,
33 /// del2
34 pub del2: f64,
35 /// del3
36 pub del3: f64,
37 /// didt
38 pub didt: f64,
39 /// dmdt
40 pub dmdt: f64,
41 /// dnodt
42 pub dnodt: f64,
43 /// domdt
44 pub domdt: f64,
45 /// argpo
46 pub argpo: f64,
47 /// argpdot
48 pub argpdot: f64,
49 /// tc
50 pub tc: f64,
51 /// gsto
52 pub gsto: f64,
53 /// xfact
54 pub xfact: f64,
55 /// xlamo
56 pub xlamo: f64,
57 /// no
58 pub no: f64,
59 /// a time
60 pub atime: f64,
61 /// eccentricity
62 pub em: f64,
63 /// argument of perigee
64 pub argpm: f64,
65 /// inclination
66 pub inclm: f64,
67 /// right ascension of ascending node
68 pub xli: f64,
69 /// mean anomaly
70 pub mm: f64,
71 /// mean motion
72 pub xni: f64,
73 /// right ascension of ascending node
74 pub nodem: f64,
75 /// rate of right ascension of ascending node
76 pub nm: f64,
77}
78
79/// Output from Dspace computation
80#[derive(Debug, Default, Clone, Copy, PartialEq)]
81pub struct DspaceOutput {
82 /// a time
83 pub atime: f64,
84 /// eccentricity
85 pub em: f64, // eccentricity
86 /// argument of perigee
87 pub argpm: f64, // argument of perigee
88 /// inclination
89 pub inclm: f64, // inclination
90 /// right ascension of ascending node
91 pub xli: f64,
92 /// mean anomaly
93 pub mm: f64, // mean anomaly
94 /// mean motion
95 pub xni: f64,
96 /// right ascension of ascending node
97 pub nodem: f64, // right ascension of ascending node
98 /// rate of right ascension of ascending node
99 pub dndt: f64,
100 /// mean motion
101 pub nm: f64, // mean motion
102}
103
104/// procedure dspace
105///
106/// this procedure provides deep space contributions to mean elements for
107/// perturbing third body. these effects have been averaged over one
108/// revolution of the sun and moon. for earth resonance effects, the
109/// effects have been averaged over no revolutions of the satellite.
110/// (mean motion)
111///
112/// author david vallado 719-573-2600 28 jun 2005
113///
114/// references
115/// - hoots, roehrich, norad spacetrack report #3 1980
116/// - hoots, norad spacetrack report #6 1986
117/// - hoots, schumacher and glover 2004
118/// - vallado, crawford, hujsak, kelso 2006
119///
120/// ## Parameters
121/// - `options`: options explaining how to compute
122/// - `tsince`: time since epoch
123///
124/// ## Returns
125/// Computed values
126pub fn dspace(options: DspaceOptions, tsince: f64) -> DspaceOutput {
127 let DspaceOptions {
128 irez,
129 d2201,
130 d2211,
131 d3210,
132 d3222,
133 d4410,
134 d4422,
135 d5220,
136 d5232,
137 d5421,
138 d5433,
139 dedt,
140 del1,
141 del2,
142 del3,
143 didt,
144 dmdt,
145 dnodt,
146 domdt,
147 argpo,
148 argpdot,
149 tc,
150 gsto,
151 xfact,
152 xlamo,
153 no,
154 mut atime,
155 mut em,
156 mut argpm,
157 mut inclm,
158 mut xli,
159 mut mm,
160 mut xni,
161 mut nodem,
162 mut nm,
163 } = options;
164
165 let fasx2 = 0.13130908;
166 let fasx4 = 2.8843198;
167 let fasx6 = 0.37448087;
168 let g22 = 5.7686396;
169 let g32 = 0.95240898;
170 let g44 = 1.8014998;
171 let g52 = 1.050833;
172 let g54 = 4.4108898;
173
174 // eslint-disable-next-line no-loss-of-precision
175 let rptim = 4.375_269_088_011_3e-3; // equates to 7.29211514668855e-5 rad/sec
176 let stepp = 720.0;
177 let stepn = -720.0;
178 let step2 = 259200.0;
179
180 let delt: f64;
181 let mut x2li: f64;
182 let mut x2omi: f64;
183 let xl: f64;
184 let mut xldot = 0.0;
185 let mut xnddt = 0.0;
186 let mut xndt = 0.0;
187 let mut xomi: f64;
188 let mut dndt = 0.0;
189 let mut ft = 0.0;
190
191 // ----------- calculate deep space resonance effects -----------
192 let theta = (gsto + tc * rptim) % TAU;
193 em += dedt * tsince;
194
195 inclm += didt * tsince;
196 argpm += domdt * tsince;
197 nodem += dnodt * tsince;
198 mm += dmdt * tsince;
199
200 // sgp4fix for negative inclinations
201 // the following if statement should be commented out
202 // if (inclm < 0.0)
203 // {
204 // inclm = -inclm;
205 // argpm = argpm - pi;
206 // nodem = nodem + pi;
207 // }
208
209 /* - update resonances : numerical (euler-maclaurin) integration - */
210 /* ------------------------- epoch restart ---------------------- */
211 // sgp4fix for propagator problems
212 // the following integration works for negative time steps and periods
213 // the specific changes are unknown because the original code was so convoluted
214
215 // sgp4fix take out atime = 0.0 and fix for faster operation
216
217 if irez != 0. {
218 // sgp4fix streamline check
219 if atime == 0.0 || tsince * atime <= 0.0 || fabs(tsince) < fabs(atime) {
220 atime = 0.0;
221 xni = no;
222 xli = xlamo;
223 }
224
225 // sgp4fix move check outside loop
226 if tsince > 0.0 {
227 delt = stepp;
228 } else {
229 delt = stepn;
230 }
231
232 let mut iretn = 381; // added for do loop
233 while iretn == 381 {
234 // ------------------- dot terms calculated -------------
235 // ----------- near - synchronous resonance terms -------
236 if irez != 2. {
237 xndt = del1 * sin(xli - fasx2)
238 + del2 * sin(2.0 * (xli - fasx4))
239 + del3 * sin(3.0 * (xli - fasx6));
240 xldot = xni + xfact;
241 xnddt = del1 * cos(xli - fasx2)
242 + 2.0 * del2 * cos(2.0 * (xli - fasx4))
243 + 3.0 * del3 * cos(3.0 * (xli - fasx6));
244 xnddt *= xldot;
245 } else {
246 // --------- near - half-day resonance terms --------
247 xomi = argpo + argpdot * atime;
248 x2omi = xomi + xomi;
249 x2li = xli + xli;
250 xndt = d2201 * sin(x2omi + xli - g22)
251 + d2211 * sin(xli - g22)
252 + d3210 * sin(xomi + xli - g32)
253 + d3222 * sin(-xomi + xli - g32)
254 + d4410 * sin(x2omi + x2li - g44)
255 + d4422 * sin(x2li - g44)
256 + d5220 * sin(xomi + xli - g52)
257 + d5232 * sin(-xomi + xli - g52)
258 + d5421 * sin(xomi + x2li - g54)
259 + d5433 * sin(-xomi + x2li - g54);
260 xldot = xni + xfact;
261 xnddt = d2201 * cos(x2omi + xli - g22)
262 + d2211 * cos(xli - g22)
263 + d3210 * cos(xomi + xli - g32)
264 + d3222 * cos(-xomi + xli - g32)
265 + d5220 * cos(xomi + xli - g52)
266 + d5232 * cos(-xomi + xli - g52)
267 + 2.0
268 * (d4410 * cos(x2omi + x2li - g44)
269 + d4422 * cos(x2li - g44)
270 + d5421 * cos(xomi + x2li - g54)
271 + d5433 * cos(-xomi + x2li - g54));
272 xnddt *= xldot;
273 }
274
275 // ----------------------- integrator -------------------
276 // sgp4fix move end checks to end of routine
277 if fabs(tsince - atime) >= stepp {
278 iretn = 381;
279 } else {
280 ft = tsince - atime;
281 iretn = 0;
282 }
283
284 if iretn == 381 {
285 xli += xldot * delt + xndt * step2;
286 xni += xndt * delt + xnddt * step2;
287 atime += delt;
288 }
289 }
290
291 nm = xni + xndt * ft + xnddt * ft * ft * 0.5;
292 xl = xli + xldot * ft + xndt * ft * ft * 0.5;
293 if irez != 1. {
294 mm = xl - 2.0 * nodem + 2.0 * theta;
295 dndt = nm - no;
296 } else {
297 mm = xl - nodem - argpm + theta;
298 dndt = nm - no;
299 }
300 nm = no + dndt;
301 }
302
303 DspaceOutput { atime, em, argpm, inclm, xli, mm, xni, nodem, dndt, nm }
304}