Numerical Propagation

The backbone of numerical propagation in pykep is based on Lagrangian coefficients for Kepler’s dynamics and Taylor numerical integration, as implemented in the Heyoka [BI21] python package, for all other cases. The state transition matrix is also available and provided, in the case of numerical integration, seamlessly via variational equations.

The main routines are listed here:

Keplerian

pykep.propagate_lagrangian(rv=[[1, 0, 0], [0, 1, 0]], tof=pi / 2, mu=1, stm=False)

Propagates (Keplerian) the state for an assigned time and computes the State Transition Matrix (if requested) using the Lagrangian coefficients.

Args:

rv (2D array-like): Cartesian components of the initial position vector and velocity [[x0, y0, z0], [v0, vy0, vz0]]. Defaults to [[1,0,0], [0,1,0]].

tof (float): time of flight. Defaults to \(\frac{\pi}{2}\).

mu (float): gravitational parameter. Defaults to 1.

stm (bool): requests the computations of the State Transition Matrix

Returns:

tuple (list, list): r and v, that is the final position and velocity after the propagation. (if stm is False) tuple (tuple (list, list), numpy.ndarray (6,6)): (r,v) and the STM. (if stm is True)

Examples:

>>> import pykep as pk
>>> import numpy as np
>>> r0 = [1,0,0]
>>> v0 = [0,1,0]
>>> tof = pi/2
>>> mu = 1
>>> [r1,v1], stm = pk.propagate_lagrangian(rv=[r0,v0], tof = tof, mu = mu, stm = True)
>>> [r1,v1] = pk.propagate_lagrangian(rv=[r0,v0], tof = tof, mu = mu, stm = False)

pykep.propagate_lagrangian_grid(rv=[[1, 0, 0], [0, 1, 0]], tofs=[pi / 2], mu=1, stm=False)

This function calls pykep.propagate_lagrangian() in a loop to provide the results on a time grid containing several epochs. The main advantage is in the API which is convenient and reminiscent of the heyoka.propagate_grid() interface.

Note that this function is not necessarily more efficient than calling pykep.propagate_lagrangian() in a loop, since there is no parallelization nor SIMD magic implemented atm. Nevertheless we offer this interface for cenvenience as it may allow more compact code.

Args:

rv (2D array-like): Cartesian components of the initial position vector and velocity [[x0, y0, z0], [v0, vy0, vz0]]. Defaults to [[1,0,0], [0,1,0]].

tof (1D array-like): time of flight. Defaults to [\(\frac{\pi}{2}\)].

mu (float): gravitational parameter. Defaults to 1.

stm (bool): requests the computations of the State Transition Matrix

Returns:

list [tuple ( list , list ) ]: For each time of flight: [r,v], that is the final position and velocity after the propagation and the flattened stm (if requested).

Examples:

>>> import pykep as pk
>>> import numpy as np
>>> r0 = [1,0,0]
>>> v0 = [0,1,0]
>>> tofs = [pi/2, pi, 3/4pi]
>>> mu = 1
>>> res = pk.propagate_lagrangian_grid(rv = [r0, v0], tofs = tofs, mu = mu, stm = True)
>>> rs = [it[0][0] for it in res]
>>> vs = [it[0][1] for it in res]
>>> stms = [it[1] for it in res]

---

Non-Keplerian

class pykep.zero_hold_kep_problem(mu=1., veff=1., tol=1e-16)

Class representing the zero_hold_kep problem. In pykep, using a term well established in electrodynamics, this is the initial value problem of a fixed inertial thrust mass-varying spacecraft orbiting a main body and described by the equations:

\[\begin{split}\left\{ \begin{array}{l} \dot{\mathbf r} = \mathbf v \\ \dot{\mathbf v} = -\frac{\mu}{r^3} \mathbf r + \frac{\mathbf T}{m} \\ \dot m = - \frac{|\mathbf T|}{I_{sp} g_0} \end{array}\right.\end{split}\]

Note

Similar and connected functionality is provided by the functions zero_hold_kep(), zero_hold_kep_var(): and zero_hold_kep_dyn():.

Args:

mu (float): central body gravitational parameter. Defaults to 1.

veff (float): propulsion system effective velocity (Isp g0). Defaults to 1.

tol (float): tolerance of the Taylor adaptive integration. Defaults to 1e-16.

Note

Units need to be consistent upon construction as well as later on when calling the propagate methods.

Examples:

>>> import pykep as pk
>>> import numpy as np
>>> mu = pk.MU_SUN
>>> veff = 3000. * pk.G0
>>> tol = 1e-14
>>> sp = pk.zero_hold_kep_problem(mu, veff, tol)
>>> sp.propagate(rvm_state = [1., 0., 0., 0., 1., 0., 1], thrust = [0., 0., 1e-8], tof = 7.32)
[0.5089647068650076, 0.8607873878989034, 0.0, -0.8607873878989032, 0.5089647068650074, 0.0, 1.0]

property mu

The central body gravity parameter.

propagate(rvm_state, thrust, tof)

zero_hold_kep problem numerical propagation. In pykep, using a term well established in electrodynamics, this is the initial value problem of a fixed inertial thrust mass-varying spacecraft orbiting a main body.

The propagation will be singular for vanishing masses (infinite acceleration) and raise an exception.

Args:

rvm_state (list (7,)): position, velocity and mass flattened into a 7D list.

thrust (list (3,)): thrust flattened into a 3D list.

tof (float): time of flight.

Returns:

list (7,) : position, velocity and mass after the numerical propagation, flattened into a 7D list.

Note

Units need to be consistent with the ones used upon constructing the instance.

Examples:

>>> import pykep as pk
>>> import numpy as np
>>> mu = pk.MU_SUN
>>> veff = 3000. * pk.G0
>>> tol = 1e-14
>>> sp = pk.zero_hold_kep_problem(mu, veff, tol)
>>> sp.propagate(rvm_state = [1., 0., 0., 0., 1., 0., 1], thrust = [0., 0., 1e-8], 7.32)
[0.5089647068650076, 0.8607873878989034, 0.0, -0.8607873878989032, 0.5089647068650074, 0.0, 1.0]

propagate_var(rvm_state, thrust, tof)

zero_hold_kep problem numerical propagation via variational equations.

In pykep, using a term well established in electrodynamics, this is the initial value problem of a fixed inertial thrust mass-varying spacecraft orbiting a main body.

The propagation will be singular for vanishing masses (infinite acceleration) and raise an exception.

It also computes the system State Transition Matrix:

\[\mathbf M = \frac{d\mathbf x_f}{d\mathbf x_0}\]

as well as the gradients of the final states with respect to the thrust direction.

\[\mathbf U = \frac{d\mathbf x_f}{d\mathbf u}\]

Args:

rvm_state (list (7,)): position, velocity and mass flattened into a 7D list.

thrust (list (3,)): thrust flattened into a 3D list.

tof (float): time of flight.

Returns:

tuple (list (7,), numpy.ndarray (7,7), numpy.ndarray (7,3)): position, velocity and mass (after propagation) flattened into a 7D list, the state transition matrix and the gradient with respect to thrust.

Note

Units need to be consistent upon construction and later on when calling the propagate methods.

Examples:

>>> import pykep as pk
>>> import numpy as np
>>> mu = pk.MU_SUN
>>> veff = 3000. * pk.G0
>>> tol = 1e-14
>>> sp = pk.zero_hold_kep_problem(mu, veff, tol)
>>> sp.propagate(state = [1., 0., 0., 0., 1., 0., 1], thrust = [0., 0., 1e-8], 7.32)

property tol

The Taylor integrator tolerance.

property veff

The effective velocity (Isp g0)

---

Taylor adaptive propagators

In addition to the above mentioned propagators, Taylor adaptive integrators are offered in pykep wrapping some of the functionalities of Heyoka [BI21] python package. Their variational version is also offered (at order one) as to be able to produce stms and, where needed, other useful quantities. Higher order variational equations can also be obtained directly using the available dynamics and using Heyoka [BI21] syntax.

Some ot the Taylor integrators are “zero hold” versions of a given dynamics, meaning that they consider a constant thrust vector in some frame and include the mass variation due to said thrust.

Some of the Taylor adaptive integrators are associated to OCPs (Optimal Control Problems) of relevance to interplanetary flight. In particular they are born when applying Pontryagin principle to dynamics of interests and result in two point boundary value problems (TPBVP). In these cases, a number of auxiliary functiond describing the control, the Hamiltonian, the switching function, etc., are also provided.

---

Two Body Problem (Kepler)

pykep.ta.get_kep(tol)

Returns a Taylor adaptive propagator (Heyoka) for the simple Keplerian dynamics retreiving one from a global cache and making a copy.

If the requested propagator was never created, a call to this function will trigger its creation, else it will return the one from a global cache, thus avoiding jitting.

In pykep, the simple keplerian dynamics is defined in Cartesian coordinates (thus it is not symplectic as not in a Hamiltonian form).

The specific dynamics used is that returned by kep_dyn().

Args:

tol (float): the tolerance of the Taylor adaptive propagator.

Returns:

hy::taylor_adaptive: The Taylor adaptive propagator.

Examples:

>>> import pykep as pk
>>> ta = pk.ta.get_kep(tol = 1e-16)
>>> ta.time = 0.
>>> ta.state[:] = [1.01238082345234, -0.0423523523454,  0.22634376321, -0.1232623614,    0.123462698209365, 0.123667064622]
>>> mu = 0.01215058560962404
>>> tof = 5.7856656782589234
>>> ta.pars[0] = mu
>>> ta.propagate_until(tof)

pykep.ta.get_kep_var(tol)

Returns a (order 1) variational Taylor adaptive propagator (Heyoka) for the simple Keplerian dynamics retreiving one from a global cache and making a copy.

If the requested propagator was never created, a call to this function will trigger its creation, else it will return the one from a global cache, thus avoiding jitting.

In pykep, the simple Keplerian dynamics is defined in Cartesian coordinates (thus it is not symplectic as not in a Hamiltonian form).

The specific dynamics used is that returned by kep_dyn().

Args:

tol (float): the tolerance of the Taylor adaptive propagator.

Returns:

hy::taylor_adaptive: The Taylor adaptive propagator.

Examples:

>>> import pykep as pk
>>> ta = pk.ta.get_kep_var(tol = 1e-16)
>>> ta.time = 0.
>>> ta.state[:] = [1.01238082345234, -0.0423523523454,  0.22634376321, -0.1232623614,    0.123462698209365, 0.123667064622]
>>> mu = 0.01215058560962404
>>> tof = 5.7856656782589234
>>> ta.pars[0] = mu
>>> ta.propagate_until(tof)

pykep.ta.kep_dyn()

The dynamics of the simple Keplerian problem (kep).

In pykep, the kep is defined in Cartesian coordinates (thus it is not symplectic as not in a Hamiltonian form).

The parameter \(\mu\) is the central body gravitational parameter.

\[\begin{split}\left\{ \begin{array}{l} \dot{\mathbf r} = \mathbf v \\ \dot{\mathbf v} = -\frac{mu}{r^3} \mathbf r \\ \end{array}\right.\end{split}\]

where \(\mu\) is the only parameter.

Returns:

list [ tuple (hy::expression, hy::expression )]: The dynamics in the form [(x, dx), …]

Circular Restricted Three Body Problem

pykep.ta.get_cr3bp(tol)

Returns a Taylor adaptive propagator (Heyoka) for the CR3BP problem retreiving one from a global cache and making a copy.

If the requested propagator was never created, a call to this function will trigger its creation, else it will return the one from a global cache, thus avoiding jitting.

In pykep, the CR3BP is defined in Cartesian coordinates (thus it is not symplectic as not in a Hamiltonian form).

The specific dynamics used is that returned by cr3bp_dyn().

Args:

tol (float): the tolerance of the Taylor adaptive propagator.

Returns:

hy::taylor_adaptive: The Taylor adaptive propagator.

Examples:

>>> import pykep as pk
>>> ta = pk.ta.get_cr3bp(tol = 1e-16)
>>> ta.time = 0.
>>> ta.state[:] = [1.01238082345234, -0.0423523523454,  0.22634376321, -0.1232623614,    0.123462698209365, 0.123667064622]
>>> mu = 0.01215058560962404
>>> tof = 5.7856656782589234
>>> ta.pars[0] = mu
>>> ta.propagate_until(tof)

pykep.ta.get_cr3bp_var(tol)

Returns a (order 1) variational Taylor adaptive propagator (Heyoka) for the CR3BP problem retreiving one from a global cache and making a copy.

If the requested propagator was never created, a call to this function will trigger its creation, else it will return the one from a global cache, thus avoiding jitting.

In pykep, the CR3BP is defined in Cartesian coordinates (thus it is not symplectic as not in a Hamiltonian form).

The specific dynamics used is that returned by cr3bp_dyn().

Args:

tol (float): the tolerance of the Taylor adaptive propagator.

Returns:

hy::taylor_adaptive: The Taylor adaptive propagator.

Examples:

>>> import pykep as pk
>>> ta = pk.ta.get_cr3bp_var(tol = 1e-16)
>>> ta.time = 0.
>>> ta.state[:] = [1.01238082345234, -0.0423523523454,  0.22634376321, -0.1232623614,    0.123462698209365, 0.123667064622]
>>> mu = 0.01215058560962404
>>> tof = 5.7856656782589234
>>> ta.pars[0] = mu
>>> ta.propagate_until(tof)

pykep.ta.cr3bp_dyn()

The dynamics of the Circular Restricted Three Body Problem (CR3BP).

In pykep, the CR3BP is defined in Cartesian coordinates (thus it is not symplectic as not in a Hamiltonian form).

The parameter \(\mu\) is defined as \(\frac{m_2}{m_1+m_2}\) where \(m_2\) is the mass of the secondary body (i.e. placed on the positive x axis).

The equations are non-dimensional with units \(L = r_{12}\) (distance between the primaries), \(M = m_1 + m_2\) (total system mass) and \(T = \sqrt{\frac{r_{12}^3}{m_1+m_2}}\) (so that the angular velocity of the primaries is \(\omega = 1\) and the period of rotation between the primaries is \(2\pi\)).

\[\begin{split}\left\{ \begin{array}{l} \dot{\mathbf r} = \mathbf v \\ \dot v_x = 2v_y + x - (1 - \mu) \frac{x + \mu}{r_1^3} - \mu \frac{x + \mu - 1}{r_2^3} \\ \dot v_y = -2 v_x + y - (1 - \mu) \frac{y}{r_1^3} - \mu \frac{y}{r_2^3} \\ \dot v_z = -(1 - \mu) \frac{z}{r_1^3} - \mu \frac{z}{r_2^3} \end{array}\right.\end{split}\]

where \(\mu\) is the only parameter.

Returns:

list [ tuple (hy::expression, hy::expression )]: The dynamics in the form [(x, dx), …]

pykep.ta.cr3bp_jacobi_C()

Jacobi constant of the CR3BP as a heyoka expression.

The Jacobi constant is a conserved quantity of the CR3BP. It is defined as:

\[C = x^2 + y^2 + 2(1 - \mu) \frac{1}{r_1} + 2\mu \frac{1}{r_2} - v_x^2 - v_y^2 - v_z^2 = 2 U - v^2\]

where \(\mu\) is the mass ratio of the primaries, \(r_1\) and \(r_2\) are the distances to the primaries.

Returns:

hy::expression: the Jacobi constant of the CR3BP as a function of \(x,y,z,v_x,v_y,v_z,\mu\).

pykep.ta.cr3bp_effective_potential_U()

Effective potential of the CR3BP as a heyoka expression. The effective potential is defined as:

\[U = \frac{1}{2} \left(x^2 + y^2 + 2(1 - \mu) \frac{1}{r_1} + 2\mu \frac{1}{r_2}\right)\]

where \(\mu\) is the mass ratio of the primaries, \(r_1\) and \(r_2\) are the distances to the primaries.

Returns:

hy::expression: the effective potential of the CR3BP as a function of \(x,y,z,v_x,v_y,v_z,\mu\).

Bicircular Problem

Introduced by Simo’ et al. in his ‘97 paper [SimoGomezJM95].

pykep.ta.get_bcp()

ta.get_cbp(tol)

Returns a Taylor adaptive propagator (Heyoka) for the Bicircular Problem (BCP) problem retreiving one from a global cache and making a copy. If the requested propagator was never created, a call to this function will trigger its creation, else it will return the one from a global cache, thus avoiding jitting.

In pykep, the BCP is defined in Cartesian coordinates (it is not symplectic nor in a Hamiltonian form). It is time-dependent and assumes the Sun is on the \(x\) axis at time zero.

The specific dynamics used is that returned by cbp_dyn().

Args:

tol (float): the tolerance of the Taylor adaptive propagator.

Returns:

hy::taylor_adaptive: The Taylor adaptive propagator.

Examples:

>>> import pykep as pk
>>> ta = pk.ta.get_bcp(tol = 1e-16)
>>> ta.time = 0.
>>> ta.state[:] = [1.01238082345234, -0.0423523523454,  0.22634376321, -0.1232623614,    0.123462698209365, 0.123667064622]
>>> mu = pk.CR3BP_MU
>>> mu_s = pk.BCP_MU_S
>>> rho_s = pk.BCP_RHO_S
>>> rho_p = pk.BCP_RHO_S
>>> tof = 5.7856656782589234
>>> ta.pars[:] = [mu, mu_s, rho_s, rho_p]
>>> ta.propagate_until(tof)

ta.get_cbp(tol)

Returns a Taylor adaptive propagator (Heyoka) for the Bicircular Problem (BCP) problem retreiving one from a global cache and making a copy. If the requested propagator was never created, a call to this function will trigger its creation, else it will return the one from a global cache, thus avoiding jitting.

In pykep, the BCP is defined in Cartesian coordinates (it is not symplectic nor in a Hamiltonian form). It is time-dependent and assumes the Sun is on the \(x\) axis at time zero.

The specific dynamics used is that returned by cbp_dyn().

Args:

tol (float): the tolerance of the Taylor adaptive propagator.

Returns:

hy::taylor_adaptive: The Taylor adaptive propagator.

Examples:

>>> import pykep as pk
>>> ta = pk.ta.get_bcp(tol = 1e-16)
>>> ta.time = 0.
>>> ta.state[:] = [1.01238082345234, -0.0423523523454,  0.22634376321, -0.1232623614,    0.123462698209365, 0.123667064622]
>>> mu = pk.CR3BP_MU
>>> mu_s = pk.BCP_MU_S
>>> rho_s = pk.BCP_RHO_S
>>> rho_p = pk.BCP_RHO_S
>>> tof = 5.7856656782589234
>>> ta.pars[:] = [mu, mu_s, rho_s, rho_p]
>>> ta.propagate_until(tof)

pykep.ta.get_bcp_var()

ta.get_cbp_var(tol)

Returns a (order 1) variational Taylor adaptive propagator (Heyoka) for the Bicircular Problem (BCP) problem retreiving one from a global cache and making a copy. If the requested propagator was never created, a call to this function will trigger its creation, else it will return the one from a global cache, thus avoiding jitting.

In pykep, the BCP is defined in Cartesian coordinates (it is not symplectic nor in a Hamiltonian form). It is time-dependent and assumes the Sun is on the \(x\) axis at time zero.

The specific dynamics used is that returned by cbp_dyn().

Args:

tol (float): the tolerance of the Taylor adaptive propagator.

Returns:

hy::taylor_adaptive: The Taylor adaptive propagator.

Examples:

>>> import pykep as pk
>>> ta = pk.ta.get_bcp_var(tol = 1e-16)
>>> ta.time = 0.
>>> ta.state[:6] = [1.01238082345234, -0.0423523523454,  0.22634376321, -0.1232623614,    0.123462698209365, 0.123667064622]
>>> mu = pk.CR3BP_MU
>>> mu_s = pk.BCP_MU_S
>>> rho_s = pk.BCP_RHO_S
>>> rho_p = pk.BCP_RHO_S
>>> tof = 5.7856656782589234
>>> ta.pars[:] = [mu, mu_s, rho_s, rho_p]
>>> ta.propagate_until(tof)

ta.get_cbp_var(tol)

Returns a (order 1) variational Taylor adaptive propagator (Heyoka) for the Bicircular Problem (BCP) problem retreiving one from a global cache and making a copy. If the requested propagator was never created, a call to this function will trigger its creation, else it will return the one from a global cache, thus avoiding jitting.

In pykep, the BCP is defined in Cartesian coordinates (it is not symplectic nor in a Hamiltonian form). It is time-dependent and assumes the Sun is on the \(x\) axis at time zero.

The specific dynamics used is that returned by cbp_dyn().

Args:

tol (float): the tolerance of the Taylor adaptive propagator.

Returns:

hy::taylor_adaptive: The Taylor adaptive propagator.

Examples:

>>> import pykep as pk
>>> ta = pk.ta.get_bcp_var(tol = 1e-16)
>>> ta.time = 0.
>>> ta.state[:6] = [1.01238082345234, -0.0423523523454,  0.22634376321, -0.1232623614,    0.123462698209365, 0.123667064622]
>>> mu = pk.CR3BP_MU
>>> mu_s = pk.BCP_MU_S
>>> rho_s = pk.BCP_RHO_S
>>> rho_p = pk.BCP_RHO_S
>>> tof = 5.7856656782589234
>>> ta.pars[:] = [mu, mu_s, rho_s, rho_p]
>>> ta.propagate_until(tof)

pykep.ta.bcp_dyn()

cbp_dyn()

The dynamics of the Bicircular Problem (BCP).

In pykep, the BCP is defined in Cartesian coordinates (it is not symplectic nor in a Hamiltonian form). It is time-dependent and assumes the Sun is on the \(x\) axis at time zero.

The equations are non-dimensional with units \(L = r_{12}\) (distance between the primaries), \(M = m_1 + m_2\) (total system mass) and \(T = \sqrt{\frac{r_{12}^3}{m_1+m_2}}\) (so that the angular velocity of the primaries is \(\omega = 1\) and the period of rotation between the primaries is \(2\pi\)). The parameter \(\mu\) is defined as \(\frac{m_2}{m_1+m_2}\) where \(m_2\) is the mass of the secondary body (i.e. placed on the positive x axis).

The equations of motion are:

\[\begin{split}\left\{ \begin{array}{l} \dot{\mathbf r} = \mathbf v \\ \dot v_x = 2v_y + x - (1 - \mu) \frac{x + \mu}{r_1^3} - \mu \frac{x + \mu - 1}{r_2^3} - \frac{\mu_s}{r_s^3} (x- \rho \cos(\omega_s t)) - \frac{\mu_s}{\rho_s^2} \cos(\omega_s t)\\ \dot v_y = -2 v_x + y - (1 - \mu) \frac{y}{r_1^3} - \mu \frac{y}{r_2^3} - \frac{\mu_s}{r_s^3} (y - \rho \sin(\omega_s t)) - \frac{\mu_s}{\rho_s^2} \sin(\omega_s t)\\ \dot v_z = -(1 - \mu) \frac{z}{r_1^3} - \mu \frac{z}{r_2^3} - \mu_s \frac{z}{r_s^3} \end{array}\right.\end{split}\]

where \(\mu, \mu_s, \rho_s\) and $omega_s$ are the system parameters. The spacecarft distance from the primaries is \(r_1, r_2\) and its distance from the Sun is \(r_s\). The Sun orbits at a distance \(\rho\) with angular velocity \(\omega_s\):.

Returns:

list [ tuple (hy::expression, hy::expression )]: The dynamics in the form [(x, dx), …]

cbp_dyn()

The dynamics of the Bicircular Problem (BCP).

In pykep, the BCP is defined in Cartesian coordinates (it is not symplectic nor in a Hamiltonian form). It is time-dependent and assumes the Sun is on the \(x\) axis at time zero.

The equations are non-dimensional with units \(L = r_{12}\) (distance between the primaries), \(M = m_1 + m_2\) (total system mass) and \(T = \sqrt{\frac{r_{12}^3}{m_1+m_2}}\) (so that the angular velocity of the primaries is \(\omega = 1\) and the period of rotation between the primaries is \(2\pi\)). The parameter \(\mu\) is defined as \(\frac{m_2}{m_1+m_2}\) where \(m_2\) is the mass of the secondary body (i.e. placed on the positive x axis).

The equations of motion are:

\[\begin{split}\left\{ \begin{array}{l} \dot{\mathbf r} = \mathbf v \\ \dot v_x = 2v_y + x - (1 - \mu) \frac{x + \mu}{r_1^3} - \mu \frac{x + \mu - 1}{r_2^3} - \frac{\mu_s}{r_s^3} (x- \rho \cos(\omega_s t)) - \frac{\mu_s}{\rho_s^2} \cos(\omega_s t)\\ \dot v_y = -2 v_x + y - (1 - \mu) \frac{y}{r_1^3} - \mu \frac{y}{r_2^3} - \frac{\mu_s}{r_s^3} (y - \rho \sin(\omega_s t)) - \frac{\mu_s}{\rho_s^2} \sin(\omega_s t)\\ \dot v_z = -(1 - \mu) \frac{z}{r_1^3} - \mu \frac{z}{r_2^3} - \mu_s \frac{z}{r_s^3} \end{array}\right.\end{split}\]

where \(\mu, \mu_s, \rho_s\) and $omega_s$ are the system parameters. The spacecarft distance from the primaries is \(r_1, r_2\) and its distance from the Sun is \(r_s\). The Sun orbits at a distance \(\rho\) with angular velocity \(\omega_s\):.

Returns:

list [ tuple (hy::expression, hy::expression )]: The dynamics in the form [(x, dx), …]

zero_hold_kep

pykep.ta.get_zero_hold_kep(tol)

Returns a Taylor adaptive propagator (Heyoka) for the zero_hold_kep dynamics retreiving one from a global cache and making a copy.

In pykep, using a term well established in electrodynamics, this is the initial value problem of a fixed inertial thrust mass-varying spacecraft orbiting a main body.

If the requested propagator was never created this will create it, else it will return the one from the global cache, thus avoiding jitting.

The dynamics is that returned by zero_hold_kep_dyn().

Args:

tol (float): the tolerance of the Taylor adaptive propagator.

Returns:

hy::taylor_adaptive: The Taylor adaptive propagator.

Examples:

>>> import pykep as pk
>>> ta = pk.ta.get_zero_hold_kep(tol = 1e-16)
>>> ta.time = 0.
>>> ta.state[:] = [1.,0.,0.,0.,1.,0.,1.]
>>> mu = 1.
>>> veff = 1.
>>> thrust = [0., 0., 0.]
>>> tof = 1.
>>> ta.pars[:] = [mu, veff] + thrust
>>> ta.propagate_until(tof)

pykep.ta.get_zero_hold_kep_var(tol)

Returns a (order 1) variational Taylor adaptive propagator (Heyoka) for the zero_hold_kep dynamics retreiving one from a global cache and making a copy.

In pykep, using a term well established in electrodynamics, this is the initial value problem of a constant inertial thrust mass-varying spacecraft orbiting a main body.

The dynamics is that returned by zero_hold_kep_dyn(): and also used in get_zero_hold_kep()

Args:

tol (float): the tolerance of the Taylor adaptive propagator.

Returns:

hy::taylor_adaptive: The Taylor adaptive propagator.

Examples:

>>> import pykep as pk
>>> ta = pk.ta.get_zero_hold_kep_var(tol = 1e-16)
>>> ta.time = 0.
>>> ta.state[:] = [1.,0.,0.,0.,1.,0.,1.]
>>> mu = 1.
>>> veff = 1.
>>> thrust = [0., 0., 1e-22]
>>> tof = 1.
>>> ta.pars[:5] = [mu, veff] + thrust
>>> ta.propagate_until(tof)

pykep.ta.zero_hold_kep_dyn()

The dynamics of a constant inertial thrust mass-varying spacecraft orbiting a main body.

\[\begin{split}\left\{ \begin{array}{l} \dot{\mathbf r} = \mathbf v \\ \dot{\mathbf v} = -\frac{\mu}{r^3} \mathbf r + \frac{\mathbf T}{m} \\ \dot m = - \frac{|\mathbf T|}{I_{sp} g_0} \end{array}\right.\end{split}\]

where \(\mu, v_{eff} = I_{sp}g_0\) and \(\mathbf T = [T_x, T_y, T_z]\) are parameters.

Returns:

list [ tuple (hy::expression, hy::expression )]: The dynamics in the form [(x, dx), …]

zero_hold_eq

pykep.ta.get_zero_hold_eq(tol: float)

Returns a Taylor adaptive propagator (Heyoka) for the zero_hold_eq dynamics retreiving one from a global cache and making a copy.

This is the initial value problem of a constant thrust (RTN frame) mass-varying spacecraft orbiting a primary. Thrust direction is fixed in the RTN frame.

If the requested propagator was never created this will create it, else it will return the one from the global cache, thus avoiding jitting.

The dynamics is that returned by zero_hold_eq_dyn().

Args:

tol (float): the tolerance of the Taylor adaptive propagator.

Returns:

hy::taylor_adaptive: The Taylor adaptive propagator.

Examples: >>> import pykep as pk >>> ta = pk.ta.get_zero_hold_eq(tol = 1e-16) >>> ta.time = 0. >>> ta.state[:] = [1.2,0.1,0.1,0.1,1.,0.123,1.] >>> mu = 1.122 >>> veff = 1.32 >>> thrust = [0.23, 0.21, 0.1] >>> tof = 1.23 >>> ta.pars[:] = [mu, veff] + thrust >>> ta.propagate_until(tof)

pykep.ta.get_zero_hold_eq_var(tol: float)

Returns a (order 1) variational Taylor adaptive propagator (Heyoka) for the zero_hold_eq dynamics retreiving one from a global cache and making a copy.

This is the initial value problem of a constant thrust (RTN frame) mass-varying spacecraft orbiting a primary. Thrust direction is fixed in the RTN frame.

Variations are only considered with repsect to initial conditions and the constant thurst \(T_r, T_t, T_n\). The variations with respect to the constant thrust differ from those with respect to the often introduced throttles \(\mathbf T = T_{max} \mathbf u\) by a factor \(T_{max}\).

The dynamics is that returned by zero_hold_eq_dyn(): and also used in get_zero_hold_qe()

Args:

tol (float): the tolerance of the Taylor adaptive propagator.

Returns:

hy::taylor_adaptive: The Taylor adaptive propagator.

Examples: >>> import pykep as pk >>> ta = pk.ta.get_zero_hold_eq_var(tol = 1e-16) >>> ta.time = 0. >>> ta.state[:] = [1.,0.,0.,0.,1.,0.,1.] >>> mu = 1. >>> veff = 1. >>> thrust = [0., 0., 1e-22] >>> tof = 1. >>> ta.pars[:5] = [mu, veff] + thrust >>> ta.propagate_until(tof)

pykep.ta.zero_hold_eq_dyn()

The dynamics in equinoctial elements of a constant thrust (RTN frame) mass-varying spacecraft orbiting a main body.

We consider the motion of a spacecraft of mass \(m\) with a position \(\mathbf{r}\) and velocity \(\mathbf{v}\) subject only to the Sun’s gravitational attraction in some inertial reference frame. The spacecraft also has an ion thruster with a specific impulse \(I_{sp}\) thrusting as \([T_r, T_t, T_n]\). We describe the spacecraft state via its mass \(m\) and the (prograde) modified equinoctial elements \(\mathbf{x}=\left[p,f,g,h,k,L\right]^T\). The dynamics are given by:

\[\begin{split}\begin{array}{l} \dot p = \frac 1m\sqrt{\frac p\mu} \frac{2p}{w} T_t \\ \dot f = \frac 1m\sqrt{\frac p\mu} \left\{T_r \sin L + \left[ (1+w)\cos L + f \right] \frac{T_t}{w} - (h\sin L-k\cos L)\frac{g\cdot T_n}{w} \right\} \\ \dot g = \frac 1m\sqrt{\frac p\mu} \left\{ - T_r\cos L + \left[ (1+w)\sin L + g \right] \frac{T_t}{w} + (h\sin L-k\cos L)\frac{f\cdot T_n}{w} \right\} \\ \dot h = \frac 1m\sqrt{\frac p\mu} \frac{s^2T_n}{2w}\cos L \\ \dot k = \frac 1m\sqrt{\frac p\mu} \frac{s^2T_n}{2w}\sin L \\ \dot L = \sqrt{\frac p\mu}\left\{\mu\left(\frac wp\right)^2 + \frac 1w\left(h\sin L-k\cos L\right)\frac{ T_n}m \right\} \\ \dot m = - \frac{|\mathbf T|}{v_{eff}} \\ \end{array}\end{split}\]

where, \(w = 1 + f\cos L + g\sin L\), \(s^2 = 1 + h^2 + k^2\) and \([T_r, T_t, T_n]\) are the radial, tangential and normal components of thrust direction. The gravitational parameter is denoted with \(\mu\) while \(v_{eff} = I_{sp}g_0\) is the effective velocity of the thruster. We rewrite the above equations in a more compact form:

\[\begin{split}\left\{ \begin{array}{l} \dot {\mathbf x} = \mathbf B(\mathbf x) \frac{\mathbf{T}}m + \mathbf D(\mathbf x) \\ \dot m = - \frac{|\mathbf T|}{v_{eff}} \\ \end{array} \right.\end{split}\]

where:

\[\sqrt{\frac \mu p} \mathbf B(\mathbf x) = \left[ \begin{array}{ccc} 0 & \frac {2p}w & 0 \ \sin L & [(1+w)\cos L + f]\frac 1w & - \frac gw (h\sin L-k\cos L) \ - \cos L & [(1+w)\sin L + g]\frac 1w & \frac fw (h\sin L-k\cos L) \ 0 & 0 & \frac 1w \frac{s^2}{2}\cos L \ 0 & 0 & \frac 1w \frac{s^2}{2}\sin L \ 0 & 0 & \frac 1w (h\sin L - k\cos L) \ \end{array} \right]\]

and:

\[\mathbf D(\mathbf x) = \left[ \begin{array}{cccccc} 0 & 0 & 0 & 0 & 0 & \sqrt{\frac{\mu}{p^3}}w^2 \end{array} \right]^T\]

and \([\mu, v_{eff}, T_x, T_y, T_z]\) are parameters.

Returns:

list [ tuple (hy::expression, hy::expression )]: The dynamics in the form [(p, dp), …]

zero_hold_cr3bp

pykep.ta.get_zero_hold_cr3bp(tol)

Returns a Taylor adaptive propagator (Heyoka) for the zero_hold_cr3bp dynamics retreiving one from a global cache and making a copy.

In pykep, using a term well established in electrodynamics, this is the initial value problem of a constant thrust mass-varying spacecraft orbiting two main bodies in circular relative motion (the circular restricted three body problem, CR3BP). Thrust direction is fixed in the rotating frame.

If the requested propagator was never created this will create it, else it will return the one from the global cache, thus avoiding jitting.

The dynamics is that returned by zero_hold_cr3bp_dyn().

Args:

tol (float): the tolerance of the Taylor adaptive propagator.

Returns:

hy::taylor_adaptive: The Taylor adaptive propagator.

Examples:

>>> import pykep as pk
>>> ta = pk.ta.get_zero_hold_cr3bp(tol = 1e-16)
>>> ta.time = 0.
>>> ta.state[:] = [1.,0.,0.,0.,1.,0.,1.]
>>> mu = 1.
>>> veff = 1.
>>> thrust = [0., 0., 0.]
>>> tof = 1.
>>> ta.pars[:] = [mu, veff] + thrust
>>> ta.propagate_until(tof)

pykep.ta.get_zero_hold_cr3bp_var(tol)

Returns a (order 1) variational Taylor adaptive propagator (Heyoka) for the zero_hold_cr3bp dynamics retreiving one from a global cache and making a copy.

In pykep, using a term well established in electrodynamics, this is the initial value problem of a constant thrust mass-varying spacecraft orbiting two main bodies in circular relative motion (the circular restricted three body problem, CR3BP). Thrust direction is fixed in the rotating frame.

The dynamics is that returned by zero_hold_cr3bp_dyn(): and also used in get_zero_hold_cr3bp()

Args:

tol (float): the tolerance of the Taylor adaptive propagator.

Returns:

hy::taylor_adaptive: The Taylor adaptive propagator.

Examples:

>>> import pykep as pk
>>> ta = pk.ta.get_zero_hold_cr3bp_var(tol = 1e-16)
>>> ta.time = 0.
>>> ta.state[:] = [1.,0.,0.,0.,1.,0.,1.]
>>> mu = 1.
>>> veff = 1.
>>> thrust = [0., 0., 1e-22]
>>> tof = 1.
>>> ta.pars[:5] = [mu, veff] + thrust
>>> ta.propagate_until(tof)

pykep.ta.zero_hold_cr3bp_dyn()

The dynamics of a fixed thrust mass-varying spacecraft orbiting two main bodies in circular relative motion (the circular restricted three body problem, CR3BP).

\[\begin{split}\left\{ \begin{array}{l} \dot{\mathbf r} = \mathbf v \\ \dot v_x = 2v_y + x - (1 - \mu) \frac{x + \mu}{r_1^3} - \mu \frac{x + \mu - 1}{r_2^3} + \frac{T_x}{m}\\ \dot v_y = -2 v_x + y - (1 - \mu) \frac{y}{r_1^3} - \mu \frac{y}{r_2^3} + \frac{T_y}{m}\\ \dot v_z = -(1 - \mu) \frac{z}{r_1^3} - \mu \frac{z}{r_2^3} + \frac{T_z}{m}\\ \dot m = - \frac{|\mathbf T|}{I_{sp} g_0} \end{array}\right.\end{split}\]

where \(\mu, v_{eff} = I_{sp}g_0\) and \(\mathbf T = [T_x, T_y, T_z]\) are parameters.

Returns:

list [ tuple (hy::expression, hy::expression )]: The dynamics in the form [(x, dx), …]

Low-thrust Pontryagin Cartesian TPBVP

pykep.ta.get_pc(tol, optimality)

Returns a Taylor adaptive propagator (Heyoka) for the TPBVP problem resulting from the application of Pontryagin Maximum Principle (PMC) to the low-trhust problem (constant maximal thrust) in Cartesian coordinates. If the requested propagator was never created, a call to this function will trigger its compilation. Otherwise, it will return the one from a global cache, thus avoiding jitting. The specific dynamics used is that returned by pc_dyn(). Both time optimal and mass optimal systems can be returned by setting the optimality parameter.

Args:

tol (float): the tolerance of the Taylor adaptive propagator.

optimality (pykep.optimality_type): the optimality principle to be used.

Returns:

hy::taylor_adaptive: The Taylor adaptive propagator.

Examples:

>>> import pykep as pk
>>> ta = pk.ta.get_pc(tol = 1e-16, optimality = pk.optimality_type.TIME)
>>> ta.time = 0.
>>> # We set the initial conditions with some arbitrary values (all costates to 1.)
>>> ta.state[:14] = [1., 0., 0., 0., 1., 0., 10., 1., 1., 1., 1., 1., 1., 1.]
>>> ta.pars[:] = [1., 0.01, 1.] # in case of TIME parameters are [mu, Tmax, Isp g0]
>>> tof = 1.2345
>>> ta.propagate_until(tof)

pykep.ta.get_pc_var(tol, optimality)

Returns a (order 1) variational Taylor adaptive propagator (Heyoka) for the TPBVP problem resulting from the application of Pontryagin Maximum Principle (PMP) to the low-thrust problem (constant maximal thrust) in Cartesian coordinates. If the requested propagator was never created, a call to this function will trigger its compilation. Otherwise, it will return the one from a global cache, thus avoiding jitting.

The specific dynamics used is that returned by pc_dyn().

Args:

tol (float): the tolerance of the Taylor adaptive propagator.

optimality (pykep.optimality_type): the optimality principle to be used.

Returns:

hy::taylor_adaptive: The Taylor adaptive propagator.

Examples:

>>> import pykep as pk
>>> ta_var = pk.ta.get_pc_var(tol = 1e-16, optimality = pk.optimality_type.TIME))
>>> ta_var.time = 0.
>>> ta_var.state[:14] = [1., 0., 0., 0., 1., 0., 10., 1., 1., 1., 1., 1., 1., 1.]
>>> ta_var.pars[:] = [1., 0.01, 1.] # in case of TIME parameters are [mu, Tmax, Isp g0]
>>> tof = 1.2345
>>> ta_var.propagate_until(tof)

pykep.ta.pc_dyn(optimality)

The augmented dynamics of the TPBVP originating when applying an indirect method to the low-thrust transfer problem in Cartesian coordinates.

The (non-augmented) dynamics is,

\[\begin{split}\left\{ \begin{array}{l} \dot{\mathbf r} = \mathbf f_r = \mathbf v \\ \dot{\mathbf v} = \mathbf f_v = -\frac{mu}{r^3}\mathbf r + c_1 \frac um \hat{\mathbf i}\\ \dot{m} = f_m = - \frac{c_1}{c_2} u \end{array} \right.\end{split}\]

The state, containing the spacecraft state and the co-states, is:

\[\mathbf x = [\mathbf r, \mathbf v, m] = [x,y,z,v_x,v_y,v_z,l_x,l_y,l_z,l_{vx},l_{vy},l_{vz},l_m]\]

While the parameters:

\[\mathbf p = [\mu, c_1, c_2, \epsilon, \lambda_0]\]

describe, respectively, the gravitational parameter, the maximum thrust, the effective velocity (product of the specific impulse \(I_{sp}\) by \(g_0\)), an homotopy parameter and a factor multiplying the instantaneous cost (in theory this can be any positive number as the solution of the problem will not change)

The controls, representing magnitude and direction of the spacecraft thrust are:

\[\mathbf u = [u, \hat{\mathbf i}]\]

The final equations of motion (in the case of mass optimality), are derived from the Hamiltonian:

\[\mathcal H(\mathbf x, \mathbf \lambda, \mathbf u) = \mathbf \lambda_r \cdot \mathbf f_r + \mathbf \lambda_v \cdot \mathbf f_v + \lambda_m f_m + \lambda_0 \frac{c_1}{c_2} \left(u + \epsilon\log(u(1-u))\right)\]

or, in the case of time optimality:

\[\mathcal H(\mathbf x, \mathbf \lambda, \mathbf u) = \mathbf \lambda_r \cdot \mathbf f_r + \mathbf \lambda_v \cdot \mathbf f_v + \lambda_m f_m + \lambda_0 \frac{c_1}{c_2}\]

by taking the derivatives:

\[\begin{split}\left\{ \begin{array}{l} \dot{\mathbf x} = \frac{\partial \mathcal H}{\partial \mathbf \lambda} \\ \dot{\mathbf \lambda} = - \frac{\partial \mathcal H}{\partial \mathbf x} \\ \end{array}\right.\end{split}\]

The equation are then made independent of the controls applying Pontryagin minimum principle which takes the general form:

\[\mathbf u^* = \textbf{argmin}_{\mathbf u \in \mathcal U} \mathcal H(\mathbf x, \mathbf \lambda, \mathbf u)\]

where \(\mathcal U\) is the set of admissible controls, i.e. the set of thrust directions and magnitudes.

and substituting in the equations the optimal controls found as function of the augmented system state.

Note

In the case of time optimality, the dynamics will not depend on neither \(\epsilon\) nor \(\lambda_0\) and the size of the parameter vector will this be three.

Args:

optimality (pykep.optimality_type): the optimality principle to be used.

Returns:

list [ tuple (hy::expression, hy::expression )]: The dynamics in the form [(x, dx), …]

Low-thrust Pontryagin Equinoctial TPBVP

pykep.ta.get_peq(tol, optimality)

Returns a Taylor adaptive propagator (Heyoka) for the TPBVP problem resulting from the application of Pontryagin Maximum Principle (PMC) to the low-trhust problem (constant maximal thrust) in Modified Equinoctial elements. If the requested propagator was never created, a call to this function will trigger its compilation. Otherwise, it will return the one from a global cache, thus avoiding jitting. The specific dynamics used is that returned by peq_dyn(). Both time optimal and mass optimal systems can be returned by setting the optimality parameter.

Args:

tol (float): the tolerance of the Taylor adaptive propagator.

optimality (pykep.optimality_type): the optimality principle to be used.

Returns:

hy::taylor_adaptive: The Taylor adaptive propagator.

Examples:

>>> import pykep as pk
>>> ta = pk.ta.get_peq(tol = 1e-16, optimality = pk.optimality_type.TIME)
>>> ta.time = 0.
>>> # We set the initial conditions with some arbitrary values (all costates to 1.)
>>> ta.state[:14] = [1., 0., 0., 0., 1., 0., 10., 1., 1., 1., 1., 1., 1., 1.]
>>> ta.pars[:] = [1., 0.01, 1.] # in case of TIME parameters are [mu, Tmax, Isp g0]
>>> tof = 1.2345
>>> ta.propagate_until(tof)

pykep.ta.get_peq_var(tol, optimality)

Returns a (order 1) variational Taylor adaptive propagator (Heyoka) for the TPBVP problem resulting from the application of Pontryagin Maximum Principle (PMP) to the low-thrust problem (constant maximal thrust) in Modified Equinoctial Elements. If the requested propagator was never created, a call to this function will trigger its compilation. Otherwise, it will return the one from a global cache, thus avoiding jitting.

The specific dynamics used is that returned by pc_dyn().

Args:

tol (float): the tolerance of the Taylor adaptive propagator.

optimality (pykep.optimality_type): the optimality principle to be used.

Returns:

hy::taylor_adaptive: The Taylor adaptive propagator.

Examples:

>>> import pykep as pk
>>> ta_var = pk.ta.get_peq_var(tol = 1e-16, optimality = pk.optimality_type.TIME))
>>> ta_var.time = 0.
>>> ta_var.state[:14] = [1., 0., 0., 0., 1., 0., 10., 1., 1., 1., 1., 1., 1., 1.]
>>> ta_var.pars[:] = [1., 0.01, 1.] # in case of TIME parameters are [mu, Tmax, Isp g0]
>>> tof = 1.2345
>>> ta_var.propagate_until(tof)

pykep.ta.peq_dyn(optimality)

The augmented dynamics of the TPBVP originating when applying an indirect method to the low-thrust transfer problem in Modified Equinoctial Elements (MEE).

The (non-augmented) dynamics is:

\[\begin{split}\left\{ \begin{array}{l} \dot p = \sqrt{\frac{p}{\mu}} \frac{2p}{w} i_t \cdot \frac{u}{m} \\ \dot f = \frac{1}{m} \sqrt{\frac{p}{\mu}} \left[ i_r \sin L + \left( (1+w)\cos L + f \right)\frac{i_t}{w} - (h\sin L - k\cos L)\frac{g i_n}{w} \right] \cdot \frac{u}{m} \\ \dot g = \frac{1}{m} \sqrt{\frac{p}{\mu}} \left[ -i_r \cos L + \left( (1+w)\sin L + g \right)\frac{i_t}{w} + (h\sin L - k\cos L)\frac{f i_n}{w} \right] \cdot \frac{u}{m} \\ \dot h = \sqrt{\frac{p}{\mu}} \frac{s^2 i_n}{2w} \cos L \cdot \frac{u}{m} \\ \dot k = \sqrt{\frac{p}{\mu}} \frac{s^2 i_n}{2w} \sin L \cdot \frac{u}{m} \\ \dot L = \sqrt{\frac{p}{\mu}} \left[ \mu\left(\frac{w}{p}\right)^2 + \frac{1}{w}(h\sin L - k\cos L) \frac{i_n}{m} \cdot \frac{u}{m} \right] \\ \dot m = -\frac{c_1}{c_2} \cdot \frac{u}{m} \end{array} \right.\end{split}\]

The auxiliary functions and variables used above are:

\[w = 1 + f \cos L + g \sin L,\quad s^2 = 1 + h^2 + k^2\]

Introducing the state, containing both the physical state and the co-states, as:

\[\mathbf x = [p, f, g, h, k, L, m, \lambda_p, \lambda_f, \lambda_g, \lambda_h, \lambda_k, \lambda_L, \lambda_m]\]

the parameter vector, containing the gravitational parameter, maximum thrust, effective exhaust velocity, a homotopy continuation parameter, and a normalizing Lagrange multiplier (see later), as:

\[\mathbf p = [\mu, c_1, c_2, \epsilon, \lambda_0]\]

and the control vector as:

\[\mathbf u = [u, \hat{\mathbf i}_\tau]\]

the dynamics can be rewritten in compact form as:

\[\begin{split}\left\{ \begin{array}{l} \dot{\mathbf x} = \frac{c_1 u(t)}{m} \mathbf B(\mathbf x) \hat{\mathbf i}_\tau + \mathbf D(\mathbf x) \\ \dot m = -c_2 u(t) \end{array} \right.\end{split}\]

Where:

\[\begin{split}\sqrt{\frac{\mu}{p}} \mathbf B(\mathbf x) = \begin{bmatrix} 0 & \frac{2p}{w} & 0 \\ \sin L & \frac{(1+w)\cos L + f}{w} & -\frac{g}{w}(h \sin L - k \cos L) \\ -\cos L & \frac{(1+w)\sin L + g}{w} & \frac{f}{w}(h \sin L - k \cos L) \\ 0 & 0 & \frac{1}{w} \frac{s^2}{2} \cos L \\ 0 & 0 & \frac{1}{w} \frac{s^2}{2} \sin L \\ 0 & 0 & \frac{1}{w} (h \sin L - k \cos L) \end{bmatrix}\end{split}\]

\[\begin{split}\mathbf D(\mathbf x) = \begin{bmatrix} 0 \\ 0 \\ 0 \\ 0 \\ 0 \\ \sqrt{\frac{\mu}{p^3}} w^2 \end{bmatrix}\end{split}\]

As in the Cartesian case, the equations of motion for the TPBVP are derived from the Hamiltonian, which depends on the chosen optimality:

• Mass Optimality:

\[\mathcal H = \boldsymbol \lambda^\top \dot{\mathbf x} + \lambda_m \dot{m} + \lambda_0 \frac{c_1}{c_2} \left( u + \epsilon \log(u(1 - u)) \right)\]

• Time Optimality:

\[\mathcal H = \boldsymbol \lambda^\top \dot{\mathbf x} + \lambda_m \dot{m} + \lambda_0 \frac{c_1}{c_2}\]

The full augmented dynamics is thus obtained by:

\[\begin{split}\left\{ \begin{array}{l} \dot{\mathbf x} = \frac{\partial \mathcal H}{\partial \boldsymbol \lambda} \\ \dot{\boldsymbol \lambda} = - \frac{\partial \mathcal H}{\partial \mathbf x} \end{array} \right.\end{split}\]

And by applying Pontryagin’s Minimum Principle:

\[\mathbf u^* = \arg\min_{\mathbf u \in \mathcal U} \mathcal H(\mathbf x, \boldsymbol \lambda, \mathbf u)\]

where \(\mathcal U\) is the admissible set of control inputs (direction and throttle).

Args:

optimality (pykep.optimality_type): the optimality principle to be used.

Returns:

list [ tuple (hy::expression, hy::expression )]: The dynamics in the form [(x, dx), …]