1 /*****************************************************************************
2 * Copyright (C) 2004-2018 The pykep development team, *
3 * Advanced Concepts Team (ACT), European Space Agency (ESA) *
4 * *
5 * https://gitter.im/esa/pykep *
6 * https://github.com/esa/pykep *
7 * *
8 * act@esa.int *
9 * *
10 * This program is free software; you can redistribute it and/or modify *
11 * it under the terms of the GNU General Public License as published by *
12 * the Free Software Foundation; either version 2 of the License, or *
13 * (at your option) any later version. *
14 * *
15 * This program is distributed in the hope that it will be useful, *
16 * but WITHOUT ANY WARRANTY; without even the implied warranty of *
17 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
18 * GNU General Public License for more details. *
19 * *
20 * You should have received a copy of the GNU General Public License *
21 * along with this program; if not, write to the *
22 * Free Software Foundation, Inc., *
23 * 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. *
24 *****************************************************************************/
25
26 #include <cmath>
27 #include <string>
28
29 #include <keplerian_toolbox/core_functions/convert_anomalies.hpp>
30 #include <keplerian_toolbox/core_functions/par2ic.hpp>
31 #include <keplerian_toolbox/exceptions.hpp>
32 #include <keplerian_toolbox/planet/tle.hpp>
33 #include <keplerian_toolbox/third_party/libsgp4/Eci.h>
34 #include <keplerian_toolbox/third_party/libsgp4/Globals.h>
35 #include <keplerian_toolbox/third_party/libsgp4/SGP4.h>
36 #include <keplerian_toolbox/third_party/libsgp4/SatelliteException.h>
37 #include <keplerian_toolbox/third_party/libsgp4/Tle.h>
38 #include <keplerian_toolbox/third_party/libsgp4/TleException.h>
39
40 namespace kep_toolbox
41 {
42 namespace planet
43 {
44
45 /**
46 * Construct a planet_tle from two strings containing the two line elements
47 * \param[in] line1 first line
48 * \param[in] line2 second line
49 */
tle(const std::string & line1,const std::string & line2)50 tle::tle(const std::string &line1, const std::string &line2)
51 try : base(),
52 m_line1(line1),
53 m_line2(line2),
54 m_tle(Tle("TLE satellite", line1, line2)),
55 m_sgp4_propagator(SGP4(m_tle)) {
56 // We read the osculating elements of the satellite
57
58 array6D keplerian_elements;
59 double mu_central_body = kMU * 1E09; // (m^3/s^2)
60 double mean_motion = m_tle.MeanMotion() * 2 * kPI / kSECONDS_PER_DAY; // [rad/s]
61 keplerian_elements[0] = std::pow(mu_central_body / (mean_motion * mean_motion), 1. / 3); // a [m]
62 keplerian_elements[1] = m_tle.Eccentricity(); // e
63 keplerian_elements[2] = m_tle.Inclination(false); // i [rad]
64 keplerian_elements[3] = m_tle.RightAscendingNode(false); // Om [rad]
65 keplerian_elements[4] = m_tle.ArgumentPerigee(false); // om [rad]
66 keplerian_elements[5] = m_tle.MeanAnomaly(false); // M [rad]
67
68 std::string year_str = m_tle.IntDesignator().substr(0, 2);
69 std::string rest = m_tle.IntDesignator().substr(2);
70 int year_int;
71 // Object that have no ID assigned yet will have an empty year_str
72 try {
73 year_int = std::stoi(year_str);
74 } catch(...) {
75 year_int = -1;
76 year_str = "xx";
77 rest = "TO BE ASSIGNED";
78 }
79
80 int prefix = (year_int > 56) ? (19) : (20);
81
82 std::string object_name(std::to_string(prefix) + year_str + std::string("-") + rest);
83 set_mu_central_body(mu_central_body);
84 set_name(object_name);
85 m_ref_mjd2000 = epoch(m_tle.Epoch().ToJulian(), epoch::JD).mjd2000();
86
87 } catch (TleException &e) {
88 // std::cout << "TleException cought in planet_tle constructor" << std::endl;
89 throw_value_error(e.what());
90 } catch (SatelliteException &e) {
91 // std::cout << "SatelliteException cought in planet_tle constructor" << std::endl;
92 throw_value_error(e.what());
93 }
94
95 /// Polymorphic copy constructor.le::clone() const
clone() const96 planet_ptr tle::clone() const
97 {
98 return planet_ptr(new tle(*this));
99 }
100
eph_impl(double mjd2000,array3D & r,array3D & v) const101 void tle::eph_impl(double mjd2000, array3D &r, array3D &v) const
102 {
103 Vector position;
104 Vector velocity;
105 double minutes_since = (mjd2000 - m_ref_mjd2000) * 24 * 60;
106
107 try {
108 Eci eci = m_sgp4_propagator.FindPosition(minutes_since);
109 position = eci.Position();
110 velocity = eci.Velocity();
111 r[0] = position.x * 1000;
112 r[1] = position.y * 1000;
113 r[2] = position.z * 1000;
114 v[0] = velocity.x * 1000;
115 v[1] = velocity.y * 1000;
116 v[2] = velocity.z * 1000;
117 } catch (SatelliteException &e) {
118 // std::cout << "SatelliteException caught while computing ephemerides" << std::endl;
119 throw_value_error(e.what());
120 } catch (DecayedException &e) {
121 // std::cout << "DecayedException caught while computing ephemerides" << std::endl;
122 throw_value_error(e.what());
123 }
124 }
125
126 /// Getter for the reference mjd2000
get_ref_mjd2000() const127 double tle::get_ref_mjd2000() const
128 {
129 return m_ref_mjd2000;
130 }
131
132 /// Getter for TLE line1
get_line1() const133 std::string tle::get_line1() const
134 {
135 return m_line1;
136 }
137
138 /// Getter for TLE line2
get_line2() const139 std::string tle::get_line2() const
140 {
141 return m_line2;
142 }
143
144 /// Setter for the epoch of the TLE (workaround for 2056/2057 bug)
set_epoch(const unsigned int year,const double day)145 void tle::set_epoch(const unsigned int year, const double day)
146 {
147 m_tle.setEpoch(year, day);
148 m_sgp4_propagator.SetTle(m_tle);
149 m_ref_mjd2000 = epoch(m_tle.Epoch().ToJulian(), epoch::JD).mjd2000();
150 }
151
152 /// Extra information streamed in human readable format
human_readable_extra() const153 std::string tle::human_readable_extra() const
154 {
155 std::ostringstream s;
156 s << "Ephemerides type: SGP4 propagator" << std::endl;
157 s << "TLE epoch: " << epoch(m_ref_mjd2000, epoch::MJD2000) << std::endl;
158 s << "TLE 1: " << m_line1 << std::endl;
159 s << "TLE 2: " << m_line2 << std::endl;
160 return s.str();
161 }
162 } // namespace planet
163 } // namespace kep_toolbox
164
165 BOOST_CLASS_EXPORT_IMPLEMENT(kep_toolbox::planet::tle)
166