1 /*
2  * Copyright 2013 Daniel Warner <contact@danrw.com>
3  *
4  * Licensed under the Apache License, Version 2.0 (the "License");
5  * you may not use this file except in compliance with the License.
6  * You may obtain a copy of the License at
7  *
8  * http://www.apache.org/licenses/LICENSE-2.0
9  *
10  * Unless required by applicable law or agreed to in writing, software
11  * distributed under the License is distributed on an "AS IS" BASIS,
12  * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
13  * See the License for the specific language governing permissions and
14  * limitations under the License.
15  */
16 
17 #include <keplerian_toolbox/third_party/libsgp4/OrbitalElements.h>
18 #include <keplerian_toolbox/third_party/libsgp4/Tle.h>
19 
OrbitalElements(const Tle & tle)20 OrbitalElements::OrbitalElements(const Tle &tle)
21 {
22     /*
23      * extract and format tle data
24      */
25     mean_anomoly_ = tle.MeanAnomaly(false);
26     ascending_node_ = tle.RightAscendingNode(false);
27     argument_perigee_ = tle.ArgumentPerigee(false);
28     eccentricity_ = tle.Eccentricity();
29     inclination_ = tle.Inclination(false);
30     mean_motion_ = tle.MeanMotion() * kTWOPI / kMINUTES_PER_DAY;
31     bstar_ = tle.BStar();
32     epoch_ = tle.Epoch();
33 
34     /*
35      * recover original mean motion (xnodp) and semimajor axis (aodp)
36      * from input elements
37      */
38     const double a1 = pow(kXKE / MeanMotion(), kTWOTHIRD);
39     const double cosio = cos(Inclination());
40     const double theta2 = cosio * cosio;
41     const double x3thm1 = 3.0 * theta2 - 1.0;
42     const double eosq = Eccentricity() * Eccentricity();
43     const double betao2 = 1.0 - eosq;
44     const double betao = sqrt(betao2);
45     const double temp = (1.5 * kCK2) * x3thm1 / (betao * betao2);
46     const double del1 = temp / (a1 * a1);
47     const double a0 = a1 * (1.0 - del1 * (1.0 / 3.0 + del1 * (1.0 + del1 * 134.0 / 81.0)));
48     const double del0 = temp / (a0 * a0);
49 
50     recovered_mean_motion_ = MeanMotion() / (1.0 + del0);
51     /*
52      * alternative way to calculate
53      * doesnt affect final results
54      * recovered_semi_major_axis_ = pow(XKE / RecoveredMeanMotion(), TWOTHIRD);
55      */
56     recovered_semi_major_axis_ = a0 / (1.0 - del0);
57 
58     /*
59      * find perigee and period
60      */
61     perigee_ = (RecoveredSemiMajorAxis() * (1.0 - Eccentricity()) - kAE) * kXKMPER;
62     period_ = kTWOPI / RecoveredMeanMotion();
63 }
64