gistools/space/propagation/
dscom.rs

1use core::f64::consts::TAU;
2use libm::{atan2, cos, sin, sqrt};
3
4/// Options for Dscom
5#[derive(Debug, Default, Clone, Copy, PartialEq)]
6pub struct DscomOptions {
7    /// Epoch
8    pub epoch: f64,
9    /// Eccentricity
10    pub ep: f64,
11    /// Argument of perigee
12    pub argpp: f64,
13    /// True anomaly
14    pub tc: f64,
15    /// Inclination
16    pub inclp: f64,
17    /// Right ascension of ascending node
18    pub nodep: f64,
19    /// Mean motion
20    pub np: f64,
21}
22
23/// Output for Dscom
24#[derive(Debug, Default, Clone, Copy, PartialEq)]
25pub struct DscomOutput {
26    /// snodm
27    pub snodm: f64,
28    /// cnodm
29    pub cnodm: f64,
30    /// sinim
31    pub sinim: f64,
32    /// cosim
33    pub cosim: f64,
34    /// sinomm
35    pub sinomm: f64,
36
37    /// cosomm
38    pub cosomm: f64,
39    /// day
40    pub day: f64,
41    /// e3
42    pub e3: f64,
43    /// ee2
44    pub ee2: f64,
45    /// em
46    pub em: f64,
47
48    /// emsq
49    pub emsq: f64,
50    /// gam
51    pub gam: f64,
52    /// peo
53    pub peo: f64,
54    /// pgho
55    pub pgho: f64,
56    /// pho
57    pub pho: f64,
58
59    /// pinco
60    pub pinco: f64,
61    /// plo
62    pub plo: f64,
63    /// rtemsq
64    pub rtemsq: f64,
65    /// se2
66    pub se2: f64,
67    /// se3
68    pub se3: f64,
69
70    /// sgh2
71    pub sgh2: f64,
72    /// sgh3
73    pub sgh3: f64,
74    /// sgh4
75    pub sgh4: f64,
76    /// sh2
77    pub sh2: f64,
78    /// sh3
79    pub sh3: f64,
80
81    /// si2
82    pub si2: f64,
83    /// si3
84    pub si3: f64,
85    /// sl2
86    pub sl2: f64,
87    /// sl3
88    pub sl3: f64,
89    /// sl4
90    pub sl4: f64,
91
92    /// s1
93    pub s1: f64,
94    /// s2
95    pub s2: f64,
96    /// s3
97    pub s3: f64,
98    /// s4
99    pub s4: f64,
100    /// s5
101    pub s5: f64,
102
103    /// s6
104    pub s6: f64,
105    /// s7
106    pub s7: f64,
107    /// ss1
108    pub ss1: f64,
109    /// ss2
110    pub ss2: f64,
111    /// ss3
112    pub ss3: f64,
113
114    /// ss4
115    pub ss4: f64,
116    /// ss5
117    pub ss5: f64,
118    /// ss6
119    pub ss6: f64,
120    /// ss7
121    pub ss7: f64,
122    /// sz1
123    pub sz1: f64,
124
125    /// sz2
126    pub sz2: f64,
127    /// sz3
128    pub sz3: f64,
129    /// sz11
130    pub sz11: f64,
131    /// sz12
132    pub sz12: f64,
133    /// sz13
134    pub sz13: f64,
135
136    /// sz21
137    pub sz21: f64,
138    /// sz22
139    pub sz22: f64,
140    /// sz23
141    pub sz23: f64,
142    /// sz31
143    pub sz31: f64,
144    /// sz32
145    pub sz32: f64,
146
147    /// sz33
148    pub sz33: f64,
149    /// xgh2
150    pub xgh2: f64,
151    /// xgh3
152    pub xgh3: f64,
153    /// xgh4
154    pub xgh4: f64,
155    /// xh2
156    pub xh2: f64,
157
158    /// xh3
159    pub xh3: f64,
160    /// xi2
161    pub xi2: f64,
162    /// xi3
163    pub xi3: f64,
164    /// xl2
165    pub xl2: f64,
166    /// xl3
167    pub xl3: f64,
168
169    /// xl4
170    pub xl4: f64,
171    /// nm
172    pub nm: f64,
173    /// z1
174    pub z1: f64,
175    /// z2
176    pub z2: f64,
177    /// z3
178    pub z3: f64,
179
180    /// z11
181    pub z11: f64,
182    /// z12
183    pub z12: f64,
184    /// z13
185    pub z13: f64,
186    /// z21
187    pub z21: f64,
188    /// z22
189    pub z22: f64,
190
191    /// z23
192    pub z23: f64,
193    /// z31
194    pub z31: f64,
195    /// z32
196    pub z32: f64,
197    /// z33
198    pub z33: f64,
199    /// zmol
200    pub zmol: f64,
201
202    /// zmos
203    pub zmos: f64,
204}
205
206/// procedure dscom
207///
208/// this procedure provides deep space common items used by both the secular
209/// and periodics subroutines.  input is provided as shown. this routine
210/// used to be called dpper, but the functions inside weren't well organized.
211///
212/// author  david vallado                  719-573-2600   28 jun 2005
213///
214/// references
215/// - hoots, roehrich, norad spacetrack report #3 1980
216/// - hoots, norad spacetrack report #6 1986
217/// - hoots, schumacher and glover 2004
218/// - vallado, crawford, hujsak, kelso  2006
219///
220/// ## Parameters
221/// - `options`: the options
222///
223/// ## Returns
224/// Computed deep space common items used by both the secular and periodics
225pub fn dscom(options: DscomOptions) -> DscomOutput {
226    let DscomOptions { epoch, ep, argpp, tc, inclp, nodep, np } = options;
227
228    let mut a1: f64;
229    let mut a2: f64;
230    let mut a3: f64;
231    let mut a4: f64;
232    let mut a5: f64;
233    let mut a6: f64;
234    let mut a7: f64;
235    let mut a8: f64;
236    let mut a9: f64;
237    let mut a10: f64;
238    let mut cc: f64;
239    let mut x1: f64;
240    let mut x2: f64;
241    let mut x3: f64;
242    let mut x4: f64;
243    let mut x5: f64;
244    let mut x6: f64;
245    let mut x7: f64;
246    let mut x8: f64;
247    let mut zcosg: f64;
248    let mut zsing: f64;
249    let mut zcosh: f64;
250    let mut zsinh: f64;
251    let mut zcosi: f64;
252    let mut zsini: f64;
253
254    let mut ss1 = 0.;
255    let mut ss2 = 0.;
256    let mut ss3 = 0.;
257    let mut ss4 = 0.;
258    let mut ss5 = 0.;
259    let mut ss6 = 0.;
260    let mut ss7 = 0.;
261    let mut sz1 = 0.;
262    let mut sz2 = 0.;
263    let mut sz3 = 0.;
264    let mut sz11 = 0.;
265    let mut sz12 = 0.;
266    let mut sz13 = 0.;
267    let mut sz21 = 0.;
268    let mut sz22 = 0.;
269    let mut sz23 = 0.;
270    let mut sz31 = 0.;
271    let mut sz32 = 0.;
272    let mut sz33 = 0.;
273    let mut s1 = 0.;
274    let mut s2 = 0.;
275    let mut s3 = 0.;
276    let mut s4 = 0.;
277    let mut s5 = 0.;
278    let mut s6 = 0.;
279    let mut s7 = 0.;
280    let mut z1 = 0.;
281    let mut z2 = 0.;
282    let mut z3 = 0.;
283    let mut z11 = 0.;
284    let mut z12 = 0.;
285    let mut z13 = 0.;
286    let mut z21 = 0.;
287    let mut z22 = 0.;
288    let mut z23 = 0.;
289    let mut z31 = 0.;
290    let mut z32 = 0.;
291    let mut z33 = 0.;
292
293    // -------------------------- constants -------------------------
294    let zes = 0.01675;
295    let zel = 0.0549;
296    let c1ss = 2.9864797e-6;
297    let c1l = 4.7968065e-7;
298    let zsinis = 0.39785416;
299    let zcosis = 0.91744867;
300    let zcosgs = 0.1945905;
301    let zsings = -0.98088458;
302
303    //  --------------------- local variables ------------------------
304    let nm = np;
305    let em = ep;
306    let snodm = sin(nodep);
307    let cnodm = cos(nodep);
308    let sinomm = sin(argpp);
309    let cosomm = cos(argpp);
310    let sinim = sin(inclp);
311    let cosim = cos(inclp);
312    let emsq = em * em;
313    let betasq = 1.0 - emsq;
314    let rtemsq = sqrt(betasq);
315
316    //  ----------------- initialize lunar solar terms ---------------
317    let peo = 0.0;
318    let pinco = 0.0;
319    let plo = 0.0;
320    let pgho = 0.0;
321    let pho = 0.0;
322    let day = epoch + 18261.5 + tc / 1440.0;
323    let xnodce = (4.523602 - 9.2422029e-4 * day) % TAU;
324    let stem = sin(xnodce);
325    let ctem = cos(xnodce);
326    let zcosil = 0.91375164 - 0.03568096 * ctem;
327    let zsinil = sqrt(1.0 - zcosil * zcosil);
328    let zsinhl = (0.089683511 * stem) / zsinil;
329    let zcoshl = sqrt(1.0 - zsinhl * zsinhl);
330    let gam = 5.8351514 + 0.001944368 * day;
331    let mut zx = (0.39785416 * stem) / zsinil;
332    let zy = zcoshl * ctem + 0.91744867 * zsinhl * stem;
333    zx = atan2(zx, zy);
334    zx += gam - xnodce;
335    let zcosgl = cos(zx);
336    let zsingl = sin(zx);
337
338    //  ------------------------- do solar terms ---------------------
339    zcosg = zcosgs;
340    zsing = zsings;
341    zcosi = zcosis;
342    zsini = zsinis;
343    zcosh = cnodm;
344    zsinh = snodm;
345    cc = c1ss;
346    let xnoi = 1.0 / nm;
347
348    let mut lsflg = 0;
349    while lsflg < 2 {
350        lsflg += 1;
351        a1 = zcosg * zcosh + zsing * zcosi * zsinh;
352        a3 = -zsing * zcosh + zcosg * zcosi * zsinh;
353        a7 = -zcosg * zsinh + zsing * zcosi * zcosh;
354        a8 = zsing * zsini;
355        a9 = zsing * zsinh + zcosg * zcosi * zcosh;
356        a10 = zcosg * zsini;
357        a2 = cosim * a7 + sinim * a8;
358        a4 = cosim * a9 + sinim * a10;
359        a5 = -sinim * a7 + cosim * a8;
360        a6 = -sinim * a9 + cosim * a10;
361
362        x1 = a1 * cosomm + a2 * sinomm;
363        x2 = a3 * cosomm + a4 * sinomm;
364        x3 = -a1 * sinomm + a2 * cosomm;
365        x4 = -a3 * sinomm + a4 * cosomm;
366        x5 = a5 * sinomm;
367        x6 = a6 * sinomm;
368        x7 = a5 * cosomm;
369        x8 = a6 * cosomm;
370
371        z31 = 12.0 * x1 * x1 - 3.0 * x3 * x3;
372        z32 = 24.0 * x1 * x2 - 6.0 * x3 * x4;
373        z33 = 12.0 * x2 * x2 - 3.0 * x4 * x4;
374
375        z1 = 3.0 * (a1 * a1 + a2 * a2) + z31 * emsq;
376        z2 = 6.0 * (a1 * a3 + a2 * a4) + z32 * emsq;
377        z3 = 3.0 * (a3 * a3 + a4 * a4) + z33 * emsq;
378
379        z11 = -6.0 * a1 * a5 + emsq * (-24.0 * x1 * x7 - 6.0 * x3 * x5);
380        z12 = -6.0 * (a1 * a6 + a3 * a5)
381            + emsq * (-24.0 * (x2 * x7 + x1 * x8) + -6.0 * (x3 * x6 + x4 * x5));
382
383        z13 = -6.0 * a3 * a6 + emsq * (-24.0 * x2 * x8 - 6.0 * x4 * x6);
384
385        z21 = 6.0 * a2 * a5 + emsq * (24.0 * x1 * x5 - 6.0 * x3 * x7);
386        z22 = 6.0 * (a4 * a5 + a2 * a6)
387            + emsq * (24.0 * (x2 * x5 + x1 * x6) - 6.0 * (x4 * x7 + x3 * x8));
388        z23 = 6.0 * a4 * a6 + emsq * (24.0 * x2 * x6 - 6.0 * x4 * x8);
389
390        z1 = z1 + z1 + betasq * z31;
391        z2 = z2 + z2 + betasq * z32;
392        z3 = z3 + z3 + betasq * z33;
393        s3 = cc * xnoi;
394        s2 = (-0.5 * s3) / rtemsq;
395        s4 = s3 * rtemsq;
396        s1 = -15.0 * em * s4;
397        s5 = x1 * x3 + x2 * x4;
398        s6 = x2 * x3 + x1 * x4;
399        s7 = x2 * x4 - x1 * x3;
400
401        //  ----------------------- do lunar terms -------------------
402        if lsflg == 1 {
403            ss1 = s1;
404            ss2 = s2;
405            ss3 = s3;
406            ss4 = s4;
407            ss5 = s5;
408            ss6 = s6;
409            ss7 = s7;
410            sz1 = z1;
411            sz2 = z2;
412            sz3 = z3;
413            sz11 = z11;
414            sz12 = z12;
415            sz13 = z13;
416            sz21 = z21;
417            sz22 = z22;
418            sz23 = z23;
419            sz31 = z31;
420            sz32 = z32;
421            sz33 = z33;
422            zcosg = zcosgl;
423            zsing = zsingl;
424            zcosi = zcosil;
425            zsini = zsinil;
426            zcosh = zcoshl * cnodm + zsinhl * snodm;
427            zsinh = snodm * zcoshl - cnodm * zsinhl;
428            cc = c1l;
429        }
430    }
431
432    let zmol = (4.7199672 + (0.2299715 * day - gam)) % TAU;
433    let zmos = (6.2565837 + 0.017201977 * day) % TAU;
434
435    //  ------------------------ do solar terms ----------------------
436    let se2 = 2.0 * ss1 * ss6;
437    let se3 = 2.0 * ss1 * ss7;
438    let si2 = 2.0 * ss2 * sz12;
439    let si3 = 2.0 * ss2 * (sz13 - sz11);
440    let sl2 = -2.0 * ss3 * sz2;
441    let sl3 = -2.0 * ss3 * (sz3 - sz1);
442    let sl4 = -2.0 * ss3 * (-21.0 - 9.0 * emsq) * zes;
443    let sgh2 = 2.0 * ss4 * sz32;
444    let sgh3 = 2.0 * ss4 * (sz33 - sz31);
445    let sgh4 = -18.0 * ss4 * zes;
446    let sh2 = -2.0 * ss2 * sz22;
447    let sh3 = -2.0 * ss2 * (sz23 - sz21);
448
449    //  ------------------------ do lunar terms ----------------------
450    let ee2 = 2.0 * s1 * s6;
451    let e3 = 2.0 * s1 * s7;
452    let xi2 = 2.0 * s2 * z12;
453    let xi3 = 2.0 * s2 * (z13 - z11);
454    let xl2 = -2.0 * s3 * z2;
455    let xl3 = -2.0 * s3 * (z3 - z1);
456    let xl4 = -2.0 * s3 * (-21.0 - 9.0 * emsq) * zel;
457    let xgh2 = 2.0 * s4 * z32;
458    let xgh3 = 2.0 * s4 * (z33 - z31);
459    let xgh4 = -18.0 * s4 * zel;
460    let xh2 = -2.0 * s2 * z22;
461    let xh3 = -2.0 * s2 * (z23 - z21);
462
463    DscomOutput {
464        snodm,
465        cnodm,
466        sinim,
467        cosim,
468        sinomm,
469
470        cosomm,
471        day,
472        e3,
473        ee2,
474        em,
475
476        emsq,
477        gam,
478        peo,
479        pgho,
480        pho,
481
482        pinco,
483        plo,
484        rtemsq,
485        se2,
486        se3,
487
488        sgh2,
489        sgh3,
490        sgh4,
491        sh2,
492        sh3,
493
494        si2,
495        si3,
496        sl2,
497        sl3,
498        sl4,
499
500        s1,
501        s2,
502        s3,
503        s4,
504        s5,
505
506        s6,
507        s7,
508        ss1,
509        ss2,
510        ss3,
511
512        ss4,
513        ss5,
514        ss6,
515        ss7,
516        sz1,
517
518        sz2,
519        sz3,
520        sz11,
521        sz12,
522        sz13,
523
524        sz21,
525        sz22,
526        sz23,
527        sz31,
528        sz32,
529
530        sz33,
531        xgh2,
532        xgh3,
533        xgh4,
534        xh2,
535
536        xh3,
537        xi2,
538        xi3,
539        xl2,
540        xl3,
541
542        xl4,
543        nm,
544        z1,
545        z2,
546        z3,
547
548        z11,
549        z12,
550        z13,
551        z21,
552        z22,
553
554        z23,
555        z31,
556        z32,
557        z33,
558        zmol,
559
560        zmos,
561    }
562}
gistools/space/propagation/dscom.rs

gistools/space/propagation/
dscom.rs
