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 <keplerian_toolbox/planet/jpl_low_precision.hpp>
27 #include <keplerian_toolbox/core_functions/convert_anomalies.hpp>
28 #include <keplerian_toolbox/core_functions/par2ic.hpp>
29 #include <keplerian_toolbox/epoch.hpp>
30 #include <keplerian_toolbox/exceptions.hpp>
31
32 namespace kep_toolbox
33 {
34 namespace planet
35 {
36
37 // Data from http://ssd.jpl.nasa.gov/txt/p_elem_t1.txt
38
39 static const double mercury_el[6] = {0.38709927, 0.20563593, 7.00497902, 252.25032350, 77.45779628, 48.33076593};
40 static const double mercury_el_dot[6] = {0.00000037, 0.00001906, -0.00594749, 149472.67411175, 0.16047689, -0.12534081};
41 static const double venus_el[6] = {0.72333566, 0.00677672, 3.39467605, 181.97909950, 131.60246718, 76.67984255};
42 static const double venus_el_dot[6] = {0.00000390, -0.00004107, -0.00078890, 58517.81538729, 0.00268329, -0.27769418};
43 static const double earth_moon_el[6] = {1.00000261, 0.01671123, -0.00001531, 100.46457166, 102.93768193, 0.0};
44 static const double earth_moon_el_dot[6] = {0.00000562, -0.00004392, -0.01294668, 35999.37244981, 0.32327364, 0.0};
45 static const double mars_el[6] = {1.52371034, 0.09339410, 1.84969142, -4.55343205, -23.94362959, 49.55953891};
46 static const double mars_el_dot[6] = {0.00001847, 0.00007882, -0.00813131, 19140.30268499, 0.44441088, -0.29257343};
47 static const double jupiter_el[6] = {5.20288700, 0.04838624, 1.30439695, 34.39644051, 14.72847983, 100.47390909};
48 static const double jupiter_el_dot[6] = {-0.00011607, -0.00013253, -0.00183714, 3034.74612775, 0.21252668, 0.20469106};
49 static const double saturn_el[6] = {9.53667594, 0.05386179, 2.48599187, 49.95424423, 92.59887831, 113.66242448};
50 static const double saturn_el_dot[6] = {-0.00125060, -0.00050991, 0.00193609, 1222.49362201, -0.41897216, -0.28867794};
51 static const double uranus_el[6] = {19.18916464, 0.04725744, 0.77263783, 313.23810451, 170.95427630, 74.01692503};
52 static const double uranus_el_dot[6] = {-0.00196176, -0.00004397, -0.00242939, 428.48202785, 0.40805281, 0.04240589};
53 static const double neptune_el[6] = {30.06992276, 0.00859048, 1.77004347, -55.12002969, 44.96476227, 131.78422574};
54 static const double neptune_el_dot[6] = {0.00026291, 0.00005105, 0.00035372, 218.45945325, -0.32241464, -0.00508664};
55 static const double pluto_el[6] = {39.48211675, 0.24882730, 17.14001206, 238.92903833, 224.06891629, 110.30393684};
56 static const double pluto_el_dot[6] = {-0.00031596, 0.00005170, 0.00004818, 145.20780515, -0.04062942, -0.01183482};
57
58 /**
59 * Construct a planet from its common name (e.g. VENUS)
60 *
61 * \param[in] name a string describing a planet
62 */
jpl_lp(const std::string & name)63 jpl_lp::jpl_lp(const std::string &name) : ref_mjd2000(epoch(2451545.0, epoch::JD).mjd2000())
64 {
65 std::map<std::string, int> mapped_planets;
66 mapped_planets["mercury"] = 1;
67 mapped_planets["venus"] = 2;
68 mapped_planets["earth"] = 3;
69 mapped_planets["mars"] = 4;
70 mapped_planets["jupiter"] = 5;
71 mapped_planets["saturn"] = 6;
72 mapped_planets["uranus"] = 7;
73 mapped_planets["neptune"] = 8;
74 mapped_planets["pluto"] = 9;
75
76 double mu_central_body = 0.;
77 double mu_self = 0.;
78 double radius = 0.;
79 double safe_radius = 0.;
80 std::string lower_case_name = name;
81 boost::algorithm::to_lower(lower_case_name);
82 switch (mapped_planets[lower_case_name]) {
83 case (1): {
84 std::copy(mercury_el, mercury_el + 6, &jpl_elements[0]);
85 std::copy(mercury_el_dot, mercury_el_dot + 6, &jpl_elements_dot[0]);
86 radius = 2440000.;
87 safe_radius = 1.1;
88 mu_self = 22032e9;
89 mu_central_body = ASTRO_MU_SUN;
90 } break;
91 case (2): {
92 std::copy(venus_el, venus_el + 6, &jpl_elements[0]);
93 std::copy(venus_el_dot, venus_el_dot + 6, &jpl_elements_dot[0]);
94 radius = 6052000.;
95 safe_radius = 1.1;
96 mu_self = 324859e9;
97 mu_central_body = ASTRO_MU_SUN;
98 } break;
99 case (3): {
100 std::copy(earth_moon_el, earth_moon_el + 6, &jpl_elements[0]);
101 std::copy(earth_moon_el_dot, earth_moon_el_dot + 6, &jpl_elements_dot[0]);
102 radius = 6378000.;
103 safe_radius = 1.1;
104 mu_self = 398600.4418e9;
105 mu_central_body = ASTRO_MU_SUN;
106 } break;
107 case (4): {
108 std::copy(mars_el, mars_el + 6, &jpl_elements[0]);
109 std::copy(mars_el_dot, mars_el_dot + 6, &jpl_elements_dot[0]);
110 radius = 3397000.;
111 safe_radius = 1.1;
112 mu_self = 42828e9;
113 mu_central_body = ASTRO_MU_SUN;
114 } break;
115 case (5): {
116 std::copy(jupiter_el, jupiter_el + 6, &jpl_elements[0]);
117 std::copy(jupiter_el_dot, jupiter_el_dot + 6, &jpl_elements_dot[0]);
118 radius = 71492000.;
119 safe_radius = 9.;
120 mu_self = 126686534e9;
121 mu_central_body = ASTRO_MU_SUN;
122 } break;
123 case (6): {
124 std::copy(saturn_el, saturn_el + 6, &jpl_elements[0]);
125 std::copy(saturn_el_dot, saturn_el_dot + 6, &jpl_elements_dot[0]);
126 radius = 60330000.;
127 safe_radius = 1.1;
128 mu_self = 37931187e9;
129 mu_central_body = ASTRO_MU_SUN;
130 } break;
131 case (7): {
132 std::copy(uranus_el, uranus_el + 6, &jpl_elements[0]);
133 std::copy(uranus_el_dot, uranus_el_dot + 6, &jpl_elements_dot[0]);
134 radius = 25362000.;
135 safe_radius = 1.1;
136 mu_self = 5793939e9;
137 mu_central_body = ASTRO_MU_SUN;
138 } break;
139 case (8): {
140 std::copy(neptune_el, neptune_el + 6, &jpl_elements[0]);
141 std::copy(neptune_el_dot, neptune_el_dot + 6, &jpl_elements_dot[0]);
142 radius = 24622000.;
143 safe_radius = 1.1;
144 mu_self = 6836529e9;
145 mu_central_body = ASTRO_MU_SUN;
146 } break;
147 case (9): {
148 std::copy(pluto_el, pluto_el + 6, &jpl_elements[0]);
149 std::copy(pluto_el_dot, pluto_el_dot + 6, &jpl_elements_dot[0]);
150 radius = 1153000.;
151 safe_radius = 1.1;
152 mu_self = 871e9;
153 mu_central_body = ASTRO_MU_SUN;
154 } break;
155 default: {
156 throw_value_error(std::string("unknown planet name: ") + name);
157 }
158 }
159 set_mu_central_body(mu_central_body);
160 set_mu_self(mu_self);
161 set_radius(radius);
162 set_safe_radius(safe_radius);
163 set_name(lower_case_name);
164 }
165
166 /// Polymorphic copy constructor.
clone() const167 planet_ptr jpl_lp::clone() const
168 {
169 return planet_ptr(new jpl_lp(*this));
170 }
171
172 /// Computes the low-precision ephemerides
eph_impl(double mjd2000,array3D & r,array3D & v) const173 void jpl_lp::eph_impl(double mjd2000, array3D &r, array3D &v) const
174 {
175 if (mjd2000 <= -73048.0 || mjd2000 >= 18263.0) {
176 throw_value_error("Ephemeris are out of range [1800-2050]");
177 }
178 // algorithm from http://ssd.jpl.nasa.gov/txt/p_elem_t1.txt downloded 2013
179 array6D elements, elements2;
180 double dt = (mjd2000 - ref_mjd2000) / 36525.0; // Number of centuries passed since J2000.0
181 for (unsigned int i = 0; i < 6; ++i) {
182 elements[i] = (jpl_elements[i] + jpl_elements_dot[i] * dt);
183 }
184 elements2[0] = elements[0] * ASTRO_AU;
185 elements2[1] = elements[1];
186 elements2[2] = elements[2] * ASTRO_DEG2RAD;
187 elements2[3] = elements[5] * ASTRO_DEG2RAD;
188 elements2[4] = (elements[4] - elements[5]) * ASTRO_DEG2RAD;
189 elements2[5] = (elements[3] - elements[4]) * ASTRO_DEG2RAD;
190 elements2[5] = m2e(elements2[5], elements2[1]);
191 par2ic(elements2, get_mu_central_body(), r, v);
192 }
193
194 /// Extra informations streamed in human readable format
human_readable_extra() const195 std::string jpl_lp::human_readable_extra() const
196 {
197 std::ostringstream s;
198 s << "Ephemerides type: JPL low-precision" << std::endl;
199 return s.str();
200 }
201 }
202 } // namespace
203
204 // Serialization code
205 BOOST_CLASS_EXPORT_IMPLEMENT(kep_toolbox::planet::jpl_lp)
206 // Serialization code (END)
207