poliastro: a Python library for interactive astrodynamics

Introduction¶

History¶

The term “astrodynamics” was coined by the American astronomer Samuel Herrick, who received encouragement from the space pioneer Robert H. Goddard, and refers to the branch of space science dealing with the motion of artificial celestial bodies Duboshin, 1973Herrick, 1971. However, the roots of its mathematical foundations go back several centuries.

Kepler first introduced his laws of planetary motion in 1609 and 1619 and derived his famous transcendental equation (1), which we now see as capturing a restricted form of the two-body problem. This work was generalized by Newton to give birth to the n-body problem, and many other mathematicians worked on it throughout the centuries (Daniel and Johann Bernoulli, Euler, Gauss). Poincaré established in the 1890s that no general closed-form solution exists for the n-body problem, since the resulting dynamical system is chaotic Battin, 1999. Sundman proved in the 1900s the existence of convergent solutions for a few restricted with $n = 3$.

In 1903 Tsiokovsky evaluated the conditions required for artificial objects to leave the orbit of the earth; this is considered as a foundational contribution to the field of astrodynamics. Tsiokovsky devised equation (2) which relates the increase in velocity with the effective exhaust velocity of thrusted gases and the fraction of used propellant.

Further developments by Kondratyuk, Hohmann, and Oberth in the early 20th century all added to the growing field of orbital mechanics, which in turn enabled the development of space flight in the USSR and the United States in the 1950s and 1960s.

The two-body problem¶

In a system of $i \in {1, ..., n}$ bodies subject to their mutual attraction, by application of Newton’s law of universal gravitation, the total force $\vec{f}_i$ affecting $m_i$ due to the presence of the other $n - 1$ masses is given by Battin (1999):

where $G = 6.67430\cdot 10^{-11} N m^{2} kg^{-2}$ is the universal gravitational constant, and $\vec{r}_{ij}$ denotes the position vector from $m_i$ to $m_j$. Applying Newton’s second law of motion results in a system of $n$ differential equations:

By setting $n = 2$ in (4) and subtracting the two resulting equalities, one arrives to the fundamental equation of the two-body problem:

where $\mu = G(m_1 + m_2) = G(M + m)$. When $m \ll M$ (for example, an artificial satellite orbiting a planet), one can consider $\mu = GM$ a property of the attractor.

Keplerian vs non-keplerian motion¶

Conveniently manipulating equation (5) leads to several properties Battin, 1999 that were already published by Johannes Kepler in the 1610s, namely:

1. The orbit always describes a conic section (an ellipse, a parabola, or an hyperbola), with the attractor at one of the two foci and can be written in polar coordinates like $r = \frac{p}{1 + e \cos{\nu}}$ (Kepler’s first law).
2. The magnitude of the specific angular momentum $h = r^2 \frac{\diff{\theta}}{\diff t}$ is constant an equal to two times the areal velocity (Kepler’s second law).
3. For closed (circular and elliptical) orbits, the period is related to the size of the orbit through $P = 2 \pi \sqrt{\frac{a^3}{\mu}}$ (Kepler’s third law).

For many practical purposes it is usually sufficient to limit the study to one object orbiting an attractor and ignore all other external forces of the system, hence restricting the study to trajectories governed by equation (5). Such trajectories are called “Keplerian”, and several problems can be formulated for them:

• The initial-value problem, which is usually called propagation, involves determining the position and velocity of an object after an elapse period of time given some initial conditions.
• Preliminary orbit determination, which involves using exact or approximate methods to derive a Keplerian orbit from a set of observations.
• The boundary-value problem, often named the Lambert problem, which involves determining a Keplerian orbit from boundary conditions, usually departure and arrival position vectors and a time of flight.

Fortunately, most of these problems boil down to finding numerical solutions to relatively simple algebraic relations between time and angular variables: for elliptic motion ($0 \le e < 1$) it is the Kepler equation, and equivalent relations exist for the other eccentricity regimes Battin, 1999. Numerical solutions for these equations can be found in a number of different ways, each one with different complexity and precision tradeoffs. In the Methods section we list the ones implemented by poliastro.

On the other hand, there are many situations in which natural and artificial orbital perturbations must be taken into account so that the actual non-Keplerian motion can be properly analyzed:

• Interplanetary travel in the proximity of other planets. On a first approximation it is usually enough to study the trajectory in segments and focus the analysis on the closest attractor, hence patching several Keplerian orbits along the way (the so-called “patched-conic approximation”) Battin, 1999. The boundary surface that separates one segment from the other is called the sphere of influence.
• Use of solar sails, electric propulsion, or other means of continuous thrust. Devising the optimal guidance laws that minimize travel time or fuel consumption under these conditions is usually treated as an optimization problem of a dynamical system, and as such it is particularly challenging Conway, 2014.
• Artificial satellites in the vicinity of a planet. This is the regime in which all the commercial space industry operates, especially for those satellites in Low-Earth Orbit (LEO).

State of the art¶

In our view, at the time of creating poliastro there were a number of issues with existing open source astrodynamics software that posed a barrier of entry for novices and amateur practitioners. Most of these barriers still exist today and are described in the following paragraphs. The goals of the project can be condensed as follows:

1. Set an example on reproducibility and good coding practices in astrodynamics.
2. Become an approachable software even for novices.
3. Offer a performant software that can be also used in scripting and interactive workflows.

The most mature software libraries for astrodynamics are arguably Orekit Orekit, 2022, a “low level space dynamics library written in Java” with an open governance model, and SPICE SPICE, 2022, a toolkit developed by NASA’s Navigation and Ancillary Information Facility at the Jet Propulsion Laboratory. Other similar, smaller projects that appeared later on and that are still maintained to this day include PyKEP Izzo et al., 2020, beyond beyond, 2022, tudatpy tudatpy, 2022, sbpy Mommert et al., 2019, Skyfield Rhodes, 2020 (Python), CelestLab (Scilab) celestlab, 2022, astrodynamics.jl (Julia) Astrodynamics.jl, n.d. and Nyx (Rust) nyx, 2021. In addition, there are some Graphical User Interface (GUI) based open source programs used for Mission Analysis and orbit visualization, such as GMAT GMAT, 2020 and gpredict gpredict, 2018, and complete web applications for tracking constellations of satellites like the SatNOGS project by the Libre Space Foundation SatNOGS, 2021.

The level of quality and maintenance of these packages is somewhat heterogeneous. Community-led projects with a strong corporate backing like Orekit are in excellent health, while on the other hand smaller projects developed by volunteers (beyond, astrodynamics.jl) or with limited institutional support (PyKEP, GMAT) suffer from lack of maintenance. Part of the problem might stem from the fact that most scientists are never taught how to build software efficiently, let alone the skills to collaboratively develop software in the open Wilson et al., 2014, and astrodynamicists are no exception.

On the other hand, it is often difficult to translate the advances in astrodynamics research to software. Classical algorithms developed throughout the 20th century are described in papers that are sometimes difficult to find, and source code or validation data is almost never available. When it comes to modern research carried in the digital era, source code and validation data is still difficult, even though they are supposedly provided “upon reasonable request” Stodden et al., 2018Gabelica et al., 2022.

It is no surprise that astrodynamics software often requires deep expertise. However, there are often implicit assumptions that are not documented with an adequate level of detail which originate widespread misconceptions and lead even seasoned professionals to make conceptual mistakes. Some of the most notorious misconceptions arise around the use of general perturbations data (OMMs and TLEs) Finkleman, 2007, the geometric interpretation of the mean anomaly Battin, 1999, or coordinate transformations Vallado et al., 2006.

Finally, few of the open source software libraries mentioned above are amenable to scripting or interactive use, as promoted by computational notebooks like Jupyter Kluyver et al., 2016.

The following sections will now discuss the various areas of current research that an astrodynamicist will engage in, and how poliastro improves their workflow.

Methods¶

Software Architecture¶

The architecture of poliastro emerges from the following set of conflicting requirements:

1. There should be a high-level API that enables users to perform orbital calculations in a straightforward way and prevent typical mistakes.
2. The running time of the algorithms should be within the same order of magnitude of existing compiled implementations.
3. The library should be written in a popular open-source language to maximize adoption and lower the barrier to external contributors.

One of the most typical mistakes we set ourselves to prevent with the high-level API is dimensional errors. Addition and substraction operations of physical quantities are defined only for quantities with the same units Drobot, 1953: for example, the operation $1~\text{km} + 100~\text{m}$ requires a scale transformation of at least one of the operands, since they have different units (kilometers and meters) but the same dimension (length), whereas the operation $1~\text{km} + 1~\text{kg}$ is directly not allowed because dimensions are incompatible (length and mass). As such, software systems operating with physical quantities should raise exceptions when adding different dimensions, and transparently perform the required scale transformations when adding different units of the same dimension.

With this in mind, we evaluated several Python packages for unit handling (see J. Goldbaum et al. (2018) for a recent survey) and chose astropy.units The Astropy Collaboration et al., 2018.

radius = 6000  # km
altitude = 500  # m

# Wrong!
distance = radius + altitude

from astropy import units as u

# Correct
distance = (radius << u.km) + (altitude << u.m)

This notion of providing a “safe” API extends to other parts of the library by leveraging other capabilities of the Astropy project. For example, timestamps use astropy.time objects, which take care of the appropriate handling of time scales (such as TDB or UTC), reference frame conversions leverage astropy.coordinates, and so forth.

One of the drawbacks of existing unit packages is that they impose a significant performance penalty. Even though astropy.units is integrated with NumPy, hence allowing the creation of array quantities, all the unit compatibility checks are implemented in Python and require lots of introspection, and this can slow down mathematical operations by several orders of magnitude. As such, to fulfill our desired performance requirement for poliastro, we envisioned a two-layer architecture:

• The Core API follows a procedural style, and all the functions receive Python numerical types and NumPy arrays for maximum performance.
• The High level API is object-oriented, all the methods receive Astropy Quantity objects with physical units, and computations are deferred to the Core API.

Most of the methods of the High level API consist only of the necessary unit compatibility checks, plus a wrapper over the corresponding Core API function that performs the actual computation.

@u.quantity_input(E=u.rad, ecc=u.one)
def E_to_nu(E, ecc):
    """True anomaly from eccentric anomaly."""
    return (
        E_to_nu_fast(
            E.to_value(u.rad),
            ecc.value
        ) << u.rad
    ).to(E.unit)

As a result, poliastro offers a unit-safe API that performs the least amount of computation possible to minimize the performance penalty of unit checks, and also a unit-unsafe API that offers maximum performance at the cost of not performing any unit validation checks.

Finally, there are several options to write performant code that can be used from Python, and one of them is using a fast, compiled language for the CPU intensive parts. Successful examples of this include NumPy, written in C Harris et al., 2020, SciPy, featuring a mix of FORTRAN, C, and C++ code Virtanen et al., 2020, and pandas, making heavy use of Cython Behnel et al., 2011. However, having to write code in two different languages hinders the development speed, makes debugging more difficult, and narrows the potential contributor base (what Julia creators called “The Two Language Problem” Bezanson et al., 2017).

As authors of poliastro we wanted to use Python as the sole programming language of the implementation, and the best solution we found to improve its performance was to use Numba, a LLVM-based Python JIT compiler Lam et al., 2015.

Usage¶

Basic Orbit and Ephem creation¶

The two central objects of the poliastro high level API are Orbit and Ephem:

• Orbit objects represent an osculating (hence Keplerian) orbit of a dimensionless object around an attractor at a given point in time and a certain reference frame.
• Ephem objects represent an ephemerides, a sequence of spatial coordinates over a period of time in a certain reference frame.

There are six parameters that uniquely determine a Keplerian orbit, plus the gravitational parameter of the corresponding attractor ($k$ or μ). Optionally, an epoch that contextualizes the orbit can be included as well. This set of six parameters is not unique, and several of them have been developed over the years to serve different purposes. The most widely used ones are:

• Cartesian elements: Three components for the position $(x, y, z)$ and three components for the velocity $(v_x, v_y, v_z)$. This set has no singularities.
• Classical Keplerian elements: Two components for the shape of the conic (usually the semimajor axis $a$ or semiparameter $p$ and the eccentricity $e$), three Euler angles for the orientation of the orbital plane in space (inclination $i$, right ascension of the ascending node Ω, and argument of periapsis ω), and one polar angle for the position of the body along the conic (usually true anomaly $f$ or ν). This set of elements has an easy geometrical interpretation and the advantage that, in pure two-body motion, five of them are fixed $(a, e, i, \Omega, \omega)$ and only one is time-dependent (ν), which greatly simplifies the analytical treatment of orbital perturbations. However, they suffer from singularities steming from the Euler angles (“gimbal lock”) and equations expressed in them are ill-conditioned near such singularities.
• Walker modified equinoctial elements: Six parameters $(p, f, g, h, k, L)$. Only $L$ is time-dependent and this set has no singularities, however the geometrical interpretation of the rest of the elements is lost Walker et al., 1985.

Here is how to create an Orbit from cartesian and from classical Keplerian elements. Walker modified equinoctial elements are supported as well.

from astropy import units as u

from poliastro.bodies import Earth, Sun
from poliastro.twobody import Orbit
from poliastro.constants import J2000

# Data from Curtis, example 4.3
r = [-6045, -3490, 2500] << u.km
v = [-3.457, 6.618, 2.533] << u.km / u.s

orb_curtis = Orbit.from_vectors(
   Earth,  # Attractor
   r, v  # Elements
)

# Data for Mars at J2000 from JPL HORIZONS
a = 1.523679 << u.au
ecc = 0.093315 << u.one
inc = 1.85 << u.deg
raan = 49.562 << u.deg
argp = 286.537 << u.deg
nu = 23.33 << u.deg

orb_mars = Orbit.from_classical(
   Sun,
   a, ecc, inc, raan, argp, nu,
   J2000  # Epoch
)

When displayed on an interactive REPL, Orbit objects provide basic information about the geometry, the attractor, and the epoch:

>>> orb_curtis
7283 x 10293 km x 153.2 deg (GCRS) orbit
around Earth (X) at epoch J2000.000 (TT)

>>> orb_mars
1 x 2 AU x 1.9 deg (HCRS) orbit
around Sun (X) at epoch J2000.000 (TT)

Similarly, Ephem objects can be created using a variety of classmethods as well. Thanks to astropy.coordinates built-in low-fidelity ephemerides, as well as its capability to remotely access the JPL HORIZONS system, the user can seamlessly build an object that contains the time history of the position of any Solar System body:

from astropy.time import Time
from astropy.coordinates import solar_system_ephemeris

from poliastro.ephem import Ephem

# Configure high fidelity ephemerides globally
# (requires network access)
solar_system_ephemeris.set("jpl")

# For predefined poliastro attractors
earth = Ephem.from_body(Earth, Time.now().tdb)

# For the rest of the Solar System bodies
ceres = Ephem.from_horizons("Ceres", Time.now().tdb)

There are some crucial differences between Orbit and Ephem objects:

• Orbit objects have an attractor, whereas Ephem objects do not. Ephemerides can originate from complex trajectories that don’t necessarily conform to the ideal two-body problem.
• Orbit objects capture a precise instant in a two-body motion plus the necessary information to propagate it forward in time indefinitely, whereas Ephem objects represent a bounded time history of a trajectory. This is because the equations for the two-body motion are known, whereas an ephemeris is either an observation or a prediction that cannot be extrapolated in any case without external knowledge. As such, Orbit objects have a .propagate method, but Ephem ones do not. This prevents users from attempting to propagate the position of the planets, which will always yield poor results compared to the excellent ephemerides calculated by external entities.

Finally, both types have methods to convert between them:

• Ephem.from_orbit is the equivalent of sampling a two-body motion over a given time interval. As explained above, the resulting Ephem loses the information about the original attractor.
• Orbit.from_ephem is the equivalent of calculating the osculating orbit at a certain point of a trajectory, assuming a given attractor. The resulting Orbit loses the information about the original, potentially complex trajectory.

Orbit propagation¶

Orbit objects have a .propagate method that takes an elapsed time and returns another Orbit with new orbital elements and an updated epoch:

>>> from poliastro.examples import iss

>>> iss
>>> 6772 x 6790 km x 51.6 deg (GCRS) ...

>>> iss.nu.to(u.deg)
<Quantity 46.59580468 deg>

>>> iss_30m = iss.propagate(30 << u.min)

>>> (iss_30m.epoch - iss.epoch).datetime
datetime.timedelta(seconds=1800)

>>> (iss_30m.nu - iss.nu).to(u.deg)
<Quantity 116.54513153 deg>

The default propagation algorithm is an analytical procedure described in Farnocchia et al. (2013) that works seamlessly in the near parabolic region. In addition, poliastro implements analytical propagation algorithms as described in Danby & Burkardt (1983)Odell & Gooding (1986)Markley (1995)Mikkola (1987)Pimienta-Penalver (2013)Charls (2022)Vallado & McClain (2007).

Natural perturbations¶

As showcased in Figure 2, at any point in a trajectory we can define an ideal Keplerian orbit with the same position and velocity under the attraction of a point mass: this is called the osculating orbit. Some numerical propagation methods exist that model the true, perturbed orbit as a deviation from an evolving, osculating orbit. poliastro implements Cowell’s method Cowell & Crommelin, 1910, which consists in adding all the perturbation accelerations and then integrating the resulting differential equation with any numerical method of choice:

The resulting equation is usually integrated using high order numerical methods, since the integration times are quite large and the tolerances comparatively tight. An in-depth discussion of such methods can be found in Hairer et al. (2009). poliastro uses Dormand-Prince 8(5,3) (DOP853), a commonly used method available in SciPy Harris et al., 2020.

There are several natural perturbations included: J2 and J3 gravitational terms, several atmospheric drag models (exponential, Jacchia (1977), Atmosphere et al. (1962), Atmosphere et al. (1976)), and helpers for third body gravitational attraction and radiation pressure as described in Curtis (2010).

@njit
def combined_a_d(
    t0, state, k, j2, r_eq, c_d, a_over_m, h0, rho0
):
    return (
        J2_perturbation(
            t0, state, k, j2, r_eq
        ) + atmospheric_drag_exponential(
            t0, state, k, r_eq, c_d, a_over_m, h0, rho0
        )
    )

def f(t0, state, k):
    du_kep = func_twobody(t0, state, k)
    ax, ay, az = combined_a_d(
        t0,
        state,
        k,
        R=R,
        C_D=C_D,
        A_over_m=A_over_m,
        H0=H0,
        rho0=rho0,
        J2=Earth.J2.value,
    )
    du_ad = np.array([0, 0, 0, ax, ay, az])

    return du_kep + du_ad

rr = propagate(
    orbit,
    tofs,
    method=cowell,
    f=f,
)

Continuous thrust control laws¶

Beyond natural perturbations, spacecraft can modify their trajectory on purpose by using impulsive maneuvers (as explained in the next section) as well as continuous thrust guidance laws. The user can define custom guidance laws by providing a perturbation acceleration in the same way natural perturbations are used. In addition, poliastro includes several analytical solutions for continuous thrust guidance laws with specific purposes, as studied in Cano Rodríguez (2017): optimal transfer between circular coplanar orbits Edelbaum, 1961Burt, 1967, optimal transfer between circular inclined orbits Edelbaum, 1961Kechichian, 1997, quasi-optimal eccentricity-only change Pollard, 1997, simultaneous eccentricity and inclination change Pollard, 2000, and agument of periapsis adjustment Pollard, 1998. A much more rigorous analysis of a similar set of laws can be found in Di Carlo & Vasile (2021).

from poliastro.twobody.thrust import change_ecc_inc

ecc_f = 0.0 << u.one
inc_f = 20.0 << u.deg
f = 2.4e-6 << (u.km / u.s**2)

a_d, _, t_f = change_ecc_inc(orbit, ecc_f, inc_f, f)

Impulsive maneuvers¶

Impulsive maneuvers are modeled considering a change in the velocity of a spacecraft while its position remains fixed. The poliastro.maneuver.Maneuver class provides various constructors to instantiate popular impulsive maneuvers in the framework of the non-perturbed two-body problem:

• Maneuver.impulse
• Maneuver.hohmann
• Maneuver.bielliptic
• Maneuver.lambert

from poliastro.maneuver import Maneuver

orb_i = Orbit.circular(Earth, alt=700 << u.km)
hoh = Maneuver.hohmann(orb_i, r_f=36000 << u.km)

Once instantiated, Maneuver objects provide information regarding total $\Delta v$ and $\Delta t$:

>>> hoh.get_total_cost()
<Quantity 3.6173981270031357 km / s>

>>> hoh.get_total_time()
<Quantity 15729.741535747102 s>

Maneuver objects can be applied to Orbit instances using the apply_maneuver method.

>>> orb_i
7078 x 7078 km x 0.0 deg (GCRS) orbit
around Earth (X)

>>> orb_f = orb_i.apply_maneuver(hoh)
>>> orb_f
36000 x 36000 km x 0.0 deg (GCRS) orbit
around Earth (X)

Targeting¶

Targeting is the problem of finding the orbit connecting two positions over a finite amount of time. Within the context of the non-perturbed two-body problem, targeting is just a matter of solving the BVP, also known as Lambert’s problem. Because targeting tries to find for an orbit, the problem is included in the Initial Orbit Determination field.

The poliastro.iod package contains izzo and vallado modules. These provide a lambert function for solving the targeting problem. Nevertheless, a Maneuver.lambert constructor is also provided so users can keep taking advantage of Orbit objects.

# Declare departure and arrival datetimes
date_launch = time.Time(
    '2011-11-26 15:02', scale='tdb'
)
date_arrival = time.Time(
    '2012-08-06 05:17', scale='tdb'
)

# Define initial and final orbits
orb_earth = Orbit.from_ephem(
    Sun, Ephem.from_body(Earth, date_launch),
    date_launch
)
orb_mars = Orbit.from_ephem(
    Sun, Ephem.from_body(Mars, date_arrival),
    date_arrival
)

# Compute targetting maneuver and apply it
man_lambert = Maneuver.lambert(orb_earth, orb_mars)
orb_trans, orb_target = ss0.apply_maneuver(
    man_lambert, intermediate=true
)

Targeting is closely related to quick mission design by means of porkchop diagrams. These are contour plots showing all combinations of departure and arrival dates with the specific energy for each transfer orbit. They allow for quick identification of the most optimal transfer dates between two bodies.

The poliastro.plotting.porkchop provides the PorkchopPlotter class which allows the user to generate these diagrams.

from poliastro.plotting.porkchop import (
    PorkchopPlotter
)
from poliastro.utils import time_range

# Generate all launch and arrival dates
launch_span = time_range(
    "2020-03-01", end="2020-10-01", periods=int(150)
)
arrival_span = time_range(
    "2020-10-01", end="2021-05-01", periods=int(150)
)

# Create an instance of the porkchop and plot it
porkchop = PorkchopPlotter(
    Earth, Mars, launch_span, arrival_span,
)

Previous code, with some additional customization, generates Figure 3.

Plotting¶

For visualization purposes, poliastro provides the poliastro.plotting package, which contains various utilities for generating 2D and 3D graphics using different backends such as matplotlib Hunter, 2007 and Plotly Plotly Technologies Inc., 2015.

Generated graphics can be static or interactive. The main difference between these two is the ability to modify the camera view in a dynamic way when using interactive plotters.

The most important classes in the poliastro.plotting package are StaticOrbitPlotter and OrbitPlotter3D. In addition, the poliastro.plotting.misc module contains the plot_solar_system function, which allows the user to visualize inner and outter both in 2D and 3D, as requested by users.

The following example illustrates the plotting capabilities of poliastro. At first, orbits to be plotted are computed and their plotting style is declared:

from poliastro.plotting.misc import plot_solar_system

# Current datetime
now = Time.now().tdb

# Obtain Florence and Halley orbits
florence = Orbit.from_sbdb("Florence")
halley_1835_ephem = Ephem.from_horizons(
    "90000031", now
)
halley_1835 = Orbit.from_ephem(
    Sun, halley_1835_ephem, halley_1835_ephem.epochs[0]
)

# Define orbit labels and color style
florence_style = {label: "Florence", color: "#000000"}
halley_style = {label: "Florence", color: "#84B0B8"}

The static two-dimensional plot can be created using the following code:

# Generate a static 2D figure
frame2D = rame = plot_solar_system(
    epoch=now, outer=False
)
frame2D.plot(florence, **florence_style)
frame2D.plot(florence, **halley_style)

As a result, Figure 4 is obtained.

The interactive three-dimensional plot can be created using the following code:

# Generate an interactive 3D figure
frame3D = rame = plot_solar_system(
    epoch=now, outer=False,
    use_3d=True, interactive=true
)
frame3D.plot(florence, **florence_style)
frame3D.plot(florence, **halley_style)

As a result, Figure 5 is obtained.

Commercial Earth satellites¶

Figure 6 gives a clear picture of the most important natural perturbations affecting satellites in LEO, namely: the first harmonic of the geopotential field $J_2$ (representing the attractor oblateness), the atmospheric drag, and the higher order harmonics of the geopotential field.

At least the most significant of these perturbations need to be taken into account when propagating LEO orbits, and therefore the methods for purely Keplerian motion are not enough. As seen above, poliastro implements a number of these perturbations already - however, numerical methods are much slower than analytical ones, and this can render them unsuitable for large scale simulations, satellite conjunction assesment, propagation in constrained hardware, and so forth.

To address this issue, semianalytical propagation methods were devised that attempt to strike a balance between the fast running times of analytical methods and the necessary inclusion of perturbation forces. One of such semianalytical methods are the Simplified General Perturbation (SGP) models, first developed in Hilton & Kuhlman (1966) and then refined in Lane & Cranford (1969) into what we know these days as the SGP4 propagator Hoots & Roehrich, 1980Vallado et al., 2006. Even though certain elements of the reference frame used by SGP4 are not properly specified Vallado et al., 2006 and that its accuracy might still be too limited for certain applications Kelso & others, 2009Lara, 2016, it is nowadays the most widely used propagation method thanks in large part to the dissemination of General Perturbations orbital data by the US 501(c)(3) CelesTrak (which itself obtains it from the 18th Space Defense Squadron of the US Space Force).

The starting point of SGP4 is a special element set that uses Brouwer mean orbital elements Brouwer, 1959 plus a ballistic coefficient based on an approximation of the atmospheric drag Lane & Cranford, 1969, and its results are expressed in a special coordinate system called True Equator Mean Equinox (TEME). Special care needs to be taken to avoid mixing mean elements with osculating elements, and to convert the output of the propagation to the appropriate reference frame. These element sets have been traditionally distributed in a compact text representation called Two-Line Element sets (TLEs) (see Figure 7 for an example). However this format is quite cryptic and suffers from a number of shortcomings, so recently there has been a push to use the Orbit Data Messages international standard developed by the Consultive Committee for Space Data Systems (CCSDS 502.0-B-2).

At the moment, general perturbations data both in OMM and TLE format can be integrated with poliastro thanks to the sgp4 Python library and the Ephem class as follows:

from astropy.coordinates import TEME, GCRS

from poliastro.ephem import Ephem
from poliastro.frames import Planes


def ephem_from_gp(sat, times):
    errors, rs, vs = sat.sgp4_array(times.jd1, times.jd2)
    if not (errors == 0).all():
        warn(
            "Some objects could not be propagated, "
            "proceeding with the rest",
            stacklevel=2,
        )
        rs = rs[errors == 0]
        vs = vs[errors == 0]
        times = times[errors == 0]

    cart_teme = CartesianRepresentation(
        rs << u.km,
        xyz_axis=-1,
        differentials=CartesianDifferential(
            vs << (u.km / u.s),
            xyz_axis=-1,
        ),
    )
    cart_gcrs = (
        TEME(cart_teme, obstime=times)
        .transform_to(GCRS(obstime=times))
        .cartesian
    )

    return Ephem(
        cart_gcrs,
        times,
        plane=Planes.EARTH_EQUATOR
    )

However, no native integration with SGP4 has been implemented yet in poliastro, for technical and non-technical reasons. On one hand, this propagator is too different from the other methods, and we have not yet devised how to add it to the library in a way that does not create confusion. On the other hand, adding such a propagator to poliastro would probably open the flood gates of corporate users of the library, and we would like to first devise a sustainability strategy for the project, which is addressed in the next section.

Future work¶

Despite the fact that poliastro has existed for almost a decade, for most of its history it has been developed by volunteers on their free time, and only in the past five years it has received funding through various Summer of Code programs (SOCIS 2017, GSOC 2018-2021) and institutional grants (NumFOCUS 2020, 2021). The funded work has had an overwhemingly positive impact on the project, however the lack of a dedicated maintainer has caused some technical debt to accrue over the years, and some parts of the project are in need of refactoring or better documentation.

Historically, poliastro has tried to implement algorithms that were applicable for all the planets in the Solar System, however some of them have proved to be very difficult to generalize for bodies other than the Earth. For cases like these, poliastro ships a poliastro.earth package, but going forward we would like to continue embracing a generic approach that can serve other bodies as well.

Several open source projects have successfully used poliastro or were created taking inspiration from it, like spacetech-ssa by IBM^[1] or mubody Bermejo Ballesteros et al., 2022. AGI (previously Analytical Graphics, Inc., now Ansys Government Initiatives) published a series of scripts to automate the commercial tool STK from Python leveraging poliastro ^[2]. However, we have observed that there is still lots of repeated code across similar open source libraries written in Python, which means that there is an opportunity to provide a “kernel” of algorithms that can be easily reused. Although poliastro.core started as a separate layer to isolate fast, non-safe functions as described above, we think we could move it to an external package so it can be depended upon by projects that do not want to use some of the higher level poliastro abstractions or drag its large number of heavy dependencies.

Finally, the sustainability of the project cannot yet be taken for granted: the project has reached a level of complexity that already warrants dedicated development effort that cannot be covered with short-lived grants. Such funding could potentially come from the private sector, but although there is evidence that several for-profit companies are using poliastro, we have very little information of how is it being used and what problems are those users having, let alone what avenues for funded work could potentially work. Organizations like the Libre Space Foundation advocate for a strong copyleft licensing model to convince commercial actors to contribute to the commons, but in principle that goes against the permissive licensing that the wider Scientific Python ecosystem, including poliastro, has adopted. With the advent of new business models and the ever increasing reliance in open source by the private sector, a variety of ways to engage commercial users and include them in the conversation exist. However, these have not been explored yet.

References¶