gistools/space/propagation/
initl.rs

1use crate::{
2    space::{
3        Method, OperationMode,
4        util::{
5            constants::{J2, X2_3, XKE},
6            time::gstime,
7        },
8    },
9    util::Date,
10};
11use core::f64::consts::TAU;
12use libm::{cos, floor, pow, round, sin, sqrt};
13
14/// Options for Initl
15pub struct InitlOptions {
16    /// eccentricity of orbit
17    pub ecco: f64,
18    /// epoch of orbit
19    pub epoch: f64,
20    /// inclination of orbit
21    pub inclo: f64,
22    /// mean motion of orbit
23    pub no: f64,
24    /// satellite number
25    pub opsmode: OperationMode,
26}
27
28/// Output for Initl
29pub struct InitlOutput {
30    /// mean motion of orbit
31    pub no: f64,
32    /// method
33    pub method: Method,
34    /// orbit period
35    pub ainv: f64,
36    /// semimajor axis
37    pub ao: f64,
38    /// con41
39    pub con41: f64,
40    /// con42
41    pub con42: f64,
42    /// cosio
43    pub cosio: f64,
44    /// cosio2
45    pub cosio2: f64,
46    /// eccentricity squared
47    pub eccsq: f64,
48    /// omeosq
49    pub omeosq: f64,
50    /// posq
51    pub posq: f64,
52    /// rp
53    pub rp: f64,
54    /// rteosq
55    pub rteosq: f64,
56    /// sinio
57    pub sinio: f64,
58    /// gsto
59    pub gsto: f64,
60}
61
62/// procedure initl
63///
64/// this procedure initializes the sgp4 propagator. all the initialization is
65/// consolidated here instead of having multiple loops inside other routines.
66///
67/// author         david vallado                  719-573-2600   28 jun 2005
68///
69/// references
70/// - hoots, roehrich, norad spacetrack report #3 1980
71/// - hoots, norad spacetrack report #6 1986
72/// - hoots, schumacher and glover 2004
73/// - vallado, crawford, hujsak, kelso  2006
74///
75/// ## Parameters
76/// - `options`: initl options
77///
78/// ## Returns
79/// Initialization params for sgp4
80pub fn initl(options: InitlOptions) -> InitlOutput {
81    let InitlOptions { ecco, epoch, inclo, opsmode, mut no } = options;
82
83    // sgp4fix use old way of finding gst
84    // ----------------------- earth constants ---------------------
85    // sgp4fix identify constants and allow alternate values
86
87    // ------------- calculate auxillary epoch quantities ----------
88    let eccsq = ecco * ecco;
89    let omeosq = 1.0 - eccsq;
90    let rteosq = sqrt(omeosq);
91    let cosio = cos(inclo);
92    let cosio2 = cosio * cosio;
93
94    // ------------------ un-kozai the mean motion -----------------
95    // let ak = (XKE / no) ** X2_3;
96    let ak = pow(XKE / no, X2_3);
97    let d1 = (0.75 * J2 * (3.0 * cosio2 - 1.0)) / (rteosq * omeosq);
98    let mut del_prime = d1 / (ak * ak);
99    let adel = ak
100        * (1.0
101            - del_prime * del_prime
102            - del_prime * (1.0 / 3.0 + (134.0 * del_prime * del_prime) / 81.0));
103    del_prime = d1 / (adel * adel);
104    no /= 1.0 + del_prime;
105
106    // let ao = (XKE / no) ** X2_3;
107    let ao = pow(XKE / no, X2_3);
108    let sinio = sin(inclo);
109    let po = ao * omeosq;
110    let con42 = 1.0 - 5.0 * cosio2;
111    let con41 = -con42 - cosio2 - cosio2;
112    let ainv = 1.0 / ao;
113    let posq = po * po;
114    let rp = ao * (1.0 - ecco);
115    let method = Method::N;
116
117    //  sgp4fix modern approach to finding sidereal time
118    let mut gsto;
119    if opsmode == OperationMode::A {
120        //  sgp4fix use old way of finding gst
121        //  count integer number of days from 0 jan 1970
122        let ts70 = epoch - 7305.0;
123        let ds70 = floor(ts70 + 1.0e-8);
124        let tfrac = ts70 - ds70;
125
126        //  find greenwich location at epoch
127        let c1 = 1.720_279_169_407_036_2e-2;
128        let thgr70 = 1.7321343856509374;
129        let fk5r = 5.075_514_194_322_695e-15;
130        let c1p2p = c1 + TAU;
131        gsto = (thgr70 + c1 * ds70 + c1p2p * tfrac + ts70 * ts70 * fk5r) % TAU;
132        if gsto < 0.0 {
133            gsto += TAU;
134        }
135    } else {
136        gsto = gstime(&Date::from_time(round(epoch + 2433281.5) as i64));
137    }
138
139    InitlOutput {
140        no,
141        method,
142        ainv,
143        ao,
144        con41,
145        con42,
146        cosio,
147        cosio2,
148        eccsq,
149        omeosq,
150        posq,
151        rp,
152        rteosq,
153        sinio,
154        gsto,
155    }
156}