1 /*****************************************************************************
2  *   Copyright (C) 2004-2018 The pagmo development team,                     *
3  *   Advanced Concepts Team (ACT), European Space Agency (ESA)               *
4  *   http://apps.sourceforge.net/mediawiki/pagmo                             *
5  *   http://apps.sourceforge.net/mediawiki/pagmo/index.php?title=Developers  *
6  *   http://apps.sourceforge.net/mediawiki/pagmo/index.php?title=Credits     *
7  *   act@esa.int                                                             *
8  *                                                                           *
9  *   This program is free software; you can redistribute it and/or modify    *
10  *   it under the terms of the GNU General Public License as published by    *
11  *   the Free Software Foundation; either version 3 of the License, or       *
12  *   (at your option) any later version.                                     *
13  *                                                                           *
14  *   This program is distributed in the hope that it will be useful,         *
15  *   but WITHOUT ANY WARRANTY; without even the implied warranty of          *
16  *   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the           *
17  *   GNU General Public License for more details.                            *
18  *                                                                           *
19  *   You should have received a copy of the GNU General Public License       *
20  *   along with this program; if not, write to the                           *
21  *   Free Software Foundation, Inc.,                                         *
22  *   59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.               *
23  *****************************************************************************/
24 
25 // NOTE: the order of inclusion in the first two items here is forced by these
26 // two issues:
27 // http://mail.python.org/pipermail/python-list/2004-March/907592.html
28 // http://mail.python.org/pipermail/new-bugs-announce/2011-March/010395.html
29 #if defined(_WIN32)
30 #include <Python.h>
31 #include <cmath>
32 #else
33 #include <Python.h>
34 #include <cmath>
35 #endif
36 
37 #if PY_MAJOR_VERSION < 3
38 #error Minimum supported Python version is 2.6.
39 #endif
40 
41 #include <boost/math/constants/constants.hpp>
42 #include <boost/python.hpp>
43 #include <boost/python/class.hpp>
44 #include <boost/python/converter/registry.hpp>
45 #include <boost/python/copy_const_reference.hpp>
46 #include <boost/python/def.hpp>
47 #include <boost/python/docstring_options.hpp>
48 #include <boost/python/enum.hpp>
49 #include <boost/python/module.hpp>
50 #include <boost/python/operators.hpp>
51 #include <boost/python/overloads.hpp>
52 #include <boost/python/register_ptr_to_python.hpp>
53 #include <boost/python/self.hpp>
54 #include <boost/utility.hpp>
55 
56 #include "../boost_python_container_conversions.h"
57 #include "../utils.h"
58 #include "core_docstrings.hpp"
59 #include <keplerian_toolbox/keplerian_toolbox.hpp>
60 
61 using namespace boost::python;
62 
par2ic_wrapper(const kep_toolbox::array6D & E,const double & mu)63 static inline tuple par2ic_wrapper(const kep_toolbox::array6D &E, const double &mu)
64 {
65     kep_toolbox::array3D r0, v0;
66     kep_toolbox::par2ic(E, mu, r0, v0);
67     return boost::python::make_tuple(r0, v0);
68 }
69 
ic2par_wrapper(const kep_toolbox::array3D & r0,const kep_toolbox::array3D & v0,const double & mu)70 static inline kep_toolbox::array6D ic2par_wrapper(const kep_toolbox::array3D &r0, const kep_toolbox::array3D &v0,
71                                                   const double &mu)
72 {
73     kep_toolbox::array6D E;
74     kep_toolbox::ic2par(r0, v0, mu, E);
75     return E;
76 }
77 
par2eq_wrapper(const kep_toolbox::array6D & E,const bool retrograde=false)78 static inline kep_toolbox::array6D par2eq_wrapper(const kep_toolbox::array6D &E, const bool retrograde = false)
79 {
80     kep_toolbox::array6D EQ;
81     kep_toolbox::par2eq(EQ, E, retrograde);
82     return EQ;
83 }
84 
eq2par_wrapper(const kep_toolbox::array6D & EQ,const bool retrograde=false)85 static inline kep_toolbox::array6D eq2par_wrapper(const kep_toolbox::array6D &EQ, const bool retrograde = false)
86 {
87     kep_toolbox::array6D E;
88     kep_toolbox::eq2par(E, EQ, retrograde);
89     return E;
90 }
91 
eq2ic_wrapper(const kep_toolbox::array6D & EQ,const double & mu,const bool retrograde=false)92 static inline tuple eq2ic_wrapper(const kep_toolbox::array6D &EQ, const double &mu, const bool retrograde = false)
93 {
94     kep_toolbox::array3D r0, v0;
95     kep_toolbox::eq2ic(EQ, mu, r0, v0, retrograde);
96     return boost::python::make_tuple(r0, v0);
97 }
98 
ic2eq_wrapper(const kep_toolbox::array3D & r0,const kep_toolbox::array3D & v0,const double & mu,const bool retrograde=false)99 static inline kep_toolbox::array6D ic2eq_wrapper(const kep_toolbox::array3D &r0, const kep_toolbox::array3D &v0,
100                                                  const double &mu, const bool retrograde = false)
101 {
102     kep_toolbox::array6D EQ;
103     kep_toolbox::ic2eq(r0, v0, mu, EQ, retrograde);
104     return EQ;
105 }
106 
closest_distance_wrapper(const kep_toolbox::array3D & r0,const kep_toolbox::array3D & v0,const kep_toolbox::array3D & r1,const kep_toolbox::array3D & v1,const double & mu)107 static inline tuple closest_distance_wrapper(const kep_toolbox::array3D &r0, const kep_toolbox::array3D &v0,
108                                              const kep_toolbox::array3D &r1, const kep_toolbox::array3D &v1,
109                                              const double &mu)
110 {
111     double min_d, ra;
112     kep_toolbox::closest_distance(min_d, ra, r0, v0, r1, v1, mu);
113     return boost::python::make_tuple(min_d, ra);
114 }
115 
propagate_lagrangian_wrapper(const kep_toolbox::array3D & r0,const kep_toolbox::array3D & v0,const double & t,const double & mu)116 static inline tuple propagate_lagrangian_wrapper(const kep_toolbox::array3D &r0, const kep_toolbox::array3D &v0,
117                                                  const double &t, const double &mu)
118 {
119     kep_toolbox::array3D r(r0), v(v0);
120     kep_toolbox::propagate_lagrangian_u(r, v, t, mu);
121     return boost::python::make_tuple(r, v);
122 }
123 
propagate_taylor_wrapper(const kep_toolbox::array3D & r0,const kep_toolbox::array3D & v0,const double & m0,const kep_toolbox::array3D & u,const double & t,const double & mu,const double & veff,const int & log10tolerance,const int & log10rtolerance)124 static inline tuple propagate_taylor_wrapper(const kep_toolbox::array3D &r0, const kep_toolbox::array3D &v0,
125                                              const double &m0, const kep_toolbox::array3D &u, const double &t,
126                                              const double &mu, const double &veff, const int &log10tolerance,
127                                              const int &log10rtolerance)
128 {
129     kep_toolbox::array3D r(r0), v(v0);
130     double m(m0);
131     kep_toolbox::propagate_taylor(r, v, m, u, t, mu, veff, log10tolerance, log10rtolerance);
132     return boost::python::make_tuple(kep_toolbox::array3D(r), kep_toolbox::array3D(v), double(m));
133 }
134 
propagate_taylor_disturbance_wrapper(const kep_toolbox::array3D & r0,const kep_toolbox::array3D & v0,const double & m0,const kep_toolbox::array3D & thrust,const kep_toolbox::array3D & disturbance,const double & t,const double & mu,const double & veff,const int & log10tolerance,const int & log10rtolerance)135 static inline tuple propagate_taylor_disturbance_wrapper(const kep_toolbox::array3D &r0, const kep_toolbox::array3D &v0,
136                                                          const double &m0, const kep_toolbox::array3D &thrust,
137                                                          const kep_toolbox::array3D &disturbance, const double &t,
138                                                          const double &mu, const double &veff,
139                                                          const int &log10tolerance, const int &log10rtolerance)
140 {
141     kep_toolbox::array3D r(r0), v(v0);
142     double m(m0);
143     kep_toolbox::propagate_taylor_disturbance(r, v, m, thrust, disturbance, t, mu, veff, log10tolerance,
144                                               log10rtolerance);
145     return boost::python::make_tuple(kep_toolbox::array3D(r), kep_toolbox::array3D(v), double(m));
146 }
147 
propagate_taylor_J2_wrapper(const kep_toolbox::array3D & r0,const kep_toolbox::array3D & v0,const double & m0,const kep_toolbox::array3D & u,const double & t,const double & mu,const double & veff,const double & J2RG2,const int & log10tolerance,const int & log10rtolerance)148 static inline tuple propagate_taylor_J2_wrapper(const kep_toolbox::array3D &r0, const kep_toolbox::array3D &v0,
149                                                 const double &m0, const kep_toolbox::array3D &u, const double &t,
150                                                 const double &mu, const double &veff, const double &J2RG2,
151                                                 const int &log10tolerance, const int &log10rtolerance)
152 {
153     kep_toolbox::array3D r(r0), v(v0);
154     double m(m0);
155     kep_toolbox::propagate_taylor_J2(r, v, m, u, t, mu, veff, J2RG2, log10tolerance, log10rtolerance);
156     return boost::python::make_tuple(kep_toolbox::array3D(r), kep_toolbox::array3D(v), double(m));
157 }
158 
propagate_taylor_s_wrapper(const kep_toolbox::array3D & r0,const kep_toolbox::array3D & v0,const double & m0,const kep_toolbox::array3D & u,const double & s,const double & mu=1.0,const double & veff=1.0,const double & c=1.0,const double & alpha=1.5,const int & log10tolerance=-10,const int & log10rtolerance=-10)159 static inline tuple propagate_taylor_s_wrapper(const kep_toolbox::array3D &r0, const kep_toolbox::array3D &v0,
160                                                const double &m0, const kep_toolbox::array3D &u, const double &s,
161                                                const double &mu = 1.0, const double &veff = 1.0, const double &c = 1.0,
162                                                const double &alpha = 1.5, const int &log10tolerance = -10,
163                                                const int &log10rtolerance = -10)
164 {
165     kep_toolbox::array3D r(r0), v(v0);
166     double m(m0);
167     double dt(0);
168     kep_toolbox::propagate_taylor_s(r, v, m, dt, u, s, mu, veff, c, alpha, log10tolerance, log10rtolerance);
169     return boost::python::make_tuple(kep_toolbox::array3D(r), kep_toolbox::array3D(v), double(m), double(dt));
170 }
171 
fb_con_wrapper(const kep_toolbox::array3D & vin_rel,const kep_toolbox::array3D & vout_rel,const kep_toolbox::planet::base & pl)172 static inline tuple fb_con_wrapper(const kep_toolbox::array3D &vin_rel, const kep_toolbox::array3D &vout_rel,
173                                    const kep_toolbox::planet::base &pl)
174 {
175     double eq(0), ineq(0);
176     kep_toolbox::fb_con(eq, ineq, vin_rel, vout_rel, pl);
177     return boost::python::make_tuple(double(eq), double(ineq));
178 }
179 
fb_prop_wrapper(const kep_toolbox::array3D & v_in,const kep_toolbox::array3D & v_pla,const double & rp,const double & beta,const double & mu)180 static inline kep_toolbox::array3D fb_prop_wrapper(const kep_toolbox::array3D &v_in, const kep_toolbox::array3D &v_pla,
181                                                    const double &rp, const double &beta, const double &mu)
182 {
183     kep_toolbox::array3D retval;
184     kep_toolbox::fb_prop(retval, v_in, v_pla, rp, beta, mu);
185     return retval;
186 }
187 
fb_vel_wrapper(const kep_toolbox::array3D & vin_rel,const kep_toolbox::array3D & vout_rel,const kep_toolbox::planet::base & pl)188 static inline double fb_vel_wrapper(const kep_toolbox::array3D &vin_rel, const kep_toolbox::array3D &vout_rel,
189                                     const kep_toolbox::planet::base &pl)
190 {
191     double dV(0);
192     kep_toolbox::fb_vel(dV, vin_rel, vout_rel, pl);
193     return dV;
194 }
195 
damon_wrapper(const kep_toolbox::array3D & v1,const kep_toolbox::array3D & v2,double tof)196 static inline tuple damon_wrapper(const kep_toolbox::array3D &v1, const kep_toolbox::array3D &v2, double tof)
197 {
198     kep_toolbox::array3D a1, a2;
199     double tau, dv;
200     kep_toolbox::damon_approx(v1, v2, tof, a1, a2, tau, dv);
201     return boost::python::make_tuple(kep_toolbox::array3D(a1), kep_toolbox::array3D(a2), double(tau), double(dv));
202 }
203 
204 #define get_constant(arg)                                                                                              \
205     static inline double get_##arg()                                                                                   \
206     {                                                                                                                  \
207         return ASTRO_##arg;                                                                                            \
208     }
209 
210 get_constant(AU);
211 get_constant(JR);
212 get_constant(MU_SUN);
213 get_constant(MU_EARTH);
214 get_constant(EARTH_VELOCITY);
215 get_constant(EARTH_J2);
216 get_constant(EARTH_RADIUS);
217 get_constant(DEG2RAD);
218 get_constant(RAD2DEG);
219 get_constant(DAY2SEC);
220 get_constant(SEC2DAY);
221 get_constant(DAY2YEAR);
222 get_constant(G0);
223 
224 #define PYKEP_REGISTER_CONVERTER(T, policy)                                                                            \
225     {                                                                                                                  \
226         boost::python::type_info info = boost::python::type_id<T>();                                                   \
227         const boost::python::converter::registration *reg = boost::python::converter::registry::query(info);           \
228         if (reg == NULL) {                                                                                             \
229             to_tuple_mapping<T>();                                                                                     \
230             from_python_sequence<T, policy>();                                                                         \
231         } else if ((*reg).m_to_python == NULL) {                                                                       \
232             to_tuple_mapping<T>();                                                                                     \
233             from_python_sequence<T, policy>();                                                                         \
234         }                                                                                                              \
235     }
236 
237 #define PYKEP_EXPOSE_CONSTANT(arg) def("_get_" #arg, &get_##arg);
238 
239 typedef std::array<double, 8> array8D;
240 typedef std::array<double, 11> array11D;
241 
BOOST_PYTHON_MODULE(core)242 BOOST_PYTHON_MODULE(core)
243 {
244     // Disable docstring c++ signature to allow sphinx autodoc to work properly
245     docstring_options doc_options;
246     doc_options.disable_signatures();
247 
248     // Register std converters to python lists if not already registered by some
249     // other module
250     PYKEP_REGISTER_CONVERTER(std::vector<double>, variable_capacity_policy)
251     PYKEP_REGISTER_CONVERTER(kep_toolbox::array3D, fixed_size_policy)
252     PYKEP_REGISTER_CONVERTER(kep_toolbox::array6D, fixed_size_policy)
253     PYKEP_REGISTER_CONVERTER(kep_toolbox::array7D, fixed_size_policy)
254     PYKEP_REGISTER_CONVERTER(array8D, fixed_size_policy)
255     PYKEP_REGISTER_CONVERTER(array11D, fixed_size_policy)
256     PYKEP_REGISTER_CONVERTER(std::vector<kep_toolbox::array3D>, variable_capacity_policy)
257     PYKEP_REGISTER_CONVERTER(std::vector<array8D>, variable_capacity_policy)
258     PYKEP_REGISTER_CONVERTER(std::vector<array11D>, variable_capacity_policy)
259 
260     // Expose the astrodynamical constants.
261     PYKEP_EXPOSE_CONSTANT(AU);
262     PYKEP_EXPOSE_CONSTANT(JR);
263     PYKEP_EXPOSE_CONSTANT(MU_SUN);
264     PYKEP_EXPOSE_CONSTANT(MU_EARTH);
265     PYKEP_EXPOSE_CONSTANT(EARTH_VELOCITY);
266     PYKEP_EXPOSE_CONSTANT(EARTH_J2);
267     PYKEP_EXPOSE_CONSTANT(EARTH_RADIUS);
268     PYKEP_EXPOSE_CONSTANT(DEG2RAD);
269     PYKEP_EXPOSE_CONSTANT(RAD2DEG);
270     PYKEP_EXPOSE_CONSTANT(DAY2SEC);
271     PYKEP_EXPOSE_CONSTANT(SEC2DAY);
272     PYKEP_EXPOSE_CONSTANT(DAY2YEAR);
273     PYKEP_EXPOSE_CONSTANT(G0);
274 
275     // Exposing enums
276     enum_<kep_toolbox::epoch::type>("_epoch_type", "Defines a julian date type exposing the corresponding c++ enum\n\n"
277                                                    "One of\n\n"
278                                                    "* pykep.epoch.epoch_type.JD\n"
279                                                    "* pykep.epoch.epoch_type.MJD\n"
280                                                    "* pykep.epoch.epoch_type.MJD2000\n")
281         .value("MJD", kep_toolbox::epoch::MJD)
282         .value("MJD2000", kep_toolbox::epoch::MJD2000)
283         .value("JD", kep_toolbox::epoch::JD);
284 
285     // Epoch class
286     class_<kep_toolbox::epoch>("epoch", pykep::epoch_doc().c_str(),
287                                init<optional<const double &, kep_toolbox::epoch::type>>())
288         .add_property("jd", &kep_toolbox::epoch::jd,
289                       "Returns the Julian Date\n\n"
290                       "Example::\n\n"
291                       "  jd = e.jd")
292         .add_property("mjd", &kep_toolbox::epoch::mjd,
293                       "Returns the Modified Julian Date\n\n"
294                       "Example::\n\n"
295                       "  jd = e.mjd")
296         .add_property("mjd2000", &kep_toolbox::epoch::mjd2000,
297                       "Returns the Modified Julian Date 2000\n\n"
298                       "Example::\n\n"
299                       "  jd = e.mjd2000")
300         .def(repr(self))
301         .def_pickle(python_class_pickle_suite<kep_toolbox::epoch>());
302 
303     // Epoch constructors helpers
304     def("epoch_from_string", &kep_toolbox::epoch_from_string, pykep::epoch_from_string_doc().c_str());
305     def("epoch_from_iso_string", &kep_toolbox::epoch_from_iso_string, pykep::epoch_from_iso_string_doc().c_str());
306 
307     // Lambert.
308     class_<kep_toolbox::lambert_problem>(
309         "lambert_problem", "Represents a multiple revolution Lambert's problem",
310         init<const kep_toolbox::array3D &, const kep_toolbox::array3D &, const double &, const double &, const int &,
311              const int &>(pykep::lambert_problem_doc().c_str(),
312                           (arg("r1") = kep_toolbox::array3D{1, 0, 0}, arg("r2") = kep_toolbox::array3D{0, 1, 0},
313                            arg("tof") = boost::math::constants::pi<double>() / 2., arg("mu") = 1., arg("cw") = false,
314                            arg("max_revs") = 0)))
315         .def("get_v1", &kep_toolbox::lambert_problem::get_v1, return_value_policy<copy_const_reference>(),
316              "Returns a sequence of vectors containing the velocities at r1 of "
317              "all computed solutions to the Lambert's Problem\n\n"
318              "Solutions are stored in order 0 rev, 1rev, 1rev, 2rev, 2rev, ...\n\n"
319              "Example (extracts v1 for the 0 revs solution)::\n\n"
320              "  v10 = l.get_v1()[0]")
321         .def("get_v2", &kep_toolbox::lambert_problem::get_v2, return_value_policy<copy_const_reference>(),
322              "Returns a sequence of vectors containing the velocities at r2 of "
323              "all computed solutions to the Lambert's Problem\n\n"
324              "Solutions are stored in order 0 rev, 1rev, 1rev, 2rev, 2rev, ...\n\n"
325              "Example (extracts v2 for the 0 revs solution)::\n\n"
326              "  v20 = l.get_v2()[0]")
327         .def("get_r1", &kep_toolbox::lambert_problem::get_r1, return_value_policy<copy_const_reference>(),
328              "Returns a vector containing the r1 defining the Lambert's "
329              "Problem\n\n"
330              "Example ::\n\n"
331              "  r1 = l.get_r1()")
332         .def("get_r2", &kep_toolbox::lambert_problem::get_r2, return_value_policy<copy_const_reference>(),
333              "Returns a vector containing the r2 defining the Lambert's "
334              "Problem\n\n"
335              "Example ::\n\n"
336              "  r2 = l.get_r2()")
337         .def("get_tof", &kep_toolbox::lambert_problem::get_tof, return_value_policy<copy_const_reference>(),
338              "Returns the time of flight defining the Lambert's Problem\n\n"
339              "Example::\n\n"
340              "  t = l.get_tof()")
341         .def("get_mu", &kep_toolbox::lambert_problem::get_mu, return_value_policy<copy_const_reference>(),
342              "Returns the gravitational parameter defining the Lambert's "
343              "Problem\n\n"
344              "Example::\n\n"
345              "  mu = l.get_mu()")
346         .def("get_x", &kep_toolbox::lambert_problem::get_x, return_value_policy<copy_const_reference>(),
347              "Returns a sequence containing the x values of all computed "
348              "solutions to the Lambert's Problem\n\n"
349              "Solutions are stored in order 0 rev, 1rev, 1rev, 2rev, 2rev, ...\n\n"
350              "Example (extracts x for the 0 revs solution)::\n\n"
351              "  x0 = l.get_x()[0]")
352         .def("get_iters", &kep_toolbox::lambert_problem::get_iters, return_value_policy<copy_const_reference>(),
353              "Returns a sequence containing the number of iterations employed to "
354              "compute each solution to the Lambert's Problem\n\n"
355              "Solutions are stored in order 0 rev, 1rev, 1rev, 2rev, 2rev, ...\n\n"
356              "Example (extracts the number of iterations employed for the 0 revs "
357              "solution)::\n\n"
358              "  p0 = l.get_iters()[0]")
359         .def("get_Nmax", &kep_toolbox::lambert_problem::get_Nmax,
360              "Returns the maximum number of revolutions allowed. The total "
361              "number of solution to the Lambert's problem will thus be n_sol = "
362              "Nmax*2 + 1\n\n"
363              "Example::\n\n"
364              "  Nmax = l.get_Nmax()\n"
365              "  n_sol = Nmax*2+1")
366         .def(repr(self))
367         .def_pickle(python_class_pickle_suite<kep_toolbox::lambert_problem>());
368 
369     // Lagrangian propagator for keplerian orbits
370     def("propagate_lagrangian", &propagate_lagrangian_wrapper, pykep::propagate_lagrangian_doc().c_str(),
371         (arg("r0") = kep_toolbox::array3D{1, 0, 0}, arg("v0") = kep_toolbox::array3D{0, 1, 0},
372          arg("tof") = boost::math::constants::pi<double>() / 2, arg("mu") = 1));
373 
374     // Taylor propagation of inertially constant thrust arcs
375     def("propagate_taylor", &propagate_taylor_wrapper, pykep::propagate_taylor_doc().c_str(),
376         (arg("r0") = kep_toolbox::array3D{1, 0, 0}, arg("v0") = arg("v0") = kep_toolbox::array3D{0, 1, 0},
377          arg("m0") = 100, arg("thrust") = kep_toolbox::array3D{0, 0, 0},
378          arg("tof") = boost::math::constants::pi<double>() / 2, arg("mu") = 1, arg("veff") = 1, arg("log10tol") = 1e-15,
379          arg("log10rtol") = 1e-15));
380 
381     // Taylor propagation of inertially constant thrust arcs with an inertially constant disturbance
382     def("propagate_taylor_disturbance", &propagate_taylor_disturbance_wrapper,
383         "pykep.propagate_taylor(r,v,m,thrust,disturbance,t,mu,veff,log10tol,log10rtol)\n\n"
384         "- r: start position, x,y,z\n"
385         "- v: start velocity, vx,vy,vz\n"
386         "- m: starting mass\n"
387         "- thrust: fixed inertial thrust, cartesian components\n"
388         "- disturbance: fixed inertial disturbance force, cartesian components\n"
389         "- tof: propagation time\n"
390         "- mu: central body gravity constant\n\n"
391         "- veff: the product (Isp g0) defining the engine efficiency \n\n"
392         "- log10tol: the logarithm of the absolute tolerance passed to taylor propagator \n\n"
393         "- log10rtol: the logarithm of the relative tolerance passed to taylor propagator \n\n"
394         "Returns a tuple containing r, v, and m the final position, velocity and mass after the propagation.\n\n"
395         "Example::\n\n"
396         "  r,v,m = propagate_taylor_disturbance([1,0,0],[0,1,0],100,[0,0,0],[1e-2,1e-3,1e-4],pi/2,1,1,-15,-15)",
397         (arg("r"), arg("v"), arg("m"), arg("thrust"), arg("disturbance"), arg("tof"), arg("mu"), arg("veff"),
398          arg("log10tol") = 1e-15, arg("log10rtol") = 1e-15));
399     // Taylor propagation of inertially constant thrust arcs under the J2
400     // influence
401     def("propagate_taylor_J2", &propagate_taylor_J2_wrapper,
402         "pykep.propagate_taylor_J2(r,v,m,thrust,t,mu,veff,J2RG2,log10tol,log10rtol)"
403         "\n\n"
404         "- r: start position, x,y,z\n"
405         "- v: start velocity, vx,vy,vz\n"
406         "- m: starting mass\n"
407         "- thrust: fixed inertial thrust, ux,uy,uz\n"
408         "- tof: propagation time\n"
409         "- mu: central body gravity constant\n\n"
410         "- veff: the product (Isp g0) defining the engine efficiency \n\n"
411         "- J2RG2: the product (J2 RE^2) defining the gravity J2 perturbance \n\n"
412         "- log10tol: the logarithm of the absolute tolerance passed to taylor "
413         "propagator \n\n"
414         "- log10rtol: the logarithm of the relative tolerance passed to taylor "
415         "propagator \n\n"
416         "Returns a tuple containing r, v, and m the final position, velocity and "
417         "mass after the propagation.\n\n"
418         "Example::\n\n"
419         "  r,v,m = "
420         "propagate_taylor([1,0,0],[0,1,0],100,[0,0,0],pi/2,1,1,1e-3,-15,-15)",
421         (arg("r"), arg("v"), arg("m"), arg("thrust"), arg("tof"), arg("mu"), arg("veff"), arg("log10tol") = 1e-15,
422          arg("log10rtol") = 1e-15));
423     // Taylor propagation of inertially constant thrust arcs (using the
424     // generalized sundmann variable)
425     def("propagate_taylor_s", &propagate_taylor_s_wrapper,
426         "pykep.propagate_taylor_s(r,v,m,u,s,mu,veff,c=1,alpha,log10tol,log10rtol]"
427         ")\n\n"
428         "- r: start position, x,y,z\n"
429         "- v: start velocity, vx,vy,vz\n"
430         "- m: starting mass\n"
431         "- u: fixed inertial thrust, ux,uy,uz\n"
432         "- s: propagation pseudo-time\n"
433         "- mu: central body gravity constant\n\n"
434         "- veff: the product (Isp g0) defining the engine efficiency \n\n"
435         "- c: constant coefficient in the generalized Sundmann transform \n\n"
436         "- alpha: r exponent in the generalized Sundmann transform \n\n"
437         "- log10tol: the logarithm of the absolute tolerance passed to taylor "
438         "propagator \n\n"
439         "- log10rtol: the logarithm of the relative tolerance passed to taylor "
440         "propagator \n\n"
441         "Returns a tuple containing r, v, m and t the final position, velocity, "
442         "mass and time after the propagation pseudo-time s\n\n"
443         "Example::\n\n"
444         "  r,v,m,t = "
445         "propagate_taylor_s([1,0,0],[0,1,0],100,[0,0,0],pi/"
446         "2,1.0,1.0,1.0,1.0,-10,-10)");
447 
448     // Fly-by helper functions
449     def("fb_con", &fb_con_wrapper,
450         "pykep.fb_con(vin,vout,pl)\n\n"
451         "- vin: cartesian coordinates of the relative hyperbolic velocity before "
452         "the fly-by\n"
453         "- vout: vout, cartesian coordinates of the relative hyperbolic velocity "
454         "after the fly-by\n"
455         "- pl: fly-by planet\n\n"
456         "Returns a tuple containing (eq, ineq). \n"
457         "  eq represents the violation of the equality constraint norm(vin)**2 = "
458         "norm(vout)**2.\n"
459         "  ineq represents the violation of the inequality constraint on the "
460         "hyperbola asymptote maximum deflection\n\n"
461         "Example::\n\n"
462         "  v2_eq, delta_ineq = fb_con(vin, vout, planet_ss('earth'))\n"
463         "  v2_eq = v2_eq / EARTH_VELOCITY**2\n");
464     def("fb_prop", &fb_prop_wrapper,
465         "pykep.fb_prop(v,v_pla,rp,beta,mu)\n\n"
466         "- v: spacecarft velocity before the encounter (cartesian, absolute)\n"
467         "- v-pla: planet inertial velocity at encounter (cartesian, absolute)\n"
468         "- rp: fly-by radius\n"
469         "- beta: fly-by plane orientation\n"
470         "- mu: planet gravitational constant\n\n"
471         "Returns the cartesian components of the spacecarft velocity after the "
472         "encounter. \n"
473         "Example::\n\n"
474         "  vout = fb_prop([1,0,0],[0,1,0],2,3.1415/2,1)\n");
475     def("fb_vel", &fb_vel_wrapper,
476         "pykep.fb_vel(vin,vout,pl)\n\n"
477         "- vin: cartesian coordinates of the relative hyperbolic velocity before "
478         "the fly-by\n"
479         "- vout: vout, cartesian coordinates of the relative hyperbolic velocity "
480         "after the fly-by\n"
481         "- pl: fly-by planet\n\n"
482         "Returns a number dV. \n"
483         "  dV represents the norm of the delta-V needed to make a given fly-by "
484         "possible. dV is necessary to fix the magnitude and the direction of "
485         "vout.\n\n"
486         "Example::\n\n"
487         "  dV = fb_vel(vin, vout, planet_ss('earth'))\n");
488 
489     // Basic Astrodynamics
490     def("closest_distance", &closest_distance_wrapper,
491         "pykep.closest_distance(r1,v1,r2,v2,mu = 1.0)\n\n"
492         "- r1: initial position (cartesian)\n"
493         "- v1: initial velocity (cartesian)\n"
494         "- r2: final position (cartesian)\n"
495         "- v2: final velocity (cartesian)\n"
496         "- mu: planet gravitational constant\n\n"
497         "Returns a tuple containing the closest distance along a keplerian orbit "
498         "defined by r1,v1,r2,v2 (it assumes the orbit actually exists) and the "
499         "apoapsis radius of the resulting orbit. \n"
500         "Example::\n\n"
501         "  d,ra = closest_distance([1,0,0],[0,1,0],[0,1,0],[-1,0,0],1.0)\n",
502         (arg("r1"), arg("v1"), arg("r2"), arg("v2"), arg("mu") = 1.0));
503 
504     def("barker", &kep_toolbox::barker,
505         "pykep.barker(r1,r2,mu = 1.0)\n\n"
506         "- r1: initial position (cartesian)\n"
507         "- r2: final position (cartesian)\n"
508         "- mu: gravity parameter\n\n"
509         "Returns the time of flight as evaluated by the Barker equation. \n"
510         "Example:: \n\n"
511         "  t = barker([1,0,0],[0,1,0],1.0)",
512         (arg("r1"), arg("r2"), arg("mu") = 1.0));
513 
514     def("ic2par", &ic2par_wrapper,
515         "pykep.ic2par(r,v,mu = 1.0)\n\n"
516         "- r: position (cartesian)\n"
517         "- v: velocity (cartesian)\n"
518         "- mu: gravity parameter\n\n"
519         "Returns the osculating keplerian elements a,e,i,W,w,E\n"
520         "E is the eccentric anomaly for e<1, the Gudermannian for e>1\n"
521         "a is the semi-major axis always a positive quantity.\n"
522         "NOTE: The routine is singular when the elements are not defined.\n"
523         "Example:: \n\n"
524         "  el = ic2par([1,0,0],[0,1,0],1.0)",
525         (arg("r"), arg("v"), arg("mu") = 1.0));
526 
527     def("par2ic", &par2ic_wrapper,
528         "pykep.par2ic(E,mu = 1.0)\n\n"
529         "- E: osculating keplerian elements a,e,i,W,w,E ( l, ND, rad, rad, rad, rad)\n"
530         "- mu: gravity parameter (l^3/s^2)\n\n"
531         "Returns cartesian elements from Keplerian elements\n"
532         "E is the eccentric anomaly for e<1, the Gudermannian for e>1\n"
533         "a is the semi-major axis always a positive quantity.\n"
534         "Example:: \n\n"
535         "  r, v = pk.par2ic([1,0.3,0.1,0.1,0.2,0.2], 1)",
536         (arg("E"), arg("mu") = 1.0));
537 
538     def("ic2eq", &ic2eq_wrapper,
539         "pykep.ic2eq(r,v,mu = 1.0, retrogade = False)\n\n"
540         "- r: position (cartesian)\n"
541         "- v: velocity (cartesian)\n"
542         "- mu: gravity parameter\n\n"
543         "- retrogade: uses the retrograde parameters. Default value is False.\n\n"
544         "Returns the modified equinoctial elements a(1-e^2),f,g,h,k,L\n"
545         "L is the true mean longitude\n"
546         "Example:: \n\n"
547         "  E = ic2eq(r = [1,0,0], v = [0,1,0], mu =1.0)",
548         (arg("r"), arg("v"), arg("mu") = 1.0, arg("retrograde") = false));
549 
550     def("eq2ic", &eq2ic_wrapper,
551         "pykep.eq2ic(EQ,mu = 1.0, retrograde = False)\n\n"
552         "- EQ: modified equinoctial elements a(1-e^2),f,g,h,k,L\n"
553         "- mu: gravity parameter (l^3/s^2)\n\n"
554         "Returns cartesian elements from Keplerian elements\n"
555         "L is the true longitude\n"
556         "Example:: \n\n"
557         "  r, v = pk.eq2ic(eq = [1,0.3,0.1,0.1,0.2,0.2], mu = 1, retrograde = False)",
558         (arg("eq"), arg("mu") = 1.0, arg("retrograde") = false));
559 
560     def("eq2par", &eq2par_wrapper,
561         "pykep.eq2par(EQ, retrograde = False)\n\n"
562         "- EQ: modified equinoctial elements a(1-e^2),f,g,h,k,L\n"
563         "- retrogade: uses the retrograde parameters. Default value is False.\n\n"
564         "Returns the osculating Keplerian elements a,e,i,W,w,E \n"
565         "E is the eccentric anomaly for e<1, the Gudermannian for e>1\n"
566         "a is the semi-major axis, always a positive quantity.\n"
567         "L is the true longitude\n"
568         "Example:: \n\n"
569         "  E = pk.eq2par(eq = [1,0.3,0.1,0.1,0.2,0.2], retrograde = False)",
570         (arg("eq"), arg("retrograde") = false));
571 
572     def("par2eq", &par2eq_wrapper,
573         "pykep.par2eq(E, retrograde = False)\n\n"
574         "- E: osculating keplerian elements a,e,i,W,w,E ( l, ND, rad, rad, rad, rad)\n"
575         "- retrogade: uses the retrograde parameters. Default value is False.\n\n"
576         "Returns the modified equinoctial elements a(1-e^2),f,g,h,k,L\n"
577         "L is the true mean longitude\n"
578         "Example:: \n\n"
579         "  E = pk.par2eq(E = [1.0,0.1,0.2,0.3,0.4,0.5], retrograde = False)",
580         (arg("E"), arg("retrograde") = false));
581 
582     def("damon", &damon_wrapper,
583         "pykep.damon(v1,v2,tof)\n\n"
584         "- v1: starting velocity relative to the departure body. This is\n"
585         "      computed from the Lambert solution linking initial and\n"
586         "      final position in tof\n"
587         "- v2: ending velocity relative to the arrival body. This is\n"
588         "      computed from the Lambert solution linking initial and\n"
589         "      final position in tof\n"
590         "- tof: time of flight\n\n"
591         "Returns:\n\n"
592         "- a1: acceleration vector in the first part of the transfer\n"
593         "- a2: acceleration vector in the second part of the transfer\n"
594         "- tau: duration of the first part of the transfer\n"
595         "- dv: total estimated DV\n\n"
596         "This function uses the model developed by Damon Landau (JPL)\n"
597         "during GTOC7 to compute an approximation to the low-thrust transfer\n\n"
598         "Example:: \n\n"
599         " a1, a2, tau, dv = pykep.damon(v1,v2,tof)",
600         (arg("v1"), arg("v2"), arg("tof")));
601 
602     def("max_start_mass", &kep_toolbox::max_start_mass,
603         "pykep.max_start_mass(a, dv, T_max, Isp)\n\n"
604         "- a: acceleration magnitude as computed from Damon's model\n"
605         "- dv: dv magnitude as computed from Damon's model\n"
606         "- T_max: maximum thrust of the spacecarft NEP propulsion system\n"
607         "- Isp: specific impulse of the spacecarft NEP propulsion system\n\n"
608         "Returns the estimated maximum mass at leg beginning\n"
609         "This function estimates the maximum\n"
610         "mass a NEP spacececraft can have in order for a given\n"
611         "arc to be convertable into a valid low-thrust transfer\n\n"
612         "Example:: \n\n"
613         " a,_,_,dv = pykep.damon(v1,v2,tof)\n"
614         " m = pykep.max_start_mass(norm(a), dv, T_max, Isp)");
615 
616     // For the three_impulses_approximation we need to distinguish the overloaded
617     // cases. We thus introduce new function pointers
618     double (*three_imp)(double, double, double, double, double, double, double, double, double)
619         = &kep_toolbox::three_impulses_approx;
620     double (*three_imp_2_pl)(const kep_toolbox::planet::base &, const kep_toolbox::planet::base &)
621         = &kep_toolbox::three_impulses_approx;
622     double (*three_imp_2_pl_2_ep)(const kep_toolbox::planet::base &, const kep_toolbox::planet::base &,
623                                   const kep_toolbox::epoch &, const kep_toolbox::epoch &)
624         = &kep_toolbox::three_impulses_approx;
625 
626     def("_three_impulses_approx", three_imp);
627     def("_three_impulses_approx", three_imp_2_pl);
628     def("_three_impulses_approx", three_imp_2_pl_2_ep);
629 }
630