∎ scaletikzpicturetowidth[1]

11institutetext: Thomas Caleb (corresponding author) 22institutetext: ISAE-
SUPAERO, 10 Avenue Édouard Belin, 54032 - 31055, Toulouse, France
OrcID: 0000-0002-4027-3340
22email<EMAIL_ADDRESS>33institutetext: Gianmario
Merisio 44institutetext: Department of Aerospace Science and Technology,
Politecnico di Milano, Via La Masa 34, 20156, Milano, Italy
OrcID: 0000-0001-8806-7952
44email<EMAIL_ADDRESS>55institutetext: Pierluigi di Lizia
66institutetext: Department of Aerospace Science and Technology, Politecnico
di Milano, Via La Masa 34, 20156, Milano, Italy
OrcID: 0000-0003-1692-3929
66email<EMAIL_ADDRESS>77institutetext: Francesco Topputo
88institutetext: Department of Aerospace Science and Technology, Politecnico
di Milano, Via La Masa 34, 20156, Milano, Italy
OrcID: 0000-0002-5369-6887
88email<EMAIL_ADDRESS>
# Stable sets mapping with Taylor differential algebra with application to
ballistic capture orbits around Mars

T. Caleb  G. Merisio  P. di Lizia  F. Topputo

(Received: date / Accepted: date)

###### Abstract

Ballistic capture orbits offer safer Mars injection at longer transfer time.
However, the search for such an extremely rare event is a computationally-
intensive process. Indeed, it requires the propagation of a grid sampling the
whole search space. This work proposes a novel ballistic capture search
algorithm based on Taylor differential algebra propagation. This algorithm
provides a continuous description of the search space compared to classical
grid sampling research and focuses on areas where the nonlinearities are the
largest.Macroscopic analyses have been carried out to obtain cartography of
large sets of solutions. Two criteria, named consistency and quality, are
defined to assess this new algorithm and to compare its performances with
classical grid sampling of the search space around Mars. Results show that
differential algebra mapping works on large search spaces, and automatic
domain splitting captures the dynamical variations on the whole domain
successfully. The consistency criterion shows that more than 87% of the search
space is guaranteed as accurate, with the quality criterion kept over 80%.

###### Keywords:

Ballistic capture Mars Stable sets Taylor differential algebra Weak stability
boundary

††journal: Celestial Mechanics and Dynamical Astronomy

##  1 Introduction

Ballistic capture (BC) allows a spacecraft to approach a planet and enter a
temporary orbit about it without requiring maneuvers in between. As part of
the low-energy transfers, it is a valuable alternative to Keplerian
approaches. Exploiting BC grants several benefits in terms of both cost
reduction (Belbruno and Miller, 1993) and mission versatility (Belbruno and
Carrico, 2000; Topputo and Belbruno, 2015), in general at the cost of longer
transfer times (Circi and Teofilatto, 2001; Ivashkin, 2002). In the past, the
BC mechanism was used to rescue Hiten (Belbruno and Miller, 1990), and to
design insertion trajectories in lunar missions like SMART-1 (Racca et al,
2002) and GRAIL (Chung et al, 2010). In the near future, BepiColombo will
exploit BC orbits to be weakly captured by Mercury (Benkhoff et al, 2021;
Schuster and Jehn, 2014). BC is an event occurring in extremely rare occasions
and requires acquiring a proper state (position and velocity) far away from
the target planet (Topputo and Belbruno, 2015). In fact, massive numerical
simulations are required to find the specific conditions that support capture
(Topputo and Belbruno, 2009) and only approximately 1 out of $10\,000$ states
lead to capture (Luo and Topputo, 2015). In a first effort to reduce the
computational burden, the variational theory for Lagrangian coherent
structures (Haller, 2011) was recently applied to find BC opportunities more
efficiently (Manzi and Topputo, 2021).

Differential algebra (DA) propagation is a worthy candidate to reduce the
computational burden for the search of BC trajectories. It consists of
propagating initial conditions (ICs), not as a single point but as an interval
around an IC. Thanks to the Taylor expansions of the flow, the state of any
point in the represented interval can be determined through convenient
polynomial evaluations (Berz, 1992, 1999). This gain in efficiency comes at
the cost of a loss of accuracy due to the finite Taylor expansion. The
automatic domain splitting (ADS) algorithm allows to represent large domains
accurately when increasing the order of the polynomial expansions fails to do
so (Wittig et al, 2015). Indeed, on poorly-defined domains, raising the order
will increase the approximation error on the edges while the error will
decrease on already well-defined areas. ADS splits the initial domain into
smaller sub-domains to reduce the approximation error when it grows above a
given tolerance.

DA propagation is increasingly used in astrodynamics. Indeed, it provides high
performances in uncertainty propagation in the two-body dynamics (Valli et al,
2013), even for highly nonlinear dynamics with large uncertainties when
exploiting ADS, as in the case of Apophis (Wittig et al, 2015) or Apollo LM-10
also known as Snoopy (Caleb and Lizy-Destrez, 2021). Other applications occur
in orbital mechanics, such as propagation of probability density functions
(Wittig et al, 2017), maximum a posteriori estimation (Servadio et al, 2022),
orbit determination (Pirovano et al, 2021; Servadio and Zanetti, 2021), and
generation and study of orbit families in the circular restricted three body
problem (CR3BP) (Di Lizia et al, 2008; Baresi et al, 2021).

DA propagation allows avoiding intensive grid sampling implied by point-wise
research of conventional algorithms for designing BC trajectories (Hyeraci and
Topputo, 2010; Luo et al, 2014), as it offers a continuous description of the
whole search-space. In addition, the computation of a Taylor expansion
provides information such as the partial derivatives of the flow up to an
arbitrary order (Wittig et al, 2017). Furthermore, the polynomial maps can be
manipulated to impose constraints on the flown trajectories (Berz, 1999; Di
Lizia et al, 2008).

The goal of this work is to use DA mapping to carry out macroscopic analyses
of the phase space about Mars to find BC trajectories. Hence, the adaptation
of the definition of BC for DA, followed by the definition of two criteria,
named consistency and quality, to assess the performances of DA-based mapping
of BC compared to point-wise mapping. Cartography of large BC sets about Mars
are computed using DA mapping. Mars is chosen without loss of generality due
to its relevance in the long-term exploration. The work proposes an
alternative classification algorithm to sort the sub-domains produced by the
ADS algorithm in newly defined capture sets. While DA mapping allows
performing good macroscopic cartography of the search space, it does not allow
to fully replace point-wise mapping. However, DA mapping has the advantage of
being continuous, so meaning that the behavior of any point in the search
space is defined. On the contrary, a point-wise mapping delivers precise
knowledge on the discrete set of propagated ICs only. As a consequence, the
behavior of the continuous space between two points is unknown and can be
challenging to interpolate due to the nonlinear dynamics.

The remainder of the paper is organized as follows. In Section 2, the
dynamical model employed is introduced, as well as the weak stability boundary
(WSB) concept, the BC mechanism, and the DA propagation. Then, the description
of the characterization process, and the mapping-assessment methodology follow
in Section 3. Results are presented and discussed in Section 4. Eventually,
conclusions are drawn in Section 5 together with the presentation of future
work.

##  2 Background

Details about dynamical model, WSB, BC phenomenon, DA propagation, and ADS
algorithms are herewith presented.

###  2.1 Dynamical model

According to the nomenclature introduced in Luo et al (2014), a target (also
referred to as central body) and a primary are defined. The target being the
body around which the motion of the spacecraft is studied (Mars in this work),
and the primary being the body around which the target revolves (the Sun).
Target and primary masses are $m_{t}$ and $m_{p}$, respectively.

####  2.1.1 Reference frames

In this work, the following reference frames are used: J2000, and RTN@$t_{i}$.

##### J2000.

Defined on the Earth’s mean equator and equinox, the J2000 is an inertial
frame determined from observations of planetary motions which was realized to
coincide almost exactly with the International Celestial Reference Frame
(ICRF) (Archinal et al, 2011). Equations of motion (EoM) are integrated in
this reference frame.

##### RTN@$t_{i}$.

The radial-tangential-normal of date frame (RTN@$t_{i}$) is an inertial frame
frozen at a prescribed epoch $t_{i}$. The frame is centered at the target. The
$x$-axis is aligned with the primary–secondary direction, the $z$-axis is
normal to the primary–secondary plane in the direction of their angular
momentum, and the $y$-axis completes the dextral orthonormal triad. ICs are
defined in this frame (Luo and Topputo, 2015).

####  2.1.2 Ephemerides

Precise states of the Sun and the major planets are retrieved from the Jet
Propulsion Laboratory (JPL)’s planetary ephemerides de440s.bsp111Data publicly
available at:
https://naif.jpl.nasa.gov/pub/naif/generic_kernels/spk/planets/de440s.bsp
[retrieved Feb 1, 2022]. (or DE440s) (Park et al, 2021). Additionally, the
ephemerides mars097.bsp of Mars (the target) and its moons are
employed222~/spk/satellites/mars097.bsp [retrieved Feb 1, 2022].. The
following generic leap seconds kernel (LSK) and planetary constants kernel
(PCK) are used: naif0012.tls, pck00010.tpc, and gm_de431.tpc333Data publicly
available at:
https://naif.jpl.nasa.gov/pub/naif/generic_kernels/lsk/naif0012.tls,
~/generic_kernels/pck/pck00010.tpc, and ~/generic_kernels/pck/gm_de431.tpc
[retrieved Feb 1, 2022]. The gm_de431.tpc PCK kernel is used because the new
version consistent with the ephemerides DE440s has not been released yet..

####  2.1.3 Equations of motion

The EoM used are those of the restricted $N$-body problem. The gravitational
attractions of the Sun, Mercury, Venus, Earth (B444Here B stands for
barycenter.), Mars (central body), Jupiter (B), Saturn (B), Phobos, and Deimos
are considered. Additionally, solar radiation pressure (SRP) is also included
and implemented as a _cannonball_ or _spherical_ model (Scheeres et al, 2011).
The assumed spacecraft specifications needed to evaluate the SRP perturbation
are collected in Table 1. They are compatible with the specifications of a 12U
deep-space CubeSat (Topputo et al, 2021).

Table 1: Assumed spacecraft specifications. Specification | Symbol | Value
---|---|---
Mass | $m$ | $24\text{\,}\mathrm{kg}$
SRP area | $A$ | $0.32\text{\,}{\mathrm{m}}^{2}$
Coefficient of reflectivity | $C_{r}$ | $1.3\text{\,}$

The EoM, written in a non-rotating Mars-centered reference frame are (Luo et
al, 2014; Merisio and Topputo, 2021)

$\ddot{\mathbf{r}}=-\frac{\mu_{t}}{r^{3}}\mathbf{r}-\sum\limits_{i\in\mathbb{P}}\mu_{i}\left(\frac{\mathbf{r}_{i}}{r_{i}^{3}}+\frac{\mathbf{r}-\mathbf{r}_{i}}{\left\lVert\mathbf{r}-\mathbf{r}_{i}\right\rVert^{3}}\right)+\frac{QA}{m}\frac{\mathbf{r}-\mathbf{r}_{\odot}}{\left\lVert\mathbf{r}-\mathbf{r}_{\odot}\right\rVert^{3}}$
(1)

where $\mu_{t}$ is the gravitational parameter of the target body;
$\mathbf{r}$ is the position vector of the spacecraft with respect to the
target and $r$ is its magnitude; $\mathbb{P}$ is a set of $N-2$ indexes each
referring to a perturbing body; $\mu_{i}$ and $\mathbf{r}_{i}$ are the
gravitational parameter and position vector with respect to the target of the
$i$-th body, respectively; $A$ is the Sun-projected area on the spacecraft for
SRP evaluation; $m$ is the spacecraft mass; $\mathbf{r}_{\odot}$ is the
position vector of the Sun with respect to the target. Lastly, $Q$ is equal to

$Q=\frac{LC_{r}}{4\pi c}$ (2)

where $C_{r}$ is the spacecraft coefficient of reflectivity,
$c=$299\,792\,458\text{\,}\mathrm{m}\text{\,}{\mathrm{s}}^{-1}$$ taken from
SPICE (Acton, 1996; Acton et al, 2018) is the speed of light in vacuum, and
$L=S_{\odot}4\pi d_{\mathrm{AU}}^{2}$ is the luminosity of the Sun. The latter
is computed from the solar
constant555https://extapps.ksc.nasa.gov/Reliability/Documents/Preferred_Practices/2301.pdf
[last accessed Feb 1, 2022].
$S_{\odot}=$1367.5\text{\,}\mathrm{W}\text{\,}{\mathrm{m}}^{-2}$$ evaluated at
$d_{\mathrm{AU}}=$149\,597\,870\,613.6889\text{\,}\mathrm{m}$$ corresponding
to $1\text{\,}\mathrm{AU}$ according to SPICE (Acton, 1996; Acton et al,
2018).

####  2.1.4 Numerical integration of EoM

The EoM in Eq. (1) are integrated in their nondimensional form to avoid ill-
conditioning (Luo et al, 2014). Nondimensionalization units are reported in
Table 2. For point-wise simulations, the numerical integration is carried with
the Dormand-Prince 8th-order embedded Runge-Kutta method (DOPRI8) propagation
scheme (Montenbruck and Gill, 2000). It is an adaptive step, 8th-order Runge-
Kutta (RK) integrator with 7th-order error control, the coefficients were
derived by Prince and Dormand (Prince and Dormand, 1981). As for the DA
simulations, the propagation scheme is Dormand-Prince of order 8(5,3) embedded
Runge-Kutta method (DOP853) (Hairer et al, 1993), which is of the same 8th-
order RK integrator family. However, the error control is performed with a
3rd-order and a 5th-order estimation. The dynamics are propagated setting the
relative tolerance to $10^{-12}$ (Luo et al, 2014).

Table 2: Nondimensionalization units. Unit | Symbol | Value | Comment
---|---|---|---
Gravity parameter | $\mathrm{MU}$ | $42\,828.376\text{\,}{\mathrm{km}}^{3}\text{\,}{\mathrm{s}}^{-2}$ | Mars’ gravity parameter $\mu_{t}$
Length | $\mathrm{LU}$ | $3396.0000\text{\,}\mathrm{km}$ | Mars’ radius $R_{\mars}$
Time | $\mathrm{TU}$ | $956.281\,42\text{\,}\mathrm{s}$ | $(\mathrm{LU}^{3}/\mathrm{MU})^{0.5}$
Velocity | $\mathrm{VU}$ | $3.551\,255\,8\text{\,}\mathrm{km}\text{\,}{\mathrm{s}}^{-1}$ | $\mathrm{LU}/\mathrm{TU}$

###  2.2 Weak stability boundary and ballistic capture mechanism

Over the years, the WSB was defined in many different ways. It was initially
identified as a fuzzy boundary region placed at approximately
$1.5\text{\times}{10}^{6}\text{\,}\mathrm{km}$ from the Earth in the Sun–Earth
direction (Belbruno, 1987; Belbruno and Miller, 1990). An algorithmic
definition followed in Belbruno (2004), later extended in García and Gómez
(2007), Topputo and Belbruno (2009), and Sousa Silva and Terra (2012). Then,
the WSB was interpreted as the intersection of three sub-sets of the phase
space (Topputo et al, 2008; Belbruno et al, 2008). The WSB concept being
closely connected to BC (Belbruno, 2004), a formal definition and a
methodology for its derivation from weakly stable and unstable sets were
finally proposed in Hyeraci and Topputo (2010). To date, despite the effort
put in numerous works (García and Gómez, 2007; Topputo et al, 2008; Belbruno
et al, 2008, 2010, 2013), both WSB and BC are still not completely understood.
Nonetheless, a connection between celestial and quantum mechanics was recently
found exploiting the WSB (Belbruno, 2020), providing a fresh perspective to
tackle the problem.

BC orbits are characterized by ICs escaping the target when integrated
backward and performing $n$ revolutions about it when propagated forward,
neither impacting or escaping the target. In forward time, particles flying on
BC orbits approach the target coming from outside its sphere of influence and
remain temporarily captured about it. After a certain time, the particle
escapes if an energy dissipation mechanism does not take place to make the
capture permanent. To dissipate energy either a breaking maneuver or the
target’s atmosphere (if available) could be used (Luo and Topputo, 2021). In
this work, BC sets for comparison purposes are derived propagating the EoM in
Eq. (1) and following the procedure in Luo et al (2014).

When searching for BC opportunities, most of the trajectories found are
spurious solutions which are typically not useful for mission design purposes
(Luo et al, 2014). Useful solutions are detected exploiting the regularity
index666In previous works this was referred to as stability index (Luo et al,
2014; Luo and Topputo, 2015, 2017). However, in Dei Tos et al (2018), the
adjustment from _stability_ to _regularity_ index was proposed to avoid
misunderstandings with the periodic orbit stability index. The same
nomenclature introduced in Dei Tos et al (2018) is used in this work. $S$ and
regularity coefficient $\Delta S_{\%}$ (Dei Tos et al, 2018). The aim is
seeking for ideal orbits that presents regular post-capture legs resulting in
$n$ revolutions about the target which are similar in orientation and shape.
Numerical experiments showed that high-quality post-capture orbits are
identified by small regularity index and coefficient (Dei Tos et al, 2018; Luo
and Topputo, 2017, 2015; Luo et al, 2014). If the regularity index and
coefficient are indicators used to qualitatively judge post-capture legs,
capture occurrence is quantitatively measured through the capture ratio
$\mathcal{R}_{\mathcal{C}}$ (Luo and Topputo, 2015). Typically, search spaces
characterized by larger capture ratio are desirable when looking for BC
orbits.

Effects on BC by gravitational attractions of many bodies besides the
primaries and SRP have been investigated in previous works (Merisio and
Topputo, 2021; Aguiar and Topputo, 2018; Luo and Topputo, 2017, 2015), with
Merisio and Topputo (2021), and Aguiar and Topputo (2018) also considering the
target non-spherical gravity perturbations. Compared to the restricted three-
body problem, the $N$-body model is more adequate for constructing ballistic
capture orbits as proved in Luo and Topputo (2015), with a particle being more
easily captured when considering additional gravitational attractions. The
latter was later confirmed in Luo and Topputo (2017), where the presence of
moons is demonstrated to increase the capture ratio
$\mathcal{R}_{\mathcal{C}}$ and improves the regularity of post-capture orbits
while accommodating larger pre-capture energies. Similarly, SRP increases the
chances of being temporarily captured about the target and regularizes the
post-capture legs (Merisio and Topputo, 2021; Aguiar and Topputo, 2018). Even
if not duly discussed, the aforementioned remarks are observed also in the
results of this work. Overall, supplementary terms in the dynamics seem to
favor the manifestation of the BC phenomenon.

####  2.2.1 Definitions of particle stability and sub-sets

A particle stability is inferred using a plane in the three-dimensional
physical space (Belbruno and Miller, 1993), according to the spatial stability
definition provided in Luo et al (2014). The following indications are used to
classify stability, see Luo et al (2014) for more details:  1) a particle
completes a revolution around the target according to Remark 1 and Eq. (5) in
Luo et al (2014);  2) a particle escapes from the target according to Remark 2
and Eq. (6) in Luo et al (2014);  3) a particle impacts with the target
according to Remark 3 and Eq. (7) in Luo et al (2014).  Consistent variants of
Eq. (7) in Luo et al (2014) can be derived to locate impacts with target’s
moons, if present.

Based on its dynamical behavior, a propagated trajectory is said to be:  i)
weakly stable(sub-set $\mathcal{W}_{i}$) if the particle performs $i$ complete
revolutions around the target without escaping or impacting with it or its
moons;  ii) unstable(sub-set $\mathcal{X}_{i}$) if the particle escapes from
the target before completing the $i$-th revolution;  iii) target–crash(sub-set
$\mathcal{K}_{i}$) if the particle impacts with the target before completing
the $i$-th revolution;  iv) moon–crash(sub-set $\mathcal{M}_{i}$) if the
particle impacts with one of the target’s moons before completing the $i$-th
revolution;  v) acrobatic(sub-set $\mathcal{D}_{i}$) if none of the previous
conditions occurs within the integration time span.  Conditions ii)-v) apply
after the particle performs $(i-1)$ revolutions around the target. The sub-
sets are defined for $i\in\mathbb{Z}\textbackslash\\{0\\}$, where the sign of
$i$ informs on the propagation direction. When $i>0$ ($i<0$) the IC is
propagated forward (backward) in time. The overall domain, union of all sub-
sets, is defined as $\Omega$. A graph clarifying the relations between sub-
sets is shown in Fig. 1. A capture set is defined as
$\mathcal{C}_{-1}^{n}:=\mathcal{W}_{n}\cap\mathcal{X}_{-1}$. Therefore, it is
the intersection between the stable set in forward time $\mathcal{W}_{n}$ and
the unstable set in backward time $\mathcal{X}_{-1}$ (Luo et al, 2014).

Backward Forward
$\Omega$$\mathcal{W}_{1}$$\mathcal{W}_{2}$$\mathcal{X}_{2}$$\mathcal{K}_{2}$$\mathcal{M}_{2}$$\mathcal{D}_{2}$$\mathcal{X}_{1}$$\mathcal{K}_{1}$$\mathcal{M}_{1}$$\mathcal{D}_{1}$$\Omega$$\mathcal{W}_{-1}$$\mathcal{W}_{-2}$$\mathcal{X}_{-2}$$\mathcal{K}_{-2}$$\mathcal{M}_{-2}$$\mathcal{D}_{-2}$$\mathcal{X}_{-1}$$\mathcal{K}_{-1}$$\mathcal{M}_{-1}$$\mathcal{D}_{-1}$

Figure 1: Sets relations according to the algorithmic definition of WSB.

###  2.3 Differential algebra

DA propagation consists of assimilating a function $f$ of $v$ variables,
contained in $C^{k+1}$, with $T_{f}^{(k)}$ the Taylor expansion of $f$ at
order $k$ (Berz, 1999). The computation of such polynomials can be performed
efficiently, and provides a representation of the function $f$ on all of its
domain. Moreover, the DA structure ensures that algebraic and functional
operations are well-defined, particularly for the numerical solving of
ordinary differential equations (Berz, 1992).

The main advantage of this method is that the polynomial map only needs to be
computed once, and then it is evaluated in an arbitrarily large number of
points. In other words, to perform a Monte-Carlo estimation with a sample of
size $S$, only one computation of the map is needed, followed by $S$
polynomial evaluations, while classic Monte-Carlo requires $S$ propagations,
see Armellin et al (2010). The DA engine used in this work is the differential
algebra core engine (DACE) 777Library available at:
https://github.com/dacelib/dace [last accessed Feb 1, 2022]., implemented by
Politecnico di Milano (Rasotto et al, 2016; Massari et al, 2018).

####  2.3.1 Automatic domain splitting

DA propagation allows reducing the approximation error by increasing the order
$k$ of the polynomial mapping. Nevertheless, increasing the order leads to an
important growth of the computational time. Thus, the introduction of ADS, see
Wittig et al (2015), and Pirovano et al (2021) for more details.

The farther a point is from the constant part of the Taylor expansion, the
higher the loss of accuracy. The ADS algorithm controls this error by dividing
the domain into halves until the accuracy is below a predetermined tolerance.
Therefore, each new sub-domain is represented by its polynomial map causing a
much smaller approximation error, at a controlled computational cost (Wittig
et al, 2015). The ADS creates a division of the initial domain separating
areas with different behaviors from one another. The ADS algorithm requires
the passing of two parameters, in addition to the order of the DA: the
tolerance, and the maximum number of splits allowed. The former being the
approximation error threshold before the splitting of a domain occurs, while
the latter being the number of times the ADS routine can be applied to a sub-
domain.

##  3 Methodology

To map BC sets on the search space using DA propagation, a novel
classification algorithm is devised and herewith presented. To be successful,
the algorithm requires an accurate revolution period estimation. Consistency
and quality criteria are defined to evaluate performances of the resulting
mapping. Finally, fresh representation methods to visualize properly the large
resulting amount of data are exposed.

The search space is chosen to maximize the capture ratio $\mathcal{R}_{c}$
based on the analysis reported in Luo and Topputo (2015). It is defined in the
Mars-centered RTN reference frame at capture epoch $t_{0}$ on December 9,
2023, at 00:45:18.363 (UTC). At that epoch, Mars’s true anomaly with respect
to the Sun is equal to $270\text{\,}\mathrm{deg}$, maximizing
$\mathcal{R}_{c}$ (Luo and Topputo, 2015). The selected plane is defined by
inclination $i=$0.6283\text{\,}\mathrm{rad}$$, and right ascension of the
ascending node $\mathrm{RAAN}=$0.6283\text{\,}\mathrm{rad}$$. That because,
according to Fig. 10 in Luo and Topputo (2015), such values maximize the
capture ratio for Mars. Sought trajectories have osculating eccentricity
$e=0.99$ (Topputo and Belbruno, 2015), and mean anomaly
$M=$0\text{\,}\mathrm{rad}$$ at the initial epoch $t_{0}$. If $R_{\mars}$ is
the radius of Mars in $\mathrm{km}$, then the search space on the plane
defined above is a circular crown centered at Mars, from radius
$R_{\mars}+$100\text{\,}\mathrm{km}$$ up to radius $5R_{\mars}$. Hence,

$\left(r_{p},\omega\right)\in\left[R_{\mars}+\text{
$100\text{\,}\mathrm{km}$},5\cdot R_{\mars}\right]\times\left(-\pi,\pi\right]$
(3)

with $r_{p}$ the radius of the periapsis, and $\omega$ the argument of the
periapsis.

The main difference between the point-wise mapping of sets performed by
GRAvity TIdal Slide (GRATIS) (Topputo et al, 2018) and the DA mapping
presented in this work is that the DA propagator does not allow to count
revolutions. It means the sub-domains cannot be classified in the sub-sets
defined in Section 2.2.1. New sets are defined to solve this problem.
Furthermore, instead of tracking revolutions geometrically, as in Luo et al
(2014), the proposed novel classification algorithm counts revolution periods.
Finally, a bridge between the two mapping methods is established.

###  3.1 Redefinition of sub-sets and classification algorithm

There is no direct method to use the definitions of Section 2.2.1 with DA.
This is because accurately counting the revolutions of a continuous set of
particles around the target cannot be performed as in Remark 1 and Eq. (5) in
Luo et al (2014). Instead of using this condition to count completed
revolutions, the number of period elapsed since epoch is used. Consequently,
$\mathcal{D}_{i}$ cannot be defined when using DA propagation. Moreover, if a
subdomain reached the minimum size allowed by the ADS algorithm, and it tries
to split again, the accuracy of the mapping of this sub-domain cannot be
guaranteed. Thus, these inconsistent sub-domains need to be ruled out from the
sub-sets, so that only the consistent ones are retained.

Therefore, the definitions of Section 2.2.1 are adapted to DA propagation as
follows:  i) inconsistent(sub-set $\mathcal{\hat{I}}_{i}$) if the sub-domain
performed the maximum number of splits allowed and tries to split again before
completing the $i$-th period;  ii) weakly stable(sub-set
$\mathcal{\hat{W}}_{i}$) if the sub-domain is consistent, and performs $i$
complete periods without escaping or impacting the target or its moons;  iii)
unstable(sub-set $\mathcal{\hat{X}}_{i}$) if the sub-domain is consistent, and
escapes from the target before completing the $i$-th;  iv) target–crash(sub-
set $\mathcal{\hat{K}}_{i}$) if the sub-domain is consistent, and impacts with
the target before completing the $i$-th period;  v) moon–crash(sub-set
$\mathcal{\hat{M}}_{i}$) if the sub-domain is consistent, and impacts with one
of the target’s moons before completing the $i$-th period.  Conditions i), and
ii)-v) apply after the particle performs $(i-1)$ periods around the target. As
in Section 2.2.1, the sub-sets are defined for
$i\in\mathbb{Z}\textbackslash\\{0\\}$ and the same considerations about
propagation direction still apply. Moreover, the domain that is consistent
after $i$ periods is defined as

$\hat{\Omega}_{i}=\Omega\textbackslash\bigcup_{j=1}^{|i|}\hat{\mathcal{I}}_{\mathrm{sgn}({i})j}.$
(4)

In Fig. 2, a graph reporting the relations between the sub-sets adapted to DA
propagation is shown. As usual, a capture set is defined in DA propagation as
$\hat{\mathcal{C}}_{-1}^{n}:=\hat{\mathcal{W}}_{n}\cap\hat{\mathcal{X}}_{-1}$.

Backward Forward
$\Omega$$\mathcal{\hat{W}}_{1}$$\mathcal{\hat{W}}_{2}$$\mathcal{\hat{X}}_{2}$$\mathcal{\hat{K}}_{2}$$\mathcal{\hat{M}}_{2}$$\mathcal{\hat{I}}_{2}$$\mathcal{\hat{X}}_{1}$$\mathcal{\hat{K}}_{1}$$\mathcal{\hat{M}}_{1}$$\mathcal{\hat{I}}_{1}$$\Omega$$\mathcal{\hat{W}}_{-1}$$\mathcal{\hat{W}}_{-2}$$\mathcal{\hat{X}}_{-2}$$\mathcal{\hat{K}}_{-2}$$\mathcal{\hat{M}}_{-2}$$\mathcal{\hat{I}}_{-2}$$\mathcal{\hat{X}}_{-1}$$\mathcal{\hat{K}}_{-1}$$\mathcal{\hat{M}}_{-1}$$\mathcal{\hat{I}}_{-1}$$\hat{\Omega}_{1}$$\hat{\Omega}_{2}$$\hat{\Omega}_{-1}$$\hat{\Omega}_{-2}$

Figure 2: Sets relations according to classification with DA and ADS.

###  3.2 Classification algorithm

The classification algorithm shown in Algorithm 1 is a while loop divided in
two parts:

  1. 1.

the DA propagation of the stable set $\mathcal{\hat{W}}_{i-1}$ from the
current period to the next one, thus from $T_{i-1}$ to $T_{i}$;

  2. 2.

the extraction of all the sub-domains to classify them in the right set.

A bridge is built between the definitions of Section 2.2.1 and the ones of
Section 3.1. The similarities are

$\hat{\Omega}_{i}\rightarrow\Omega,\quad\mathcal{\hat{I}}_{i}\rightarrow\emptyset,\quad\mathcal{\hat{K}}_{i}\rightarrow\mathcal{K}_{i},\quad\mathcal{\hat{M}}_{i}\rightarrow\mathcal{M}_{i},\quad\mathcal{\hat{X}}_{i}\rightarrow\mathcal{X}_{i},\quad\mathcal{\hat{W}}_{i}\rightarrow\mathcal{W}_{i},\quad\mathcal{\hat{C}}_{-1}^{i}\rightarrow\mathcal{C}_{-1}^{i},\quad\emptyset\rightarrow\mathcal{D}_{i}.$
(5)

Note that the inconsistent set $\mathcal{\hat{I}}_{i}$ is mapped to the empty
set due to the inability to compute these sub-domains with accuracy. Moreover,
the acrobatic set $\mathcal{D}_{i}$ is not mapped by the DA classification
algorithm. Nonetheless, this region of the search space represents a small
fraction of the total for long enough propagation times (Luo et al, 2014).

Set either $i=1$ (forward propagation) or $i=-1$ (backward propagation);

Set the search space $\mathcal{\hat{W}}_{0}$ $=\Omega$;

Set capture epoch $T_{0}=t_{0}$;

Set maximum number of periods $n_{\mathrm{max}}$;

while  _$|i|\leq n_{\mathrm{max}}$_ do

Propagate $\mathcal{\hat{W}}_{i-\mathrm{sgn}({i})}$ from
$T_{i-\mathrm{sgn}({i})}$ to $T_{i}$;

Extract inconsistent sub-domains in $\mathcal{\hat{I}}_{i}$;

Extract crash sub-domains in $\mathcal{\hat{K}}_{i}$;

Extract moon-crash sub-domains in $\mathcal{\hat{M}}_{i}$;

Extract escaped sub-domains in $\mathcal{\hat{X}}_{i}$;

Retain remaining sub-domains as $\mathcal{\hat{W}}_{i}$;

Set $i=i+\mathrm{sgn}({i})$;

end while

Algorithm 1 Classification algorithm.

###  3.3 Propagation time span

Since it is not possible to track revolutions geometrically, they are tracked
by counting the revolution periods. Thus, the need to estimate them with
fidelity. These periods are determined by a least-square regression on data
from point-wise computations issued by GRATIS 888The time regressions have
been derived from ICs propagated taking into account also the gravitational
attractions of Uranus (B), and Neptune (B), later discarded due to their
negligible influence., with respect to the radius of the periapsis $r_{p}$.
Two regression shapes are chosen:

  1. 1.

a square root shape

$T_{\sqrt{\cdot}}=A+B\sqrt{r_{p}},\ \text{with }A=-0.19329,\ B=2.10555;$ (6)

  2. 2.

a logarithmic shape

$T_{\ln{(\cdot)}}=A+B\ln{(r_{p})},\ \text{with }A=-7.35410,\ B=1.44254.$ (7)

Fig. 3 shows these regressions for 2 and 6 periods, compared to the dataset
generated by GRATIS.

(a) After 2 revolutions.

(b) After 6 revolutions.

Figure 3: Regression of the revolution periods as a function of the radius of
the periapsis $r_{p}$. Data from point-wise computations issued by GRATIS.

For the rest of this work, the logarithmic shape is chosen over the square
root one, due to a better fit. Nevertheless, the regression does not fit well
for 2 revolutions at large $r_{p}$, due to the high concentration of points at
small radii, as shown by Fig. 4. Moreover, Fig. 5 highlights how the
distribution of the revolution times is widely spread, with respect to
$r_{p}$, as the whiskers represent the minimum and maximum values, with the
median in the middle. Therefore, the regressions displayed in Eqs. (6) and (7)
make up for a strong hypothesis on the behavior on the revolution periods.

(a) After 2 revolutions.

(b) After 6 revolutions.

Figure 4: Histogram of the periapsis $r_{p}$. Data from point-wise
computations issued by GRATIS.

(a) After 2 revolutions.

(b) After 6 revolutions.

Figure 5: Distribution of the revolution periods as a function of the radius
of the periapsis $r_{p}$. Data from point-wise computations issued by GRATIS.

###  3.4 Consistency and quality criteria

Once the classification of the search space is performed, tools need to be
developed to evaluate it. Two criteria are defined to do so. The first is the
consistency criterion, which assesses the parts of the domain where the DA
mapping cannot be trusted. The second is the quality criterion, which assesses
the performances of the DA mapping with respect to point-wise reference
mapping carried out with GRATIS.

##### Consistency criterion.

The consistency is a number assigned to each sub-domain. After the
classification algorithm is performed, each sub-domain is given a consistency
of either 1 if it does not belong to an inconsistent set
$\mathcal{\hat{I}}_{i}$, or 0 if it does. Consequently, the default
consistency of a sub-domain is 1. Yet, if a sub-domain reaches the minimum
size allowed and tries to split again, it is set to 0. The main interest is to
compute the global consistency on all of the search space. It can be computed
by performing the mean of the consistencies, pondered by the size of each sub-
domain. The global consistency will then represent the ratio of the mapping
that is not inconsistent. Therefore, it is the proportion of the mapping where
ADS guarantees the accuracy of the mapping.

##### Quality criterion.

The quality criterion is a value assigned to each set issued from GRATIS’s
classification. It represents the proportion of the set from GRATIS that is
well-mapped by the DA classification algorithm. To be computed, the quality
criterion requires a sample propagated point-wise with GRATIS. The quality
$q_{\mathcal{A}}$ of a GRATIS set $\mathcal{A}$ is computed as follows:

  1. 1.

for each point $x\in\mathcal{A}$, if $x$ belongs to the set
$\mathcal{\hat{A}}$, according to the bridge between DA and point-wise
classification exposed in Eq. (5), then the quality of that point $q_{x}=1$,
otherwise $q_{x}=0$;

  2. 2.

the quality $q_{\mathcal{A}}$ of a set $\mathcal{A}$ is the mean of all the
values of $\left(q_{x}\right)_{x\in\mathcal{A}}$. In other words

$q_{\mathcal{A}}=\frac{\\#\left(\mathcal{A}\cap\mathcal{\hat{A}}\right)}{\\#\left(\mathcal{A}\right)}.$
(8)

In the case where $\\#\left(\mathcal{A}\right)=0$, then
$\\#\left(\mathcal{A}\cap\mathcal{\hat{A}}\right)=0$. Thus, by convention
$q_{\mathcal{A}}=1$;

  3. 3.

the confidence intervals on $q_{\mathcal{A}}$ are evaluated according to
Robert and Casella (2004), and Hanley and Lippman-Hand (1983).

The quality represents the probability for a point from the point-wise mapping
to be mapped correctly with the DA mapping. Note that since
$\mathcal{\hat{I}}_{i}$ is mapped to $\emptyset$, the quality criterion of
$\Omega_{i}$ is less than or equal to the consistency criterion after $i$
revolutions. A schematic representation clarifying the meaning of the quality
criterion is shown in Fig. 6.

$\hat{\mathcal{A}}$$\times$$\mathcal{A}$$\circ$$\mathcal{B}$$\hat{\mathcal{B}}$$\times$$\mathcal{B}$$\times$$\mathcal{B}$$\hat{\mathcal{C}}$$\circ$$\mathcal{A}$$\times$$\mathcal{C}$$\circ$$\mathcal{B}$$\times$$\mathcal{C}$Legend:Sub-
domain classified by DA$\times$GRATIS IC accurately mapped by DA$\circ$GRATIS
IC inaccurately mapped by
DA$q_{\mathcal{A}}=\frac{1}{2}=0.5$$q_{\mathcal{B}}=\frac{2}{4}=0.5$$q_{\mathcal{C}}=\frac{2}{2}=1$$\\#\left(\mathcal{D}\right)=0\implies
q_{\mathcal{D}}=1$ Figure 6: Example of computation of the quality criterion
for a given mapping.

###  3.5 Representation methods

The representation of the results produced by this methodology raises two
issues. Firstly, due to ADS, the largest fully-split sub-domains have a size
of $\approx$50\text{\,}\mathrm{km}$$, while the typical length of the search
space is $\approx$3000\text{\,}\mathrm{km}$$. Thus, it is impossible to see
them or smaller sub-domains on a global representation. To solve this problem,
the resolution of the display is downgraded to a smaller one. Furthermore,
instead of plotting each sub-domain independently, the density of sub-domains
per pixel is preferred. It allows detecting areas with a high number of sub-
domains even on a global visualization.

Secondly, the Cartesian representation of the search space tends to shrink
details for low values of the radius of the periapsis $r_{p}$, while the size
of structures located at larges values of $r_{p}$ is amplified. However, areas
with a larger amount of data to visualize are located at small $r_{p}$. Thus,
results are visualized on the $r_{p}\times\omega$ plane, which is the search
space in Keplerian coordinates.

##  4 Results

In this work, the results from two simulations with initial grids having
different resolution are presented. The order of the DA mappings is 20 to
provide good accuracy on consistent domains and delay the triggering of ADS.
Therefore, increasing the order lowered the computational time and enhanced
consistency. Nevertheless, a larger order would not increase the approximation
performances, see Wittig et al (2015). The first simulation relies on a coarse
grid providing a low-resolution mapping where the search space is initially
divided into $32\times 32$ domains and the maximum number of splits allowed to
the ADS algorithm is 9. On the contrary, the second is computed on a finer
grid returning a high-resolution mapping. In the latter, the search space is
divided into $128\times 128$ domains and 10 maximum splits are allowed.
Results are visualized in the low-resolution mapping in Cartesian coordinates.
Sub-domains densities and last step epochs are computed on the whole search
space. Then, the DA classification outcome is rendered for qualitative
analysis. Finally, a quantitative analysis of those mappings exploiting
consistency and quality criteria follows.

###  4.1 Mapping

Fig. 7 represents the low-resolution search space in Cartesian coordinates,
before and after the propagation with ADS (left and right, respectively).
Before propagation, the search space is divided regularly into small domains
for parallelization (see Fig. 7(a)). Conversely, an irregular sub-divisions
highlighting dynamical changes is shown in Fig. 7(b). This visualization shows
the necessity to use different visualization methods. Indeed, it is hard to
analyze close to zones where ADS creates a large number of domains, although
these are the most interesting regions.

(a) Before propagation with ADS.

(b) After propagation with ADS.

Figure 7: Cartesian visualization of the mapping of the search space.

####  4.1.1 Density of sub-domains and last step epochs

Fig. 8 shows the density of sub-domains for the low-resolution and high-
resolution mappings (left and right, respectively). They mostly differ because
the high-resolution mapping provides more variations. Indeed, the density is
either maximal (in grey) or minimal (in white) in the low-resolution one. On
the contrary, the high-resolution mapping provides more shades, highlighting
several degrees of nonlinearity accurately mapped. In addition, the low-
resolution mapping created $7.2\times 10^{4}$ sub-domains. It is fewer
propagations than the point-wise cartography of Section 4.1.2 for the
representation of the same search space. Meanwhile, the number of sub-domains
in the high-resolution mapping is $1.2\times 10^{6}$, thus, more than ten
times the number of point-wise propagations. Therefore, there is no
computational benefit in using high-resolution mapping instead of point-wise
propagations. However, as opposed to a point-wise mapping, the sub-domains are
heterogeneously distributed, and are mainly located in highly nonlinear
regions. For instance, the vicinity of Mars and areas that lead to collisions
have a large sub-domains density compared to the majority of the search space.
Furthermore, the sub-domain density draws the boundaries between different
sets, due to their different dynamical behaviors. On that account, DA mapping
provides useful macroscopic information on the search space, while grid point-
wise mapping does not allow highlighting highly nonlinear areas or boundaries.

Moreover, epochs of last steps carried out by the propagation scheme are shown
in Fig. 9, in percent of the overall propagation time span. The integration of
a sub-domain stops either after a collision occurs or when the sub-domain
itself is declared inconsistent. Note that at low $r_{p}$, the lightest zones
on the low-resolution mapping (Fig. 9(a)) become completely white when
visualized in high resolution (Fig. 9(b)). Therefore, propagation of these
areas is now completed and sub-domains becomes consistent on a finer grid. In
these regions, the consistency criterion is expected to rise when the mapping
resolution increases.

(a) Low resolution.

(b) High resolution.

Figure 8: Density of sub-domains, in percent of the overall sum of sub-domains
after 6 periods.

(a) Low resolution.

(b) High resolution.

Figure 9: Last step epochs, in percent of the overall propagation time span
corresponding to 6 periods.

####  4.1.2 Classification

Fig. 10 shows the results of the classification algorithm after two
revolutions for the low-resolution mapping (left), the high-resolution mapping
(middle), and from point-wise propagation using GRATIS (right). The latter
used as a reference and derived through point-wise propagation from a sample
of $10^{5}$ points. Each color represents a different set, following the
bridge between sets from DA to point-wise mapping established in Section 3.2.
Thus, a set $\mathcal{A}$ from GRATIS classification is colored the same way
as $\mathcal{\hat{A}}$ from DA classification. Since inconsistent sets
$\mathcal{\hat{I}}_{i}$ and acrobatic sets $\mathcal{D}_{i}$ have not
corresponding sets in point-wise and DA mappings, respectively, they share the
same colors.

(a) Low resolution.

(b) High resolution.

(c) GRATIS.

Figure 10: Results of the classification algorithm after 2 revolutions
compared to GRATIS.

While these three mappings have a similar macroscopic look, they differ in
many ways. First of all, the sizes of the inconsistent sets decrease
dramatically from the low-resolution mapping to the high-resolution one.
Moreover, these gains in consistency mainly benefit the representation of the
stable set $\mathcal{W}_{2}$. Furthermore, the escape sets
$\mathcal{\hat{X}}_{1}$ and $\mathcal{\hat{X}}_{2}$ are poorly mapped compared
to the results of GRATIS on $\mathcal{X}_{1}$ and $\mathcal{X}_{2}$. These
issues in the mapping are mostly due to errors in the approximation of the
revolution period.

Fig. 11 presents the results of the classification algorithm after six
revolutions for the low-resolution mapping (left), the high-resolution mapping
(middle), and from point-wise propagation using GRATIS (right). Also in this
case the latter mapping is used as reference and it is derived through point-
wise propagation from a sample of $10^{5}$ points as before. These three
charts show simplified mappings. For the sake of clarity, all sets
$\left(\mathcal{A}_{j}\right)_{j\in\llbracket 1,6\rrbracket}$ are united under
the name $\mathcal{A}$, defined as

$\mathcal{A}=\bigcup_{j=1}^{6}{\mathcal{A}_{j}}.$ (9)

This representation method is adopted for all sets but $\mathcal{W}_{6}$, and
$\mathcal{\hat{W}}_{6}$.

(a) Low resolution.

(b) High resolution.

(c) GRATIS.

Figure 11: Results of the classification algorithm after 6 revolutions
compared to GRATIS.

The same observations can be done for Figs. 10 and 11. However, Fig. 11 shows
a strong predominance of the inconsistent set, especially at low resolution.
The high-resolution mapping improves these results but still struggles to
handle long accurate propagations on such a large search space. Nevertheless,
these results show that DA mapping with ADS can highlight changes of behavior
all over the search space. Moreover, DA mapping provides a continuous
description of the domain. The dynamical behavior is the same everywhere,
providing consistency. Meanwhile, point-wise computations deliver discrete
information challenging to interpolate due to nonlinearities.

###  4.2 Quantitative performance analysis

After performing a qualitative analysis on the resulting DA mappings, the two
criteria developed deliver quantitative data on the computed mappings.

####  4.2.1 Consistency criterion

In Fig. 12, the global consistency rate for both the low- and high-resolution
mappings after several revolution periods are shown in bright and dark colors,
respectively. The consistency at zero revolution is 100% since mappings only
represent ICs. It appears that the consistency remains high, even after six
revolutions, since it stays above 87%. It means that the size of the
inconsistent sets is 13% or less on the total search space. Moreover, the
consistency of the high-resolution mapping is around 5 points above the one of
the low-resolution one. It demonstrates the gain of accuracy delivered by
these additional computations.

Figure 12: Consistency criterion for several number of revolutions and the two
mapping resolutions.

####  4.2.2 Quality criterion

Fig. 13 represents the quality criterion for both low- and high-resolution
mappings on the left and right, respectively. Each group of bars corresponds
to a different collection of sets. For instance, $\mathcal{X}_{i}$ for
$i\in\llbracket 1,6\rrbracket$. Differently, each color represents a different
number of periods, from 1 (light) to 6 (dark). A bar’s height represents the
quality of that set in percent, while whiskers inform on the confidence
interval at 95% on such value. These estimations of the quality criterion are
computed against the sample of size $10^{5}$ obtained with GRATIS.

The quality of acrobatic sets $\mathcal{D}_{i}$ is not shown since they are
not mapped by the DA algorithm. In fact, their qualities are automatically set
to 0%. In addition, the quality of moon-crash sets $\mathcal{M}_{i}$ is always
100%. Indeed, among the $10^{5}$ propagated trajectories, none of them crashed
on Mars’ moons. The quality of these empty sets is then $100\%$ by convention.

(a) Low resolution.

(b) High resolution.

Figure 13: Quality criterion in percent of several sets and revolution
periods. Note that $\Omega_{i}$ is used to refer to the quality criteria of
$\Omega$ after $i$ periods.

The quality criterion is always higher for the high-resolution mapping than
for the low-resolution one. Especially for stable sets $\mathcal{W}_{i}$ and
capture sets $\mathcal{C}_{-1}^{i}$ since a larger part of these sets becomes
inconsistent during the low-resolution propagation, as shown by Fig. 11.
Furthermore, capture sets $\mathcal{C}_{-1}^{i}$ are mapped in a better or
equal way compared to stable sets $\mathcal{W}_{i}$. This is due to the
positioning of $\mathcal{C}_{-1}^{i}$ farther from Mars (Merisio and Topputo,
2021), than most of $\mathcal{W}_{i}$. Thus, as highlighted in Figs. 10 and
11, the consistency of this region is higher, which enables better quality.
However, designing a specific mission with ballistic capture requires much
more data than these two criteria and is beyond the scope of this work.
Moreover, this work does not provide any additional property on the identified
solutions, except that they belong to capture sets.

Except for those two groups of sets, only slight improvements are remarkable
when switching from low resolution to high resolution. Indeed, apart from
$\mathcal{X}_{1}$, escape sets $\mathcal{X}_{i}$ and crash sets
$\mathcal{K}_{i}$ fail to be accurately mapped, even though these regions of
the search space are consistent. Indeed, the proposed classification algorithm
does not track revolutions geometrically but uses period approximations.
Consequently, the classification implies major approximations. Nevertheless,
the poor mapping of these sets hardly impacts the quality of sets
$\mathcal{W}_{i}$ due to their small sizes. Therefore, the global quality
never drops below 80% in low resolution. The overall quality can be improved
at a high computational cost, compared to the point-wise mapping, thanks to
high-resolution mapping. Such a mapping improves the global classification
results by at least 4. These poor performances are mainly due to the inability
of the DA mapping engine to track revolutions geometrically instead of
temporally. Thus, DA mapping as presented in this work cannot outperform
point-wise mapping in terms of classification precision, but it exceeds point-
wise mapping concerning computational cost.

##  5 Conclusion

This article presented a methodology to build cartography of ballistic capture
sets using differential algebra mapping. It proposes an alternative
classification algorithm to sort the sub-domains produced by the automatic
domain splitting algorithm in newly defined capture sets. Moreover, instead of
tracking revolutions with geometrical methods, as in Luo et al (2014), they
are counted temporally, using regression on point-wise data. This work
establishes a bridge between the point-wise mapping and DA mapping, allowing
us to assess the performances of this new method. Hence, the introduction of
two criteria. The first one is the consistency criterion, which represents the
proportion of the search space where the mapping accuracy is guaranteed by the
ADS algorithm. The second one is the quality criterion, which represents the
success rate of the DA classification algorithm compared to the point-wise
one.

Results show that DA mapping of BC sets could be performed on large search
spaces, as ADS captured the dynamical variations on the whole domain.
Furthermore, the consistency criterion shows that more than 87% of the search
space is guaranteed by the ADS algorithm as accurate, even for low-resolution
mappings. Moreover, the quality criterion demonstrates that the global error
rate of DA mapping is below 20%. However, some small sets are either not or
poorly mapped using revolution period regression, mainly due to the incapacity
to count revolutions via geometrical information. This phenomenon does not
disappear when the resolution of the mapping is increased. Nonetheless, DA
mapping of BC sets can prove useful in various situations. For instance, it
can be used for fast mapping of a plane, to then target a restricted domain of
the search space. Either with a dense point-wise mapping or with a high-
resolution DA mapping of the restricted area of interest.

The data sets supporting the findings of this study are available on
Zenodo999https://zenodo.org/ [last accessed Feb 1, 2022] with the identifiers
https://doi.org/10.5281/zenodo.6103720 and
https://doi.org/10.5281/zenodo.5940101.

###### Acknowledgements.

T.C. would like to thank the Department of Aerospace Science and Technology
(DAER) of Politecnico di Milano for the host and the warm welcome, as well as
the Toulouse graduate School of Aerospace Engineering (TSAE) for the funding
of his stay at Politecnico di Milano. G.M. and F.T. would like to acknowledge
the European Research Council (ERC) since part of this work has received
funding from the ERC under the European Union’s Horizon 2020 research and
innovation programme (grant agreement No. 864697).

## References

  * Acton (1996) Acton C (1996) Ancillary data services of NASA’s navigation and ancillary information facility. Planetary and Space Science 44(1):65–70, doi:10.1016/0032-0633(95)00107-7
  * Acton et al (2018) Acton C, Bachman N, Semenov B, Wright E (2018) A look towards the future in the handling of space science mission geometry. Planetary and Space Science 150:9–12, doi:10.1016/j.pss.2017.02.013
  * Aguiar and Topputo (2018) Aguiar G, Topputo F (2018) A technique for designing Earth–Mars low-thrust transfers culminating in ballistic capture. In: 7th International Conference on Astrodynamics Tools and Techniques (ICATT), pp 1–8
  * Archinal et al (2011) Archinal B, A’Hearn M, Bowell E, Conrad A, Consolmagno G, Courtin R, Fukushima T, Hestroffer D, Hilton J, Krasinsky G, et al (2011) Report of the IAU working group on cartographic coordinates and rotational elements: 2009\. Celestial Mechanics and Dynamical Astronomy 109(2):101–135, doi:10.1007/s10569-010-9320-4
  * Armellin et al (2010) Armellin R, Di Lizia P, Bernelli-Zazzera F, Berz M (2010) Asteroid close encounters characterization using differential algebra: the case of Apophis. Celestial Mechanics and Dynamical Astronomy 107(4):451–470, doi:10.1007/s10569-010-9283-5
  * Baresi et al (2021) Baresi N, Fu X, Armellin R (2021) A high-order taylor polynomials approach for continuing trajectories in restricted three-body problems. In: 2020 AAS/AIAA Astrodynamics Specialist Conference (South Lake Tahoe, California), 09/08/2020 - 13/08/2020), pp 1–19
  * Belbruno (1987) Belbruno E (1987) Lunar capture orbits, a method of constructing Earth Moon trajectories and the lunar GAS mission. In: 19th International Electric Propulsion Conference, American Institute of Aeronautics and Astronautics, pp 1–9, doi:10.2514/6.1987-1054
  * Belbruno (2004) Belbruno E (2004) Capture Dynamics and Chaotic Motions in Celestial Mechanics. Princeton University Press, doi:10.1515/9780691186436
  * Belbruno (2020) Belbruno E (2020) Relation between solutions of the Schrödinger equation with transitioning resonance solutions of the gravitational three-body problem. Journal of Physics Communications 4(1):015,012, doi:10.1088/2399-6528/ab693f
  * Belbruno and Carrico (2000) Belbruno E, Carrico J (2000) Calculation of weak stability boundary ballistic lunar transfer trajectories. In: Astrodynamics Specialist Conference, p 4142, doi:10.2514/6.2000-4142
  * Belbruno and Miller (1990) Belbruno E, Miller J (1990) A ballistic lunar capture trajectory for the Japanese spacecraft Hiten. Tech. rep., Jet Propulsion Laboratory, IOM 312/904-1731-EAB
  * Belbruno and Miller (1993) Belbruno E, Miller J (1993) Sun-perturbed Earth-to-Moon transfers with ballistic capture. Journal of Guidance, Control, and Dynamics 16(4):770–775, doi:10.2514/3.21079
  * Belbruno et al (2008) Belbruno E, Topputo F, Gidea M (2008) Resonance transitions associated to weak capture in the restricted three-body problem. Advances in Space Research 42(8):1330–1351, doi:10.1016/j.asr.2008.01.018
  * Belbruno et al (2010) Belbruno E, Gidea M, Topputo F (2010) Weak stability boundary and invariant manifolds. SIAM Journal on Applied Dynamical Systems 9(3):1061–1089, doi:10.1137/090780638
  * Belbruno et al (2013) Belbruno E, Gidea M, Topputo F (2013) Geometry of weak stability boundaries. Qualitative Theory of Dynamical Systems 12(1):53–66, doi:10.1007/s12346-012-0069-x
  * Benkhoff et al (2021) Benkhoff J, Murakami G, Baumjohann W, Besse S, Bunce E, Casale M, Cremosese G, Glassmeier K-H, Hayakawa H, Heyner D, et al (2021) BepiColombo-mission overview and science goals. Space Science Reviews 217(8):1–56, doi:10.1007/s11214-021-00861-4
  * Berz (1992) Berz M (1992) High-order computation and normal form analysis of repetitive systems. In: AIP Conference Proceedings, pp 1–36, doi:10.1063/1.41975
  * Berz (1999) Berz M (1999) Modern Map Methods in Particle Beam Physics. Elsevier, doi:10.1016/s1076-5670(08)70227-1
  * Caleb and Lizy-Destrez (2021) Caleb T, Lizy-Destrez S (2021) Can uncertainty propagation solve the mysterious case of Snoopy? In: Vasile M, Quagliarella D (eds) Advances in Uncertainty Quantification and Optimization Under Uncertainty with Aerospace Applications, Springer International Publishing, pp 109–128, doi:10.1007/978-3-030-80542-5_8
  * Chung et al (2010) Chung M, Hatch S, Kangas J, Long S, Roncoli R, Sweetser T (2010) Trans-lunar cruise trajectory design of GRAIL (Gravity Recovery and Interior Laboratory) mission. In: AIAA/AAS Astrodynamics Specialist Conference, p 8384, doi:10.2514/6.2010-8384
  * Circi and Teofilatto (2001) Circi C, Teofilatto P (2001) On the dynamics of weak stability boundary lunar transfers. Celestial Mechanics and Dynamical Astronomy 79(1):41–72, doi:10.1023/A:1011153610564
  * Dei Tos et al (2018) Dei Tos D, Russell R, Topputo F (2018) Survey of Mars ballistic capture trajectories using periodic orbits as generating mechanisms. Journal of Guidance, Control, and Dynamics 41(6):1227–1242, doi:10.2514/1.g003158
  * Di Lizia et al (2008) Di Lizia Pierluigi, Armellin R, Lavagna M (2008) Application of high order expansions of two-point boundary value problems to astrodynamics. Celestial Mechanics and Dynamical Astronomy 102(4):355–375, doi:10.1007/s10569-008-9170-5
  * García and Gómez (2007) García F, Gómez G (2007) A note on weak stability boundaries. Celestial Mechanics and Dynamical Astronomy 97(2):87–100, doi:10.1007/s10569-006-9053-6
  * Hairer et al (1993) Hairer E, Nørsett S, Wanner G (1993) Solving Ordinary Differential Equations I. Springer Berlin Heidelberg, doi:10.1007/978-3-540-78862-1
  * Haller (2011) Haller G (2011) A variational theory of hyperbolic lagrangian coherent structures. Physica D: Nonlinear Phenomena 240(7):574–598, doi:10.1016/j.physd.2010.11.010
  * Hanley and Lippman-Hand (1983) Hanley J, Lippman-Hand A (1983) If nothing goes wrong, is everything all right? Journal of American Medical Association (JAMA) 249(13):1743, doi:10.1001/jama.1983.03330370053031
  * Hyeraci and Topputo (2010) Hyeraci N, Topputo F (2010) Method to design ballistic capture in the elliptic restricted three-body problem. Journal of guidance, control, and dynamics 33(6):1814–1823, doi:10.2514/1.49263
  * Ivashkin (2002) Ivashkin V (2002) On trajectories of Earth–Moon flight of a particle with its temporary capture by the Moon. In: Doklady Physics, Springer, vol 47, pp 825–827, doi:10.1134/1.1526433
  * Luo and Topputo (2015) Luo Z-F, Topputo F (2015) Analysis of ballistic capture in Sun–planet models. Advances in Space Research 56(6):1030–1041, doi:10.1016/j.asr.2015.05.042
  * Luo and Topputo (2017) Luo Z-F, Topputo F (2017) Capability of satellite-aided ballistic capture. Communications in Nonlinear Science and Numerical Simulation 48:211–223, doi:10.1016/j.cnsns.2016.12.021
  * Luo and Topputo (2021) Luo Z-F, Topputo F (2021) Mars orbit insertion via ballistic capture and aerobraking. Astrodynamics 5(2):167–181, doi:10.1007/s42064-020-0095-4
  * Luo et al (2014) Luo Z-F, Topputo F, Bernelli-Zazzera F, Tang G (2014) Constructing ballistic capture orbits in the real Solar System model. Celestial Mechanics and Dynamical Astronomy 120(4):433–450, doi:10.1007/s10569-014-9580-5
  * Manzi and Topputo (2021) Manzi M, Topputo F (2021) A flow-informed strategy for ballistic capture orbit generation. Celestial Mechanics and Dynamical Astronomy 133(11-12):1–16, doi:10.1007/s10569-021-10048-2
  * Massari et al (2018) Massari M, Di Lizia P, Cavenago F, Wittig A (2018) Differential algebra software library with automatic code generation for space embedded applications. In: 2018 AIAA Information Systems-AIAA Infotech @ Aerospace, American Institute of Aeronautics and Astronautics, pp 1–13, doi:10.2514/6.2018-0398
  * Merisio and Topputo (2021) Merisio G, Topputo F (2021) Characterization of ballistic capture corridors aiming at autonomous ballistic capture at Mars. In: 2021 AIAA/AAS Astrodynamics Specialist Conference, pp 1–21
  * Montenbruck and Gill (2000) Montenbruck O, Gill E (2000) Satellite Orbits Models, Methods and Applications. Springer, doi:10.1007/978-3-642-58351-3
  * Park et al (2021) Park R, Folkner W, Williams J, Boggs D (2021) The JPL planetary and lunar ephemerides DE440 and DE441. The Astronomical Journal 161(3):105, doi:10.3847/1538-3881/abd414
  * Pirovano et al (2021) Pirovano L, Armellin R, Siminski J, Flohrer T (2021) Differential algebra enabled multi-target tracking for too-short arcs. Acta Astronautica 182:310–324, doi:10.1016/j.actaastro.2021.02.023
  * Prince and Dormand (1981) Prince P, Dormand J (1981) High order embedded Runge-Kutta formulae. Journal of computational and applied mathematics 7(1):67–75, doi:10.1016/0771-050x(81)90010-3
  * Racca et al (2002) Racca G, Marini A, Stagnaro L, Van Dooren J, Di Napoli L, Foing B, Lumb R, Volp J, Brinkmann J, Grünagel R, et al (2002) SMART-1 mission description and development status. Planetary and space science 50(14-15):1323–1337, doi:10.1016/S0032-0633(02)00123-X
  * Rasotto et al (2016) Rasotto M, Morselli A, Wittig A, Massari M, Di Lizia P, Armellin R, Yabar Valles C, Urbina Ortega C (2016) Differential algebra space toolbox for nonlinear uncertainty propagation in space dynamics. In: 6th International Conference on Astrodynamics Tools and Techniques (ICATT), International Conference on Astrodynamics Tools and Techniques, pp 1–11
  * Robert and Casella (2004) Robert C, Casella G (2004) Monte Carlo Statistical Methods. Springer New York, doi:10.1007/978-1-4757-4145-2
  * Scheeres et al (2011) Scheeres D, Rosengren A, McMahon J (2011) The dynamics of high area-to-mass ratio objects in earth orbit: The effect of solar radiation pressure. In: Proceedings of the AAS/AIAA Space Flight Mechanics Meeting, number AAS, pp 11–178
  * Schuster and Jehn (2014) Schuster A, Jehn R (2014) Influence of the Mercury gravity field on the orbit insertion strategy of BepiColombo. Aerospace Science and Technology 39:546–551, doi:10.1016/j.ast.2014.06.003
  * Servadio and Zanetti (2021) Servadio S, Zanetti R (2021) Differential algebra-based multiple Gaussian particle filter for orbit determination. Journal of Optimization Theory and Applications 191:1–27, doi:10.1007/s10957-021-01934-8
  * Servadio et al (2022) Servadio S, Zanetti R, Armellin R (2022) Maximum a posteriori estimation of Hamiltonian systems with high order Taylor polynomials. The Journal of the Astronautical Sciences 69(2):511–536, doi:10.1007/s40295-022-00304-4
  * Sousa Silva and Terra (2012) Sousa Silva P, Terra M (2012) Applicability and dynamical characterization of the associated sets of the algorithmic weak stability boundary in the lunar sphere of influence. Celestial Mechanics and Dynamical Astronomy 113(2):141–168, doi:10.1007/s10569-012-9409-z
  * Topputo and Belbruno (2009) Topputo F, Belbruno E (2009) Computation of weak stability boundaries: Sun–Jupiter system. Celestial Mechanics and Dynamical Astronomy 105(1-3):3–17, doi:10.1007/s10569-009-9222-5
  * Topputo and Belbruno (2015) Topputo F, Belbruno E (2015) Earth–Mars transfers with ballistic capture. Celestial Mechanics and Dynamical Astronomy 121(4):329–346, doi:10.1007/s10569-015-9605-8
  * Topputo et al (2008) Topputo F, Belbruno E, Gidea M (2008) Resonant motion, ballistic escape, and their applications in astrodynamics. Advances in Space Research 42(8):1318–1329, doi:10.1016/j.asr.2008.01.017
  * Topputo et al (2018) Topputo F, Dei Tos D, Mani K, Ceccherini S, Giordano C, Franzese V, Wang Y (2018) Trajectory design in high-fidelity models. In: 7th International Conference on Astrodynamics Tools and Techniques (ICATT), pp 1–9
  * Topputo et al (2021) Topputo F, Wang Y, Giordano C, Franzese V, Goldberg H, Perez-Lissi F, Walker R (2021) Envelop of reachable asteroids by M-ARGO CubeSat. Advances in Space Research 67(12):4193–4221, doi:10.1016/j.asr.2021.02.031
  * Valli et al (2013) Valli M, Armellin R, Di Lizia P, Lavagna M (2013) Nonlinear mapping of uncertainties in celestial mechanics. Journal of Guidance, Control, and Dynamics 36(1):48–63, doi:10.2514/1.58068
  * Wittig et al (2015) Wittig A, Di Lizia P, Armellin R, Makino K, Bernelli-Zazzera F, Berz M (2015) Propagation of large uncertainty sets in orbital dynamics by automatic domain splitting. Celestial Mechanics and Dynamical Astronomy 122(3):239–261, doi:10.1007/s10569-015-9618-3
  * Wittig et al (2017) Wittig A, Colombo C, Armellin R (2017) Long-term density evolution through semi-analytical and differential algebra techniques. Celestial Mechanics and Dynamical Astronomy 128(4):435–452, doi:10.1007/s10569-017-9756-x