poliastro: a Python library for interactive astrodynamics

—Space is more popular than ever, with the growing public awareness of interplanetary scientiﬁc missions, as well as the increasingly large number of satellite companies planning to deploy satellite constellations. Python has become a fundamental technology in the astronomical sciences, and it has also caught the attention of the Space Engineering community. One of the requirements for designing a space mission is studying the trajectories of satellites, probes, and other artiﬁcial objects, usually ignoring non-gravitational forces or treating them as perturbations: the so-called n-body problem. However, for preliminary design studies and most practical purposes, it is sufﬁcient to consider only two bodies: the object under study and its attractor. Even though the two-body problem has many analytical solutions, orbit propagation (the initial value problem) and targeting (the boundary value problem) remain computationally intensive because of long propagation times, tight tolerances, and vast solution spaces. On the other hand, astrodynamics researchers often do not share the source code they used to run analyses and simulations, which makes it challenging to try out new solutions. This paper presents poliastro, an open-source Python library for interactive astrodynamics that features an easy-to-use API and tools for quick visualization. poliastro implements core astrodynamics algorithms (such as the resolution of the Kepler and Lambert problems) and leverages numba, a Just-in-Time compiler for scientiﬁc Python, to optimize the running time. Thanks to Astropy, poliastro can perform seamless coordinate frame conversions and use proper physical units and timescales. At the moment, poliastro is the longest-lived Python library for astrodynamics, has contributors from all around the world, and several New Space companies and people in academia use it.


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 ( [Dub73], [Her71]).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 closedform solution exists for the n-body problem, since the resulting dynamical system is chaotic [Bat99].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.∆v = v e ln m 0 m f (2) 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 ∈ 1, ..., n bodies subject to their mutual attraction, by application of Newton's law of universal gravitation, the total force f i affecting m i due to the presence of the other n − 1 masses is given by [Bat99]: where G = 6.67430 • 10 −11 N m 2 kg −2 is the universal gravitational constant, and r i j 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 µ = G(m 1 + m 2 ) = G(M + m).When m M (for example, an artificial satellite orbiting a planet), one can consider µ = GM a property of the attractor.

Keplerian vs non-keplerian motion
Conveniently manipulating equation 5 leads to several properties [Bat99] 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 = p 1+e cos ν (Kepler's first law).
2) The magnitude of the specific angular momentum h = r 2 d θ dt 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π a 3 µ (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 ≤ e < 1) it is the Kepler equation, and equivalent relations exist for the other eccentricity regimes [Bat99].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") [Bat99].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 [Con14].

•
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 [noa22c] 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 [WAB + 14], 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" [SSM18] [GBP22].
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) [Fin07], the geometric interpretation of the mean anomaly [Bat99], or coordinate transformations [VCHK06].
Finally, few of the open source software libraries mentioned above are amenable to scripting or interactive use, as promoted by computational notebooks like Jupyter [KRKP + 16].
The following sections will now discuss the various areas of current research that an astrodynamicist will engage in, and how poliastro improves their workflow.

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 [Dro53]: for example, the operation 1 km + 100 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 km + 1 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 [JGAZJT + 18] for a recent survey) and chose astropy.units[TPWS + 18].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.timeobjects, 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.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 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 [LPS15].

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, Ω, ω) 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 [WIO85].
Here is how to create an Orbit from cartesian and from classical Keplerian elements.Walker modified equinoctial elements are supported as well.Similarly, Ephem objects can be created using a variety of classmethods as well.Thanks to astropy.coordinatesbuilt-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.ephemimport 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 .propagatemethod, 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 .propagatemethod that takes an elapsed time and returns another Orbit with new orbital elements and an updated epoch: The default propagation algorithm is an analytical procedure described in [FCM13] that works seamlessly in the near parabolic region.In addition, poliastro implements analytical propagation algorithms as described in [  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 [CC10], 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 [HNW09].poliastro uses Dormand-Prince 8(5,3) (DOP853), a commonly used method available in SciPy [HMvdW + 20].
There are several natural perturbations included: J2 and J3 gravitational terms, several atmospheric drag models (exponential, [Jac77]

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 [CR17]: optimal transfer between circular coplanar orbits [Ede61] [Bur67], optimal transfer between circular inclined orbits [Ede61] [Kec97], quasi-optimal eccentricity-only change [Pol97], simultaneous eccentricity and inclination change [Pol00], and agument of periapsis adjustment [Pol98].A much more rigorous analysis of a similar set of laws can be found in [DCV21].Maneuver objects can be applied to Orbit instances using the apply_maneuver method.

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.iodpackage contains izzo and vallado modules.These provide a lambert function for solving the targeting problem.Nevertheless, a Maneuver.lambertconstructor is also provided so users can keep taking advantage of Orbit objects.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.porkchopprovides the PorkchopPlotter class which allows the user to generate these diagrams.Previous code, with some additional customization, generates figure 3.

Plotting
For visualization purposes, poliastro provides the poliastro.plottingpackage, which contains various utilities for generating 2D and 3D graphics using different backends such as matplotlib [Hun07] and Plotly [Inc15].

Launch date
The most important classes in the poliastro.plottingpackage are StaticOrbitPlotter and OrbitPlotter3D.In addition, the poliastro.plotting.miscmodule 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: The interactive three-dimensional plot can be created using the following code:

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 starting point of SGP4 is a special element set that uses Brouwer mean orbital elements [Bro59] plus a ballistic coefficient based on an approximation of the atmospheric drag [LC69], 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 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).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.earthpackage, 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 [BBVPFSC22].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.corestarted as a separate layer to isolate fast, nonsafe 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.

@
parts.Successful examples of this include NumPy, written in C [HMvdW + 20], SciPy, featuring a mix of FORTRAN, C, and C++ code [VGO + 20], and pandas, making heavy use of Cython [BBC + 11].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" [BEKS17]).

Fig. 3 :
Fig. 3: Porkchop plot for Earth-Mars transfer arrival energy showing latest missions to the Martian planet.

Fig. 4 :
Fig. 4: Two-dimensional view of the inner Solar System, Florence, and Halley.

Fig. 5 :
Fig. 5: Three-dimensional view of the inner Solar System, Florence, and Halley.
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: