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]

---

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

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

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

pykep.ta.get_zoh_kep(tol: float)

Returns a Taylor adaptive propagator (Heyoka) for the zoh_kep_dyn() dynamics retrieving one from a global cache if available.

This solves the initial value problem of a constant thrust mass-varying spacecraft. Thrust direction is fixed in the inertial frame.

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
veff = 1.32
controls = [0.022, 0.023, -0.21, 0.1]
tof = 1.23

Create propagator and propagate:

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

pykep.ta.get_zoh_kep_var(tol: float)

Returns a (order 1) variational Taylor adaptive propagator (Heyoka) for the zoh_kep_dyn() dynamics retrieving one from a global cache if available.

Note

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

Args:

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

Returns:

hy::taylor_adaptive: The Taylor adaptive variational propagator.

Examples:

Import and setup:

import pykep as pk
veff = 1.32
controls = [0.022, 0.023, -0.21, 0.1]
tof = 1.23

Create propagator and propagate:

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

pykep.ta.zoh_kep_dyn()

The dynamics in Cartesian coordinates of a constant thrust mass-varying spacecraft orbiting a main central body with unitary gravitational parameter \(\mu = 1\). The electric propulsion is characterized by its effective exhaust velocity \(v_{eff} = I_{sp} g_0\).

The dynamics are given by:

\[\begin{split}\left\{ \begin{array}{l} \dot{\mathbf{r}} = \mathbf{v} \\ \dot{\mathbf{v}} = -\frac{\mathbf{r}}{r^3} + \frac{T}{m} \mathbf{\hat{i}} \\ \dot{m} = -c T \end{array} \right.\end{split}\]

where \(c = \frac{1}{v_{eff}}\), \(\mathbf{\hat{i}} = [i_x, i_y, i_z]\).

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 [(p, dp), …]

Zero Hold Keplerian Propagator in Equinoctial elements

pykep.ta.get_zoh_eq(tol: float)

Returns a Taylor adaptive propagator (Heyoka) for the zoh_eq_dyn() dynamics retrieving one from a global cache if available.

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

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
veff = 1.32
controls = [0.022, 0.023, -0.21, 0.1]
tof = 1.23

Create propagator and propagate:

ta = pk.ta.get_zoh_eq(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(tof)

pykep.ta.get_zoh_eq_var(tol: float)

Returns a (order 1) variational Taylor adaptive propagator (Heyoka) for the zoh_eq_dyn() dynamics retrieving one from a global cache if available.

Note

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

Args:

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

Returns:

hy::taylor_adaptive: The Taylor adaptive variational propagator.

Examples:

Import and setup:

import pykep as pk
veff = 1.32
controls = [0.022, 0.023, -0.21, 0.1]
tof = 1.23

Create propagator and propagate:

ta = pk.ta.get_zoh_eq_var(tol=1e-16)
ta.time = 0.
ta.state[:] = [1., 0., 0., 0., 1., 0., 1.]
ta.pars[:] = controls + [1 / veff]
ta.propagate_until(tof)

pykep.ta.zoh_eq_dyn()

The dynamics in equinoctial elements of a constant thrust (RTN frame) mass-varying spacecraft orbiting a main central body with unitary gravitational parameter \(\mu = 1\). The electric propulsion is characterized by its effective exhaust velocity \(v_{eff} = I_{sp} g_0\).

The dynamics are given by:

\[\begin{split}\left\{ \begin{array}{l} \dot p = \frac 1m\sqrt{p} \frac{2p}{w} T_t \\ \dot f = \frac 1m\sqrt{p} \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{p} \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{p} \frac{s^2T_n}{2w}\cos L \\ \dot k = \frac 1m\sqrt{p} \frac{s^2T_n}{2w}\sin L \\ \dot L = \sqrt{p}\left\{\left(\frac wp\right)^2 + \frac 1w\left(h\sin L-k\cos L\right)\frac{ T_n}m \right\} \\ \dot m = - c T \\ \end{array} \right.\end{split}\]

where \(c = \frac{1}{v_{eff}}\), \(w = 1 + f\cos L + g\sin L\), \(s^2 = 1 + h^2 + k^2\).

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

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

In compact matrix form:

\[\dot{\mathbf{x}} = \mathbf{B}(\mathbf{x}) \frac{\mathbf{T}}{m} + \mathbf{D}(\mathbf{x})\]

with:

\[\begin{split}\sqrt{\frac{p}{1}} \mathbf{B}(\mathbf{x}) = \begin{bmatrix} 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{bmatrix}\end{split}\]

and:

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

Returns:

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

Zero Hold CR3BP Propagator in Cartesian coordinates

pykep.ta.get_zoh_cr3bp(tol: float)

Returns a Taylor adaptive propagator (Heyoka) for the zoh_cr3bp_dyn() dynamics retrieving one from a global cache if available.

This solves the initial value problem of a constant thrust mass-varying spacecraft in the circular restricted three-body problem (CR3BP). Thrust direction is fixed in the rotating frame.

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

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
mu = pk.CR3BP_MU_EARTH_MOON
veff = 1.
controls = [0.001, 1., 0., 0.]
tof = 1.

Create propagator and propagate:

ta = pk.ta.get_zoh_cr3bp(tol=1e-16)
ta.time = 0.
ta.state[:] = [1., 0., 0., 0., 1., 0., 1.]
ta.pars[:] = controls + [1 / veff, mu]
ta.propagate_until(tof)

pykep.ta.get_zoh_cr3bp_var(tol: float)

Returns a (order 1) variational Taylor adaptive propagator (Heyoka) for the zoh_cr3bp_dyn() dynamics retrieving one from a global cache if available.

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 and the thrust parameters \([T, i_x, i_y, i_z]\).

Args:

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

Returns:

hy::taylor_adaptive: The Taylor adaptive variational propagator.

Examples:

Import and setup:

import pykep as pk
mu = pk.CR3BP_MU_EARTH_MOON
veff = 1.
controls = [0.001, 1., 0., 0.]
tof = 1.

Create propagator and propagate:

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, mu]
ta.propagate_until(tof)

pykep.ta.zoh_cr3bp_dyn()

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

The dynamics are given by:

\[\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} + T\frac{i_x}{m} \\ \dot{v}_y = -2v_x + y - (1-\mu)\frac{y}{r_1^3} - \mu\frac{y}{r_2^3} + T\frac{i_y}{m} \\ \dot{v}_z = -(1-\mu)\frac{z}{r_1^3} - \mu\frac{z}{r_2^3} + T\frac{i_z}{m} \\ \dot{m} = -c T \end{array} \right.\end{split}\]

where:

• \(r_1 = ||\mathbf{r} - \mathbf{r}_1||\), \(r_2 = ||\mathbf{r} - \mathbf{r}_2||\)

• \(\mathbf{r}_1 = [-\mu, 0, 0]\), \(\mathbf{r}_2 = [1-\mu, 0, 0]\)

• \(c = \frac{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 form [(x, dx), …]

Zero Hold Solar Sail Propagator in Cartesian coordinates

pykep.ta.get_zoh_ss(tol: float)

Returns a Taylor adaptive propagator (Heyoka) for the zoh_ss_dyn() dynamics retrieving one from a global cache if available.

This solves the initial value problem of a solar sailing spacecraft under simple hypothesis for the sail characteristics.

Sail acceleration direction is fixed in the RTN frame.

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
sail_c = 0.01
controls = [0.022, 0.023]
tof = 1.23

Create propagator and propagate:

ta = pk.ta.get_zoh_ss(tol=1e-16)
ta.time = 0.
ta.state[:] = [1.2, 0.1, 0.1, 0.1, 1., 0.123, 1.]
ta.pars[:] = controls + [sail_c]
ta.propagate_until(tof)

pykep.ta.get_zoh_ss_var(tol: float)

Returns a (order 1) variational Taylor adaptive propagator (Heyoka) for the zoh_ss_dyn() dynamics retrieving one from a global cache if available.

Note

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

Args:

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

Returns:

hy::taylor_adaptive: The Taylor adaptive variational propagator.

Examples:

Import and setup:

import pykep as pk
sail_c = 0.01
controls = [0.022, 0.023]
tof = 1.23

Create propagator and propagate:

ta = pk.ta.get_zoh_ss_var(tol=1e-16)
ta.time = 0.
ta.state[:7] = [1.2, 0.1, 0.1, 0.1, 1., 0.123]
ta.pars[:] = controls + [sail_c]
ta.propagate_until(tof)

pykep.ta.zoh_ss_dyn()

The dynamics in Cartesian coordinates of a spacecraft equipped with a solar sail orbiting a main central body with unitary gravitational parameter \(\mu = 1\). The body is assumed to be a star with flux \(C\) at a reference distance R = 1. These conventions define naturally the nondimenional units of the system, which in the solar system, for example, are \(MU = \mu_{Sun}\) and \(L = 1 AU\) (so that C is the flux at 1AU).

The sail attitude is controlled through two angles (\(\alpha\), \(\beta\)), defining the direction of the sail acceleration: they are commonly named cone and clock angle. The angles are assumed to be fixed, hence creating an acceleration having constant direction in the RTN reference frame \([\mathbf i_R, \mathbf i_T, \mathbf i_N]\).

The dynamics are given in non dimensional units (:matby:

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

where the acceleration coming from the solar sail is:

\[\mathbf a_{ss} = T \left( \cos\alpha \mathbf i_R + \sin\alpha\sin\beta \mathbf i_T + \sin\alpha\cos\beta \mathbf i_N\right)\]

and:

\[T = c (1. / r^2) \cos^2\alpha\]

The constant \(c\) (in acceleration units) contains all physical parameters defining the sailcraft and is defined as \(c = 2 \frac{CA}M\), where \(C\) is central star flux at the reference unitary distance, A the sail area and M its mass.

The cone angle must be constrained in \(\alpha \in [-\frac{\pi}2, \frac{\pi}2]\) for the dynamics to be feasible (outside these bound the corresponding sail acceleration is unfeasible)

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

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

Note

with respect to other Zero Order Hold integrators the ``zoh_ss``has a smaller dimension (no mass equation) and only two controls. It CANNOT thus be used as Taylor integrator in zoh() as it does not meet its requirements. A dedicated leg, zoh_ss() is instead to be used.

Returns:

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

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