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 = np.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 convenience 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]].

tofs (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 = [np.pi/2, np.pi, 3*np.pi/4]
>>> 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]

---

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 of 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. These are intended to be used to build “zero hold” trajectory legs and direct methods to optimally control the spacecraft.

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 retrieving 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 retrieving 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.

Note

Variations are only considered with respect to initial conditions.

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 retrieving 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 retrieving 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.

Note

Variations are only considered with respect to initial conditions.

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(tol)

Returns a Taylor adaptive propagator (Heyoka) for the Bicircular Problem (BCP) dynamics retrieving one from a global cache (making a copy).

If the requested propagator was never created, this creates it; otherwise returns cached version (avoiding re-jitting).

The BCP is defined in Cartesian coordinates (non-symplectic, non-Hamiltonian, time-dependent) with Sun on the \(x\) axis at \(t=0\).

Args:

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

Returns:

hy::taylor_adaptive: The Taylor adaptive propagator.

Examples:

Import and setup:

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_EARTH_MOON
mu_s = pk.BCP_MU_S
rho_s = pk.BCP_RHO_S
rho_p = pk.BCP_RHO_S
tof = 5.7856656782589234

Create propagator and propagate:

ta.pars[:] = [mu, mu_s, rho_s, rho_p]
ta.propagate_until(tof)

ta.get_bcp(tol)

Returns a Taylor adaptive propagator (Heyoka) for the Bicircular Problem (BCP) dynamics retrieving one from a global cache (making a copy).

If the requested propagator was never created, this creates it; otherwise returns cached version (avoiding re-jitting).

The BCP is defined in Cartesian coordinates (non-symplectic, non-Hamiltonian, time-dependent) with Sun on the \(x\) axis at \(t=0\).

Args:

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

Returns:

hy::taylor_adaptive: The Taylor adaptive propagator.

Examples:

Import and setup:

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_EARTH_MOON
mu_s = pk.BCP_MU_S
rho_s = pk.BCP_RHO_S
rho_p = pk.BCP_RHO_S
tof = 5.7856656782589234

Create propagator and propagate:

ta.pars[:] = [mu, mu_s, rho_s, rho_p]
ta.propagate_until(tof)

pykep.ta.get_bcp_var(tol: float)

Returns a (order 1) variational Taylor adaptive propagator (Heyoka) for the Bicircular Problem (BCP) dynamics retrieving one from a global cache (making a copy).

If the requested propagator was never created, this creates it; otherwise returns cached version (avoiding re-jitting).

Note

Variations are only considered with respect to initial conditions.

The BCP is defined in Cartesian coordinates (non-symplectic, non-Hamiltonian, time-dependent) with Sun on the \(x\) axis at \(t=0\).

Args:

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

Returns:

hy::taylor_adaptive: The variational Taylor adaptive propagator.

Examples:

Import and setup:

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_EARTH_MOON
mu_s = pk.BCP_MU_S
rho_s = pk.BCP_RHO_S
rho_p = pk.BCP_RHO_S
tof = 5.7856656782589234

Create propagator and propagate:

ta.pars[:] = [mu, mu_s, rho_s, rho_p]
ta.propagate_until(tof)

ta.get_bcp_var(tol:float)

Returns a (order 1) variational Taylor adaptive propagator (Heyoka) for the Bicircular Problem (BCP) dynamics retrieving one from a global cache (making a copy).

If the requested propagator was never created, this creates it; otherwise returns cached version (avoiding re-jitting).

Note

Variations are only considered with respect to initial conditions.

The BCP is defined in Cartesian coordinates (non-symplectic, non-Hamiltonian, time-dependent) with Sun on the \(x\) axis at \(t=0\).

Args:

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

Returns:

hy::taylor_adaptive: The variational Taylor adaptive propagator.

Examples:

Import and setup:

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_EARTH_MOON
mu_s = pk.BCP_MU_S
rho_s = pk.BCP_RHO_S
rho_p = pk.BCP_RHO_S
tof = 5.7856656782589234

Create propagator and propagate:

ta.pars[:] = [mu, mu_s, rho_s, rho_p]
ta.propagate_until(tof)

pykep.ta.bcp_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 spacecraft 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), …]

bcp_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 spacecraft 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 Keplerian Propagator in Cartesian coordinates

pykep.ta.get_zoh_kep(tol)

Returns a Taylor adaptive propagator (Heyoka) for the zero-order-hold Keplerian low-thrust dynamics, retrieving one from a global cache and making a copy.

If the requested propagator was never created, this creates it; otherwise returns a copy of the cached version (avoiding re-jitting while keeping independent runtime state).

The specific dynamics used is that returned by zoh_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
>>> veff = 1.32
>>> controls = [0.022, 0.023, -0.21, 0.1]
>>> ta = pk.ta.get_zoh_kep(tol = 1e-16)
>>> ta.time = 0.
>>> ta.state[:] = [1.2, 0.1, 0.1, 0.1, 1., 0.123, 1.]
>>> ta.pars[:] = controls + [1. / veff]
>>> ta.propagate_until(1.23)

pykep.ta.get_zoh_kep_var(tol)

Returns a (order 1) variational Taylor adaptive propagator (Heyoka) for the zero-order-hold Keplerian low-thrust dynamics, retrieving one from a global cache and making a copy.

If the requested propagator was never created, this creates it; otherwise returns a copy of the cached version (avoiding re-jitting while keeping independent runtime state).

Note

Variations are only considered with respect to initial conditions and the thrust parameters \([T, i_x, i_y, i_z]\).

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

Args:

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

Returns:

hy::taylor_adaptive: The variational Taylor adaptive propagator.

Examples:

>>> import pykep as pk
>>> veff = 1.32
>>> controls = [0.022, 0.023, -0.21, 0.1]
>>> ta = pk.ta.get_zoh_kep_var(tol = 1e-16)
>>> ta.time = 0.
>>> ta.state[:7] = [1.2, 0.1, 0.1, 0.1, 1., 0.123, 1.]
>>> ta.pars[:] = controls + [1. / veff]
>>> ta.propagate_until(1.23)

pykep.ta.zoh_kep_dyn()

The dynamics in Cartesian coordinates of a constant-thrust mass-varying spacecraft orbiting a central body with unitary gravitational parameter \(\mu = 1\).

The electric propulsion is characterized by its effective exhaust velocity \(v_{eff} = I_{sp} g_0\), with \(c = 1 / v_{eff}\).

The state is: \([x, y, z, v_x, v_y, v_z, m]\)

The system parameters are (in this order): \([T, i_x, i_y, i_z] + [c]\).

Returns:

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

Zero Hold Keplerian Propagator in Equinoctial elements

pykep.ta.get_zoh_eq(tol)

Returns a Taylor adaptive propagator (Heyoka) for the zero-order-hold equinoctial low-thrust dynamics, retrieving one from a global cache and making a copy.

If the requested propagator was never created, this creates it; otherwise returns a copy of the cached version (avoiding re-jitting while keeping independent runtime state).

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

Args:

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

Returns:

hy::taylor_adaptive: The Taylor adaptive propagator.

pykep.ta.get_zoh_eq_var(tol)

Returns a (order 1) variational Taylor adaptive propagator (Heyoka) for the zero-order-hold Equinoctial low-thrust dynamics, retrieving one from a global cache and making a copy.

If the requested propagator was never created, this creates it; otherwise returns a copy of the cached version (avoiding re-jitting while keeping independent runtime state).

Note

Variations are only considered with respect to initial conditions and the thrust parameters \([T, i_r, i_t, i_n]\).

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

Args:

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

Returns:

hy::taylor_adaptive: The variational Taylor adaptive propagator.

pykep.ta.zoh_eq_dyn()

The dynamics in equinoctial elements of a constant-thrust mass-varying spacecraft orbiting a central body with unitary gravitational parameter \(\mu = 1\).

The electric propulsion is characterized by its effective exhaust velocity \(v_{eff} = I_{sp} g_0\), with \(c = 1 / v_{eff}\).

The state is: \([p, f, g, h, k, L, m]\)

The system parameters are (in this order): \([T, i_r, i_t, i_n] + [c]\).

Returns:

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

Zero Hold CR3BP Propagator in Cartesian coordinates

pykep.ta.get_zoh_cr3bp(tol)

Returns a Taylor adaptive propagator (Heyoka) for the zero-order-hold CR3BP low-thrust dynamics, retrieving one from a global cache and making a copy.

If the requested propagator was never created, this creates it; otherwise returns a copy of the cached version (avoiding re-jitting while keeping independent runtime state).

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

Args:

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

Returns:

hy::taylor_adaptive: The Taylor adaptive propagator.

pykep.ta.get_zoh_cr3bp_var(tol)

Returns a (order 1) variational Taylor adaptive propagator (Heyoka) for the zero-order-hold CR3BP low-thrust dynamics, retrieving one from a global cache and making a copy.

If the requested propagator was never created, this creates it; otherwise returns a copy of the cached version (avoiding re-jitting while keeping independent runtime state).

Note

Variations are only considered with respect to initial conditions and the thrust parameters \([T, i_x, i_y, i_z]\).

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

Args:

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

Returns:

hy::taylor_adaptive: The variational Taylor adaptive propagator.

pykep.ta.zoh_cr3bp_dyn()

The dynamics of a fixed-thrust mass-varying spacecraft in the circular restricted three-body problem (CR3BP), with thrust fixed in the rotating frame.

The state is: \([x, y, z, v_x, v_y, v_z, m]\)

The system parameters are (in this order): \([T, i_x, i_y, i_z] + [c, \mu]\).

Returns:

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

Zero Hold Solar Sail Propagator in Cartesian coordinates

pykep.ta.get_zoh_ss(tol)

Returns a Taylor adaptive propagator (Heyoka) for the zero-order-hold solar-sail dynamics, retrieving one from a global cache and making a copy.

If the requested propagator was never created, this creates it; otherwise returns a copy of the cached version (avoiding re-jitting while keeping independent runtime state).

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

Args:

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

Returns:

hy::taylor_adaptive: The Taylor adaptive propagator.

pykep.ta.get_zoh_ss_var(tol)

Returns a (order 1) variational Taylor adaptive propagator (Heyoka) for the zero-order-hold solar-sail dynamics, retrieving one from a global cache and making a copy.

If the requested propagator was never created, this creates it; otherwise returns a copy of the cached version (avoiding re-jitting while keeping independent runtime state).

Note

Variations are only considered with respect to initial conditions and the sail-angle parameters \([\alpha, \beta]\).

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

Args:

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

Returns:

hy::taylor_adaptive: The variational Taylor adaptive propagator.

pykep.ta.zoh_ss_dyn()

The dynamics in Cartesian coordinates of a solar-sail spacecraft around a central body with unitary gravitational parameter \(\mu = 1\).

The state is: \([x, y, z, v_x, v_y, v_z]\)

The system parameters are (in this order): \([\alpha, \beta] + [c]\).

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.

Note

Variations are considered with respect to the initial conditions on the costates and to the parameters \(\epsilon\) and \(\lambda_0\). In the time optimal case \(\epsilon\) is still considered a parameter (for consistency) but it is not used.

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.

Note

Variations are considered with respect to the initial conditions on the costates and to the parameters \(\epsilon\) and \(\lambda_0\). In the time optimal case \(\epsilon\) is still considered a parameter (for consistency) but it is not used.

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), …]