Thesis Patricia Egger

Ecole Polytechnique Fédérale de
Lausanne
&
Massachussetts Institute of
Technology
Master thesis
Efficient computational approaches

to trajectory prediction for future
Earth-asteroid impact mitigation
Patricia Egger Under the supervision of:

pegger@mit.edu Prof. Marco Picasso
patricia.egger@epfl.ch
and
Sciper No : 184951 Prof. Olivier de Weck
July 14, 2014

Abstract
The problem that we will study in this thesis is that of predicting trajec-
tories of asteroids within our solar system. This thesis is a component of a
MIT PhD thesis that aims at determining the utility of precursor missions
in asteroid deflection campaigns.
The main contributions of this thesis are (i) the retracing of the historical
development of the equations of motion used in asteroid and planetary tra-
jectory propagation, (ii) the development of a propagator, PAT2 , through
validation of the numerical integrator ode113 built-in to MATLAB with er-
ror validation against NASA’s HORIZONS system and a short discussion
on why symplectic solvers will not be considered in this context, (iii) the
rescaling of the problem in order to reduce numerical noise and (iv) the ap-
plication of this propagator to five asteroids with very different orbits.
One of the main advantages and key differences with most available prop-
agators is that PAT2 only requires the initial conditions of the bodies (i.e.
the state vector y at t = 0, y0 ). In fact, there is no need to store previously
computed planetary ephemerides as is the case in HORIZONS, for example.
This method could produce bigger numerical errors compared to the latter
because the propagated planetary orbits will also include some errors. How-
ever, this is a major simplification for the user and furthermore allows for
the study of planetary trajectory propagation as well as that for asteroids.
The asteroids we consider are chosen because of the differences in their or-
bits’ characteristics, i.e. in order to observe the behavior of PAT2 relative
to important physical parameters. Using this propagator on the notorious
asteroid Apophis, we obtain a maximal numerical error of 500 kilometers
after 10 years. Results for the other asteroids studied range from 600 kilo-
meters to 2.5 thousand kilometers after 10 years. Runtimes range from less
than 7 minutes to about 28 minutes.
As the tradeoff between accuracy, or error, and runtime is essential in many
applications of the N-body problem, we present some data that will help
in choosing the optimal set-up for the propagation of an asteroid given its
physical parameters as well as the runtime and accuracy requirements.
i
Acknowledgements
There are a lot of people that I would like to thank for my experience
researching and writing this thesis. First, I would like express my very
great appreciation for my two supervisors, Professor Picasso from EPFL in
Switzerland and Professor de Weck from MIT in the United States. They
have allowed me to pursue a dream that all young engineering students
share: that of studying in the world’s greatest school. Both of them have
been very supportive of my work and helped me through some important
obstacles. A second special thank you goes to Professor de Weck, a truly
inspirational man and a great role model, for providing me with unimag-
inable opportunities. I would also like to thank my student mentor, PhD
candidate Sung Wook Paek, for introducing me to the world of asteroids and
helping me in my exploration of the world of astronautics. Another thank
you goes to Paul Chodas for not only giving me an inside look at asteroid
and planetary ephemerides at NASA, but also for allowing me to discover
JPL in a front row seat.
I would also like to extend my thanks to all the people from room 33-409
- Andrew, David, Ioana, Koki, Margaret, Narek, Paul, Roi, Sreeja, Sydney
and Takuto from SERG and Dani, Íñigo, Marc, Morgan and Peter from
SAL. They made my few months at MIT that much more fun. Dani and
Marc, thank you for taking me under your wing and Íñigo, thank you for the
spanish lessons. On the Swiss side, I would like to acknowledge my friends
Angelina and Charlotte for their never ending support as well as Alessandro,
Christoph and Francesco for their invaluable help throughout my Bachelor’s
and Master’s degrees at EPFL.
Finally, I would like to thank my family - my mother, father, brothers and
my better half, who have been there for me through the highs and lows, al-
ways encouraging me to reach for the stars (or asteroids, as the case may be).
Thanks so much to all of you! I will never forget this experience.
ii
Contents
1 Introduction and Motivation 1

2 Physical and Mathematical Background 5
2.1 The N-body problem . . . . . . . . . . . . . . . . . . . . . . . 6
2.1.1 Equations of motion . . . . . . . . . . . . . . . . . . . 7
2.1.2 Hamiltonian formulations . . . . . . . . . . . . . . . . 12
2.1.3 Derivation of the equations of motion . . . . . . . . . 13
2.2 Numerical Integration . . . . . . . . . . . . . . . . . . . . . . 16
2.2.1 Geometric numerical integration . . . . . . . . . . . . 16
2.2.2 Accuracy requirements . . . . . . . . . . . . . . . . . . 17
3 PAT2 : Propagator for Asteroid Trajectories Tool 19
3.1 Mathematical and Physical Model . . . . . . . . . . . . . . . 20
3.2 Rescaling . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 22
3.3 Numerical Integration Scheme . . . . . . . . . . . . . . . . . . 23
3.4 Output . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 26
4 Benchmark comparisons 27
4.1 Possible Candidates . . . . . . . . . . . . . . . . . . . . . . . 27
4.2 Benchmarking Results . . . . . . . . . . . . . . . . . . . . . . 30
5 Case Studies 31
5.1 Apophis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 33
5.2 Icarus . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 47
5.3 2007 FT3 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 51
5.4 2009 VZ39 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 52
5.5 2008 FF5 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 57
6 Conclusions and Future Work 61
Appendices 66
A JPL Visit Report 67
B Chebyshev Polynomial Interpolation 71
C How Asteroids Get Their Names 72
iii
List of Tables
1.1 Typical impact types, frequency and conesquences . . . . . . . . . 2
3.1 Positions, velocities and accelerations in units of seconds and

kilometers . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 22
3.2 Positions, velocities and accelerations in units of sidereal years
and hundred thousand kilometers (105 km) . . . . . . . . . . 22
3.3 Number of function evaluations for given tolerances . . . . . . 24
3.4 Error in the last time step for Apophis . . . . . . . . . . . . . 25
4.1 Benchmark comparisons . . . . . . . . . . . . . . . . . . . . . 30
5.1 Some of asteroid Apophis’ relevant physical data . . . . . . . 33

5.2 Error in the last time step for Apophis . . . . . . . . . . . . . 35
5.3 Error in the position of Apophis using the Newtonian equa-
tions of motion . . . . . . . . . . . . . . . . . . . . . . . . . . 37
5.4 Some of asteroid Icarus’ relevant physical data . . . . . . . . 47
5.5 Error in the last time step for Icarus . . . . . . . . . . . . . . 48
5.6 Some of asteroid 2007 FT3’s relevant physical data . . . . . . 51
5.7 Error in the last time step for 2007 FT3 . . . . . . . . . . . . 52
5.8 Some of asteroid 2009 VZ39’s relevant physical data . . . . . 54
5.9 Error in the last time step for 2009 VZ39 . . . . . . . . . . . 55
5.10 Number of observations used in the statistical fit and solution
number . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 56
5.11 Some of asteroid 2008 FF5’s relevant physical data . . . . . . 57
5.12 Error in the last time step for 2008 FF5 . . . . . . . . . . . . 58
6.1 Summary of results . . . . . . . . . . . . . . . . . . . . . . . . 63
iv
List of Figures
1.2 Differences between orbits of a typical NEA and PHA . . . . 3
2.2 Illustration of the ecliptic plane . . . . . . . . . . . . . . . . . 7

2.3 The asteroid belt . . . . . . . . . . . . . . . . . . . . . . . . . 10
3.2 Flow chart describing the inner workings of PAT2 . . . . . . . 21
4.2 Flow chart describing the inner workings of HORIZONS . . . 28
5.2 Orbits of the asteroids considered for the case studies . . . . 32

5.3 Orbit of Apophis . . . . . . . . . . . . . . . . . . . . . . . . . 34
5.4 Error in the position of Apophis . . . . . . . . . . . . . . . . 35
5.5 Component-wise error in the position of Apophis . . . . . . . 36
5.6 Normalized simulated total energy . . . . . . . . . . . . . . . 38
5.7 Ratio of the simulated potential over kinetic energy . . . . . . 39
5.8 Normalized total energy from HORIZONS . . . . . . . . . . . 40
5.9 Ratio of the potential over kinetic energy from HORIZONS . 41
5.10 Error in the position of the Sun . . . . . . . . . . . . . . . . . 43
5.11 Component-wise error in the position of the Sun . . . . . . . 43
5.12 Error in the position of Jupiter . . . . . . . . . . . . . . . . . 44
5.13 Component-wise error in the position of Jupiter . . . . . . . . 44
5.14 Error in the position of Mercury . . . . . . . . . . . . . . . . 45
5.15 Component-wise error in the position of Mercury . . . . . . . 45
5.16 Error in the position of Venus . . . . . . . . . . . . . . . . . . 46
5.17 Component-wise error in the position of Venus . . . . . . . . 46
5.18 Orbit of Icarus . . . . . . . . . . . . . . . . . . . . . . . . . . 48
5.19 Error in the position of Icarus . . . . . . . . . . . . . . . . . . 49
5.20 Points on Icarus’ orbit with biggest and smallest errors . . . . 50
5.21 Component-wise error in the position of Icarus . . . . . . . . 50
5.22 Orbit of 2007 FT3 . . . . . . . . . . . . . . . . . . . . . . . . 51
5.23 Error in the position of 2007 FT3 . . . . . . . . . . . . . . . . 53
5.24 Component-wise error in the position of 2007 FT3 . . . . . . 53
5.25 Orbit of 2009 VZ39 . . . . . . . . . . . . . . . . . . . . . . . . 54
5.26 Error in the position of 2009 VZ39 . . . . . . . . . . . . . . . 55
v
LIST OF FIGURES
5.27 Component-wise error in the position of 2009 VZ39 . . . . . . 56

5.28 Orbit of 2008 FF5 . . . . . . . . . . . . . . . . . . . . . . . . 58
5.29 Error in the position of 2008 FF5 . . . . . . . . . . . . . . . . 59
5.30 Points on 2008 FF5’s orbit with biggest errors . . . . . . . . . 59
5.31 Component-wise error in the position of 2008 FF5 . . . . . . 60
6.2 Tolerance versus runtime tradeoff. . . . . . . . . . . . . . . . . . 62

6.3 Tolerance versus numerical error tradeoff . . . . . . . . . . . . 62
vi
Glossary and acronyms
Aphelion: Defined as the point on a planet’s orbit where it is at

maximal distance from the Sun.
Asteroid: Also known as a minor planet, this is a large rocky body in

orbit around the Sun. Asteroids are believed to be ancient remnants of
the earliest years of the formation of our solar system. In comparison,
comets are a lot less dense than asteroids as they are mainly made of
ice and rock. This is why comets do not pose as much of a threat as
asteroids: they usually do not make it past our atmosphere.
AU (Astronomical Unit): Is the mean distance between the Earth

and the Sun. 1 AU = 149,597,870.700 kilometers.
Barycenter: Of a body, is its center of mass. Since the Sun makes

up about 99.8% of the mass of the entire solar system [20], the solar
system barycenter nearly coincides with that of the Sun.
Deflection mission: A technological solution for deviating the tra-

jectory of a potentially hazardous object so as to avoid it colliding with
the Earth.
∆V : Delta V: Change is velocity (denoted V ) needed to modify a

body’s trajectory. This can be understood as the amount of effort that
is needed to move a body from one trajectory to another.
Ecliptic, or ecliptic plane: Defined as the imaginary plane contain-

ing the Earth’s orbit around the Sun. The planetary bodies of our
solar system all lie approximately in this plane.
Ephemeris: Gives the positions of naturally occurring astronomical

objects as well as artificial satellites in the sky at a given time or times.
Ephemerides: Plural of the ephemeris.
Hamiltonian: Mathematical formalism to describe the equations of

motion of a physical system.
vii
LIST OF FIGURES
ICRF (International Celestial Reference Frame): With its ori-

gin at the solar system barycenter and "space fixed" axis directions, the
ICRF is meant to represent the most appropriate coordinate system
for expressing reference data on the positions and motions of celestial
objects.
NEO - respectively NEA (Near-Earth Object - resp. Near-

Earth Asteroid): These are objects that have a closest approach
distance with the Earth’s orbit of less than 0.3 AU.
Orbital eccentricity: Of a celestial object, is a real and positive

number, e, that determines the amount by which its orbit around the
center body deviates from a perfect circle. An orbit with e = 0 is
therefore a circle. When 0 < e < 1, we have an ellipse, when e = 1 a
parabola and e > 1 a hyperbola.
Perihelion: Defined as the point on a planet’s orbit where it is closest

to the Sun.
PHO - respectively PHA (Potentially Hazardous Object -

resp. Potentially Hazardous Asteroid): These are NEAs whose
Minimum Orbit Intersection Distance (MOID) with the Earth is 0.05
AU and whose diameter is greater than 140 meters (minimum diameter
needed in order to penetrate the Earth’s atmosphere).
PN formalism: Post-Newtonian formalism - In situations where

gravity is weak and the objects are moving slowly compared to the
speed of light, the theories of general relativity and that of Newton’s
gravity lead to very similar descriptions of the motion of bodies. This
can be seen as taking the Newtonian description and adding successive
correction terms that take into account the effects of general relativity.
The method for successively adding these correction terms is called the
PN formalism.
PPN formalism (Parametrized Post-Newtonian formalism):

Is a version of the PN formalism that explicitly details the parameters
in which a general theory of gravity can differ from Newtonian gravity.
Propagator: Simulates the trajectory of one or more objects over

time given these bodies’ initial positions and velocities.
State vector: Vector whose components are the positions and veloc-
ities of a given body. As the positions and velocities depend on time,
the state vector is time-dependent.
viii
Section 1
Introduction and Motivation
Figure 1.1: Orbits of over 1000 known Potentially Hazardous Asteroids (PHAs).
These are over 140 meters in diameter and will pass within 7.5 million kilometers
of Earth – about 20 times the distance to the Moon.
† Image taken from http://apod.nasa.gov/apod/ap130812.html
The first asteroid ever discovered was Ceres in 1801 [21]. Since then,
around 600’000 asteroids have been discovered in our solar system. Recently,
these objects have become a source of scientific research due to a few events
that occurred in the past years. In fact, in 1908, an asteroid entered into
the Earth’s atmosphere and exploded in the sky above Siberia. This has
been documented as the largest impact event on Earth. More recently, in
2013, the Chelyabinsk meteor entered Earth’s atmosphere and exploded in
air causing injuries to almost 1’500 people.
1
SECTION 1. INTRODUCTION AND MOTIVATION
Event type Diameter Fatalities Impact interval

High alt. beak-up < 50 m 0 annual
Regional > 140 m 50’000 5’000 years
Large sub-global > 300 m 500’000 25’000 years
High global > 5 km > 2 billion 6 million years
Extinction-class > 10 km 6 billion 100 million years
Table 1.1: Typical impact types, frequency and conesquences. [1]
Table 1.1 shows a summary of collision event types with their associated
object diameters, number of fatalities and typical impact intervals.
In the scientific community, it is widely believed that a large asteroid col-
lided with our planet approximately 65 million years ago, resulting in the
extinction of dinosaurs. It is not unreasonable to assume that if such an
event were to take place again, it could lead to the demise of the human
race. Hence, the question may not be if an event such as the one respon-
sible for the extinction of dinosaurs will happen, but rather when it will
happen.
Assuming then that we will one day be faced with an Earth-asteroid colli-
sion, it would be prudent to have a mitigation plan, enabling us to avoid
catastrophic consequences and perhaps even save our species.
In view of this, PhD candidate Sung Wook Paek from MIT’s AeroAstro
department decided to focus his doctoral research on asteroid mitigation
missions. More specifically, he will prove the utility of precursor missions
under high initial uncertainties concerning potentially hazardous asteroids.
In fact, there is great uncertainty surrounding certain asteroid parameters
that greatly impact our ability to predict their trajectories. These uncer-
tainties relate to mass, density, shape, etc. and make mitigation missions
tricky. Because it is difficult to reduce these uncertainties with remote ob-
servations from Earth, we can consider precursor missions whose goals are
to obtain valuable information about the potentially hazardous body that
could then be used for an effective and optimized deflection process. There-
fore, Sung Wook’s goal will first be to prove the utility of a precursor mission
and second, he will model and optimize the details of a two-stage mitigation
campaign consisting of a precursor mission and a mitigation mission.
Now, in order to simulate a mitigation mission, it is necessary to estimate fu-
ture trajectories of the potentially hazardous objects. Moreover, it is crucial
to know whether such a deflection mission is necessary. Therefore one must
be able to predict an asteroid’s passage through what is called a keyhole.
2
Figure 1.2: Differences between orbits of a typical near-Earth asteroid (blue) and
a potentially hazardous asteroid, or PHA (orange). PHAs are a subset of the near-
Earth asteroids (NEAs). They have the closest orbits to Earth’s orbit, coming
within about 8 million kilometers, and they are large enough to survive passage
through Earth’s atmosphere and cause damage on a regional, or greater, scale .
† Image taken from
http://www.nasa.gov/mission_pages/WISE/multimedia/gallery/neowise/pia15628.html
Simply put, these are small regions in space with the interesting property
that a body passing through one of these results in it colliding with Earth
(keyholes will be discussed further in Section 2.2.2). In turn, the study of
keyhole passages requires a precise trajectory propagator.
The goal of this thesis will therefore be to develop a high fidelity and fast-
running orbit propagator that can be used both for determining whether or
not an asteroid will pass through a keyhole and for simulating a deflection
mission. Of course, there are propagators that already exist, some of which
are available for public use. However, some of these public propagators do
not make their source code available, making the available options limited
for use in the context of deflection missions. We will explore some of the
ready-made propagators in this thesis.
One of the main motivations for this thesis is to obtain an orbit prop-
agator where the user can explicitly tradeoff accuracy versus computation
time. This tradeoff will be discussed in the case studies in Section 5. In fact,
depending on the work at hand, one might be more interested in obtaining a
very precise ephemeris with a long runtime whereas others might find more
use in a less precise ephemeris whose runtime is shorter. We want to pro-
3
vide a propagator that is easy to use and that is flexible in the equations,
number of bodies, relativity theory and output and, of course, that is open
source. It should not be required to have advanced knowledge in aerospace
engineering, physics or mathematics in order to use the tool. However, it
should be possible to modify it for those with the knowledge and interest.
The gap that is being addressed with this thesis is that of providing a simple
and flexible orbit propagator that allows the user to propagate and study
not only an asteroid’s trajectory but also all the other bodies that are con-
sidered in the physical model and that only requires initial conditions for
the bodies. Results for five asteroids are given in order to help the user in
chosing the numerical tolerance given desired runtimes and accuracy and
will procvide ball park values for the errors.
This thesis will focus in particular on the asteroid Apophis. Apophis was
discovered only ten years ago, in 2004, and was put in the spotlight when
initial observations and computations indicated a high probability (around
2.7%) that it would collide with Earth in 2029 [22]. We will also study four
other asteroids from NASA’s Sentry table [25] that lists the bodies with
potential future Earth impact events.
The remaining part of this thesis is structured as follows. Section 2

provides an overview of the physical and mathematical models that will be
used throughout this research. The equations of motion are described as
well as the variables that come into consideration. Section 3 describes some
aspects of geometric numerical integration that are relevant to the given
problem as well as the requirements for deflection missions. In Section 4, a
new propagator will be introduced along with its mathematical foundation.
Section 5 will be dedicated to testing several different propagators on a
benchmark problem, specifically that of the asteroid Apophis. Using the
results from Section 5, the best available propagator will be selected and
Section 6 will give a detailed analysis of two case studies using the chosen
propagator. Finally, Section 7 will contain conclusions to the problem in
general and specific to asteroid deflection missions. Future work will be
briefly described.
4
Section 2
Physical and Mathematical

Background
Figure 2.1: Illustration of the N-body problem.

† Image taken from
http://www.lactamme.polytechnique.fr/images/NCOR.U1.2048.D/display.html
This section will provide an overview of the physical and mathematical

formulation used to model the movement of bodies in our solar system.
When considering N celestial objects (the Sun, planets, moons, et.), this is
referred to as the N-body problem. Simpler versions of this are the 2-body
problem (e.g. the Earth-Sun system) and the 3-body problem (e.g. the
Earth-Moon-Sun system). By including more bodies into the model, i.e. by
5
SECTION 2. PHYSICAL AND MATHEMATICAL BACKGROUND
increasing N, we obtain more precise trajectories.
2.1 The N-body problem

The point mass N-body problem may be stated:
Given the positions and velocities of bodies with known mass at some initial
time t0 , find the positions of the bodies at any other time t > t0 .
The N-body problem is chaotic [2]. Simply put, even if the present
determines the future, the approximate present does not approximately de-
termine the future. In terms of the N-body problem, this means that small
perturbations in mass, trajectories or other physical properties will lead to
enormous and unpredictable changes in the system. This is a purely physical
phenomenon. However, chaos has a numerical counterpart called instabil-
ity. Similarly to the physical meaning, numerical instability means that by
modifying the input data slightly, we will obtain very different numerical re-
sults. While instability is a negative feature of any problem, it does remove
the burden of attempting to identify an integrator that preserves inherent
stability. In fact, when dealing with physically non chaotic (or stable) sys-
tems, we try to choose a numerical scheme that will preserve the stability.
However, if the physical system is chaotic to start with, there is no point in
using a numerically stable integration scheme.
Even though the N-body problem is unstable, a convergent expansion

does exist for a general N as proved by Wang in 1991 [3]. This result is far
from trivial. In fact, in 1912, i.e. almost 80 years before Wang generalized
the result, Sundman developed a method that produced a global analytical
solution of the 3 body problem. Unfortunately, these expansions cannot, as
a practical matter, be applied to real problems as they require millions of
terms to be computed even for short times. Therefore solving the N-body
problem requires numerical integration of the equations of motion. More-
over, in the abstract of his paper "The existence of global solution of the
N-body problem" [3], Wang writes:
"The convergence of my power series is admittedly unsatisfactory, and so

the present result is of limited value for practical calculation."
It is noteworthy that in his proof, Wang considers bodies as point masses

(as do we) but does not take into account the relativistic terms that will be
discussed in 2.1.1. Furthermore, there is no mention of the uniqueness of the
solution, i.e. the problem is not well posed according to the mathematical
definition. Considering that it took almost 80 years and a new mathematical
6
Figure 2.2: Illustration of the ecliptic plane. This plane is defined as the imaginary
plane containing the Earth’s orbit around the Sun. The planetary bodies of our
solar system all nearly lie in this plane. Originally, the ecliptic plane was defined
using the path of the Sun around the Earth.
† Image taken from http://www.herongyang.com/astrology_horoscope/ecliptic_plane
transformation in order to prove the existence of a solution to the N-body

problem with no relativistic terms, I will not attempt to extend this proof
to the equations of motion that include these extra terms. In summary, we
do not actually know whether or not there is a solution to the problem that
we will be solving and if it does exist we do not know whether or not it is
unique.
We will therefore be solving the N-body problem numerically in order to
determine a given asteroid’s trajectory in space.
2.1.1 Equations of motion

One of the physical models currently used in trajectory propagation includes
the mutual Newtonian gravitational accelerations and their relativistic cor-
rections, which is a modified form of the Einstein-Infeld-Hoffmann equation.
7
For each body i, the point mass acceleration is given by:

X µj (rj − ri ) 2(β + γ) X µk 2β − 1 X µk
r̈i = 3 1− −
j6=i
rij c2 r
k6=i ik
c2 k6=j rjk
2 2
vi vj 2(1 + γ)

+γ + (1 + γ) − ṙi · ṙj
c c c2
" #2
3 (ri − rj ) · ṙj 1

− 2 3 + 2 (rj − ri ) · r̈j (2.1)
2c rij 2c
1 X µj
+ 2 3 {[ri − rj ] · [(2 + 2γ)ṙi − (1 + 2γ)ṙj ]} (ṙi − ṙj )
c j6=i rij
(3 + 4γ) X µj r̈j
+ ,
2c2 j6=i rij
where ri , ṙi and r̈i are the solar-system-barycentric position, velocity and
acceleration vectors of body i, i.e. the positions and velocities of body i using
the center of mass of the solar system as the center of the reference frame.
This means that we are considering an inertial reference frame centered at
the center of mass of the solar system. The constant c is the speed of light,
µi = Gmi , where G is the gravitational constant and mi is the mass of body
i. Furthermore, rij = |rj − ri | is the distance between bodies i and j and
vi = |ṙi |. Finally, β and γ are PPN parameters measuring the nonlinearity
in superposition of gravity and the space curvature produced by unit rest
mass, respectively. In general relativity, β = γ = 1.
For the rest of this paper, the word relativistic will refer to considering
the relativistic terms in the equations of motion. In contrast, the word
Newtonian will be used when these terms are not included. In equation
(2.1), the relativistic terms are written in gray. They are all the terms in c12 .
Because acceleration terms appear both on the right and left hand side of
equation 2.1, this problem is implicit. In fact, if
 
r1
 . 
 .. 
 
 r 
 N 
y= ,
 ṙ1 
 .. 
 
 . 
ṙN
8
is the state vector of the N bodies considered, then

 
ṙ1
 . 
 .. 
 
 ṙ 
y0 =  N 
 
 r̈1 
 .. 
 
 . 
r̈N
is the vector containing velocities and accelerations and the N-body problem
can be formulated as
f (t, y, y 0 ) = 0 (2.2)
where we omit the time-dependency to simplify notation.
However, it is possible to exploit the structure of the equations of motion
in a more elegant manner by putting all the acceleration terms to the left
hand side of the equality sign. We then obtain:
M (y)y 0 = f (t, y) (2.3)
where M is called the mass matrix and depends on the state vector y.
In contrast to the implicit problem (2.2), an explicit problem could be writ-
ten as
y 0 = f (t, y),
which is equivalent to (2.3) when we replace f (t, y) with M −1 f (t, y). Of
course, this requires an invertible mass matrix M . If the matrix is singular,
a different approach must be used. We will not be discussing this case fur-
ther in the context of this thesis.
In the context of numerical integration, explicit problems are preferred.

In fact, when dealing with an implicit problem, one must solve an implicit
sub-problem at each time step, which is computationally expensive. For
this reason and for comparison purposes, we introduce slightly modified
equations:
X µj (rj − ri ) 2(β + γ) X µk 2β − 1 X µk
r̈i = 3 1− −
j6=i
rij c2 r
k6=i ik
c2 k6=j rjk
2 2
vi vj 2(1 + γ)

+γ + (1 + γ) − ṙi · ṙj
c c c2
" #2
3 (ri − rj ) · ṙj
− 2 3 (2.4)
2c rij
1 X µj
+ 2 3 {[ri − rj ] · [(2 + 2γ)ṙi − (1 + 2γ)ṙj ]} (ṙi − ṙj ),
c j6=i rij
9
where we simply omit the two terms in (2.1) where the acceleration r̈j ap-
pears. Again, the terms in gray refer to the relativistic corrections of general
relativity theory. This way, the N-body problem can be re-written as:
y 0 = f (t, y).
For the sake of comparison, we also introduce the Newtonian equations of

motion:
X µj (rj − ri )
r̈i = 3 . (2.5)
j6=i
rij
Now, in order to obtain the positions of each body at a given time t > t0 ,
we need to integrate the equations of motion.
In what follows, the N bodies considered will be the Sun, the 8 planets,
Pluto, the Earth’s Moon along with the asteroids Ceres, Pallas, Vesta and
the asteroid to be studied. Ceres, Pallas and Vesta are the most massive
asteroids in our solar system, accounting for about 46% of the total mass of
the asteroid belt. They are often referred to as the "Big 3". The asteroid
belt is illustrated in Figure 2.3.
Figure 2.3: The asteroid belt lies in the region between Mars and Jupiter. The
Trojan asteroids lie in Jupiter’s orbit, in two distinct regions in front of and behind
the planet.
† Image taken from https://solarsystem.nasa.gov/multimedia/display.cfm?IM_ID=850
Additional forces may be included in the physical model in order to

obtain more accurate numerical results. Some of these forces are:
10
• Solar radiation pressure or SRP:

The Sun exerts electromagnetic radiation on the bodies in the solar
system. This radiation is then absorbed and/or reflected by the body
receiving the radiation which impacts this body’s motion. Solar radi-
ation pressure affects mostly small bodies.
• Drag:
Several types of drag can influence a celestial body’s trajectory. For
example, magnetic drag and air drag.
• Solar oblateness:
The Sun is not exactly spherical. Hence its gravitational field is not
entirely symmetric as it is modeled using point masses. In order to
correct for this, one needs to take into account the non-sphericality of
the Sun. This is also valid for every one of the celestial bodies, none
of which are perfectly spherical.
• Perturbations from the 300 most massive asteroids in the asteroid belt:
Just as the planets, Sun and moons attract bodies in space, massive
asteroids have a gravitational pull and therefore influence the trajec-
tories of bodies that pass at a certain distance of them.
• Perturbations from the less massive asteroids in the belt, modeled by

a homogeneous ring:
The same argument as the previous one holds for the less massive
asteroids. However, as there is a very large number of these objects
that are relatively evenly spaced in the asteroid ring, we can imagine
a homogeneous asteroid ring that would account for the gravitational
pull of all these smaller asteroids combined.
• Perturbations from the 21 largest Trans-Neptunian Objects (or TNOs):

A TNO is any minor planet in the solar system that orbits the Sun
with a greater semi-major axis than Neptune. As for the most mas-
sive asteroids, TNOs can be taken into account for a more precise
propagator.
• Perturbations from other, less massive, TNOs, modeled by a homoge-

neous ring:
Similar to the less massive asteroids in the asteroid belt.
Each one of these forces on its own accounts for a very weak gravitational
effect. However, taken together they would contribute to the accuracy of
the numerical solution. However, for the sake of asteroid deflection missions,
these terms are not necessary and therefore are not included in the model.
They are mainly of interest when building an extremely precise planetary
ephemeris propagator.
11
2.1.2 Hamiltonian formulations

When dealing with geometric numerical integration, one usually tries to
write the problem at hand using the Hamiltonian formulation. This is a
mathematical formalism describing the equations of motion of a dynamic
system. Usually, the Hamiltonian H represents the total energy of the sys-
tem at hand, i.e.
H = T + V,
where T and V represent the kinetic and potentiel energy, respectively. More
generally, suppose that H(p, q) : R2d → R is a smooth function of its argu-
ments p, q ∈ Rd . Then the dynamical system t 7→ (p(t), q(t)) defined by the
equations 
∂H(p(t),q(t))
 ṗk (t) = −

 ∂qk
 q̇k (t) = ∂H(p(t),q(t))



∂p k
for k = 1, . . . , d, is called a Hamiltonian system. Usually, the variables q

and p refer to the generalized positions and momenta (or velocities) of the
system, respectively. Therefore, we can write the traditional state vector as
y = (p, q)T . With this, we can write the Hamiltonian equivalently as
ẏ = J −1 ∇H(y), with J = 0 Id
−Id 0 .
Note that it is not always possible to find a Hamiltonian formulation for a

given problem and as it turns out, it seems extremely complicated to find
one for the problem that we will be considering because of the relativistic
terms included in the equations of motion. Furthermore, seeing as we want
to be able to modify the equations of motion by adding perturbing forces to
model deflection processes, it is better to not assume a Hamiltonian system,
allowing for more flexibility in the perturbing forces.
The reason why this Hamiltonian is so important for geometric numerical
integration and therefore the reason why mathematicians and engineers use
it whenever it is possible is because it has been found that by using this
formalism, some very nice results are true and can be used. These results
are linked to the notion of symplectic solvers.
A linear map A : R2d → R2d is said to be symplectic if
AT JA = J,
with J defined above. Similarly, a differentiable map g : U → R2d is said to

be symplectic if
T
g 0 (p, q) Jg 0 (p, q) = J ∀(p, q) ∈ U.
Symplectic maps have the property of preserving volumes. Now, if Φh (y0 ) is

the numerical flow of a one-step method (see 2.2 for more about numerical
12
methods and the definition of one-step methods), i.e. y1 = Φh (y0 ), then the
method is said to be symplectic if
∂ΦTh ∂Φh
(y0 )J (y0 ) = J ∀y0 .
∂y0 ∂y0
It can be shown that symplectic solvers nearly preserve the Hamiltonian of a
system. Hence, when the Hamiltonian represents the total energy and when
using a symplectic solver, the numerical total energy is preserved just as it
is in the physical one (when conservation of energy applies). This is usually
a very attractive property.
We will not be considering symplectic solvers. In fact, the perturb-

ing forces that will be included in order to model the asteroid deflections
break the conservation of energy and in fact traditional Lagrangian and
Hamiltonian mechanics cannot be used with non-conservative systems. Note
that even if a Hamiltonian were available for the N-body problem, it would
almost certainly not be separable, i.e. it would not be possible to write
H = T (v) + U (r) where the kinetic energy T only depends on the velocities
and the potential energy U only depends on the positions. As a consequence,
it would not be possible to use an explicit method of integration, unless a
splitting is used. In a splitting, we write r̈ = (A + B)r where A is the
Hamiltonian part and B is whatever cannot be written in this formulation
(this can be seen as a perturbation from a Hamiltonian system). Then we
can use a symplectic method on A.
2.1.3 Derivation of the equations of motion

The equations of motion given in equation (2.1) are to be understood as a
perturbed version of the Newtonian equations of motion shown in equation
(2.5). The perturbations, which will often be referred to as the relativistic
terms, account for effects from general relativity. The Newtonian, or unper-
turbed equations come from Newton’s second law of motion F = ma and
his law of universal gravitation F = GM r2
m
. Below, we will give a short de-
scription of how the relativistic terms came to be included in the equations
of motion.
In 1957, Infeld suggested deriving the post-Newtonian equations for the

N-body problem from the variational principle of the Theory of General
Relativity, or GRT field equations, using the linearized mass tensor. We
will not give extensive details here. Instead, we refer the reader to [4] and
[17] for more details.
As a starting point, we need the GRT field equations, also known as the
Einstein Field Equations, or EFE. In fact, the basic idea of GRT is that
the properties of spacetime are determined by motion and distribution of
13
masses and, conversely, motion and distribution of masses are determined

by the field metric. The EFE equation is:
1 8πG
Rµν − Rg µν + g µν Λ = 4 T µν , (2.6)
2 c
where Rµν is the Ricci curvature tensor, R is the scalar curvature, g µν the
metric tensor, Λ the cosmological constant (whose value is the energy den-
sity of the vacuum of space), c the speed of light, G Newton’s gravitational
constant and T µν the stress-energy-momentum tensor. The indices µ and ν
are the spacetime dimensions. They can each take 4 values: 0,1,2 or 3 where
0 represents time, 1 the x variable, 2 the y variable and 3 the z variable.
Note that the term g µν Λ is very small so it is almost always left out (it is
only useful when dealing with large cosmological scales).
These equations describe the fundamental interaction of gravitation as a
result of spacetime being curved by matter and energy. They are used to
determine the geometry of spacetime resulting from the presence of mass-
energy and linear momentum. The solutions of EFE, are the components of
the metric tensor g µν . Exact solutions for the EFE can only be determined
under simplifying assumptions such as symmetry.
The left hand side of the EFE refers to the curvature of spacetime while the
right hand side has to do with mass and energy. This is the basic idea of
GRT: the properties of spacetime, i.e. the spacetime metric are governed
by the motion and distribution of masses and, conversely, the motion and
distribution of masses are determined by the field metric. Hence the EFE
can be understood as follows:
Mass tells spacetime how to curve and curved spacetime tells mass how to
move.
Now, let
√
Z
Sg = −gRdΩ,
where dΩ = cdtd3 x. We want to compute the variation of this integral.

Using the variation of the Ricci tensor, the definition of the determinant g
and the definition of the Christoffel symbols, it can be shown that
1 √
Z
δSg = − Rµν − Rg µν −gδg µν dΩ. (2.7)
2
Notice that the inside of the parenthesis coincides with the left-hand side of
the EFE equation (2.6). Now consider the action integral
√
Z
Sm = ρ ∗ (c2 + Π) −gdΩ,
14
with ρ∗ being an invariant density satisfying the equation of continuity and

Π being the potential compressional energy. Using the variation of ρ∗, we
find
1 2 √
Z
δSm = c T µν −gδg µν dΩ (2.8)
2
Using (2.6) along with (2.7) and (2.8), we see that the EFE follow from the
variational equation
δ(2c−2 κSm − Sg ) = 0. (2.9)
c4
Multiplying (2.9) by 16πG , we can rewrite it as
Z " #
c4 √ ds
δ −gJ + c2 (1 + c−2 Π)ρ 0 dΩ = 0
16πG dx
Using the fact that dΩ = cdtd3 x is the elementary 4-volume, Infeld and
Plebansky [4] show that this principle can be re-written
Z
δ Ldt = 0
with Z " #
c4 √ ds
L=− −gJ + c2 (1 + c−2 Π)ρ 0 d3 x
16πG dx
to treat L as the Lagrangian of the post-Newtonian N-body problem. The
integral is to be understood as the sum of all the integrals taken over the
volumes of the bodies. After a few pages of math, one comes to the following
result, also known as the Einstein-Infeld-Hoffman, or EIH equation:
X 1 1 X mi mj 1 1 1 X mi mj

L= mi vi2 + G + 2 mi (vi2 )2 + G
i
2 2 j6=i rij c 8 4 j6=i rij
1 1 X mi mj (mi + mj )

× 3vi2 + 3vj2 − 7vi · vj − (vi rij )(vj rij ) 2 − G2 2
rij 4 j6=i
rij
1 X X 1 1 1

− G2 mi mj mk + + .
6 j6=i k6=i,j
rij rik rji rjk rki rkj
In 1971, Frank B. Estabrook from NASA’s Jet Propulsion Laboratory (JPL)

used this Lagrangian in order to derive the conservative PPN N-body La-
15
grangian [9]:
1X 1 X 1 + 2γ X X µi µj 2
L= µi vi2 + 2 µi (vi2 )2 + (vi + vj2 )
2 i 8c i 4c2 i j6=i
rij
3 + 4γ X X µi µj
− ṙi · ṙj (2.10)
4c2 i j6=i
rij
1 X X µi µ j 1 X X µ i µj
− 2 3 {(rj − ri ) · ṙj } {(rj − ri ) · ṙi } +
4c i j6=i rij 2 i j6=i rij
2β − 1 X X µi µj (µi + µj ) 2β − 1 X X X µi µj µk
− 2 − ,
4c2 i j6=i
rij 2c2 r r
i j6=i k6=j ij jk
This Lagrangian will in turn yield the equations of motion given in (2.1).
2.2 Numerical Integration

In numerical integration of dynamical systems, we look for the solution y(t)
to problems of the form
(
y 0 (t) = f (t, y(t))
y(0) = y0
with f : R × Rn → Rn and y0 ∈ Rn for some n ≥ 1. These are also known

as Initial Value Problems, or IVPs because of the necessity of providing the
initial conditions, or initial values, y0 .
Given this problem and a time step h ≥ 0, the fundamental theorem of
calculus reads:
Z t0 +h
y(t0 + h) = y(t0 ) + f (τ, y(τ ))dτ.
t0
It is the integral that must be approximated numerically in order to obtain

the solution at the end time t0 + h.
The solution at time t given the initial conditions y0 is denoted
y(t, 0, y0 ) = φt (y0 ),
and is called the exact flow of the problem. Similarly, the numerical flow is
Φh (y0 ) with h the time step.
2.2.1 Geometric numerical integration

Considering the nature of the problem and its solution, geometric integra-
tion is most appropriate. In fact, geometric integration deals with integra-
tors that preserve some geometric properties of the exact solution. In fact,
16
physical problems often have properties of conservation or symmetry. In

some cases, these properties are fundamental and therefore it is desirable to
obtain a numerical method that preserves them. However, not all problems
have these inherent properties. It is therefore important to understand the
physics of the system modeled in order to know whether or not we want
the numerical scheme to preserve certain properties. Often, we look at re-
versibility, symplecticity and conservation of first integrals. These can be
mathematically formulated as follows:
• Reversibility:
Φ−h (Φh (y0 )) = y0 ,
i.e. we should be able to go forward in time with a step size h and
then back to where we started (with a step size −h) and arrive exactly
where we started.
• Symplecticity:
∂ΦTh ∂Φh
(y0 )J (y0 ) = J ∀y0 .
∂y0 ∂y0
It is difficult to give a geometric or physical explanation of this prop-
erty without introducing manifolds and symplectic geometry. For this
reason, we will simply say that it is a mathematical property of nu-
merical integrators that can be used in simple classical mechanics.
• Conservation of first integrals:
I(y1 ) = I(y0 )
where y1 = Φh (y0 ) and I is a first integral of the system. First integrals

d
are characterized by the simple equation dt I(y(t)) = 0, meaning that
I is constant along the trajectory y(t).
It turns out that the preservation of these properties makes for a favorable
long-term behavior.
However, implementing solvers that preserve symplecticity, for example, re-
duces the flexibility in equations of motion, as discussed in Section 2.1.2.
2.2.2 Accuracy requirements

In the context of this thesis, it is important to introduce the concept of
keyholes. Keyholes are defined as specific regions in space with a very par-
ticular property. In fact, if an object were to pass through a keyhole, it
would eventually collide with Earth. The smallest keyholes are less than
1 kilometer in diameter and the largest ones are about 1000 kilometers in
diameter. This is the reason behind the need for the high accuracy of the
propagator stated, as stated in Section 1. If we want to know whether or
17
not an asteroid will be passing through a keyhole, we need it’s position to be

known with an error less than 1000 kilometers (assuming a large keyhole),
which is quite a small distance in astronomical terms.
Furthermore, since we would like to be able to study slow push methods
of deflection, the timespan of the numerical integration needs to be long
enough (at least 10 years which is about the time necessary to complete
slow push methods).
We remark that the only accuracy requirement that remains when all exter-
nal sources of uncertainty are removed is that of the size of a keyhole. This
means that, assuming no propagation or measurement errors, it should be
possible to avoid an asteroid passing through a keyhole by deflecting it as
much as the diameter of the keyhole. For example, Apophis’ 2036 keyhole
is approximately 700 meters in size and therefore a 1 kilometer deflection
would move Apophis (whose diameter is about 325 meters) completely out
of this keyhole.
18
Section 3
PAT2: Propagator for

Asteroid Trajectories Tool
Figure 3.1: An intercept satellite racing toward an asteroid. Studies by Iowa

State University’s Bong Wie indicate such a vehicle could blast apart asteroids
that threaten Earth.
† Image taken from http://www.news.iastate.edu/news/2013/03/06/asteroiddeflection
As stated previously, we want to explore different options in terms of

possible propagators.
For this reason, we build a new propagator, PAT2 (pronounced pat squared)
using MATLAB. Similarly to other propagators, PAT2 requires the user to
input a time span as well as initial positions and velocities of the bodies
19
SECTION 3. PAT2 : PROPAGATOR FOR ASTEROID
TRAJECTORIES TOOL
considered and outputs an ephemeris giving the positions (and velocities)

of all these objects at any time in the specified timespan. The flow chart
illustrating the inner working of PAT2 is shown in Figure 3.2, where the
arrows show the direction of flow between the pink ellipses representing the
input data, the green rectangles the processes and the purple ellipses the
output data. The timespan and initial conditions are necessary inputs. All
other inputs have default values but can be easily modified, if necessary.
3.1 Mathematical and Physical Model

The mathematical and physical foundation for PAT2 is that described in
Section 2.1.1, with the equations of motion including all relativistic terms.
The bodies included in the model are the Sun, the 8 planets, Pluto, the
Moon, the Big 3 (Ceres, Pallas and Vesta) as well as the asteroid whose
trajectory we would like to study. Note that the initial conditions can be
taken to be the positions and velocities of the barycenters of each body
rather than the body centers (except for planet Earth because its moon is
explicitly included in the model). For bodies with no moons, this will not
change anything. However, for planets such as Jupiter, who has 50 moons or
Saturn who has 53 moons, there will be a difference between the barycenter
of the planet-moon system and the center of the planet. Considering these
barycenters instead of the body centers can somehow compensate for the
fact that the moons are not included in the physical model. Note that if we
consider the barycenter of a planet-moon system, we must adapt the mass,
i.e. if we want to model the Jupiter-moons system, we must use the total
mass of the system which is the sum of the planet and its moons.
In contrast to most propagators out there, PAT2 only uses data for initial
conditions and treats all bodies the same way. In fact, most of the other
asteroid propagators available only propagate the trajectory of one single
body, the asteroid. The positions and velocities of all other bodies in the
model are queried from a previously generated ephemeris. The reason for
going with the first philosophy (of propagating all bodies rather than just the
asteroid) is simple: in the latter situation, each time step requires a query
from an ephemeris file. This requires access to such a file. Moreover, the
code must be modifiable in order to include perturbing forces, e.g. kinetic
impactors or slow push methods and therefore the equations of motion must
be available and they must be easy to modify.
Of course, choosing to propagate 15 bodies instead of just 1 will make for a
considerably longer runtime. However, as long as runtime stays reasonably
short when propagating the 15 bodies, this method is acceptable.
Note that is it also possible to modify the PPN parameters β and γ . As a
reminder, β and γ measure the nonlinearity in superposition of gravity and
20
TRAJECTORIES TOOL
β and γ Init. state vector Timespan
Initialize
model
Dist. unit Rescale Time unit
Numerical
integration
Output
Asteroid Total
State Numerical Numerical

vectors error energy
Planets Maximum Time t Ratio
Figure 3.2: Flow chart describing the inner workings of PAT2 . The pink ellipses
represent the input data, the green rectangles the processes and the purple ellipses
the output data. The arrows show the direction of flow. The timespan and initial
conditions are necessary inputs. All other inputs have default values but can be
easily modified.
21
TRAJECTORIES TOOL
Maximum Minimum
Positions 4.8 · 109 8.1 · 103
Velocities 3.7 · 101 2.4 · 10−4
Table 3.1: Maximal and minimal (absolute) values for the positions, velocities and
accelerations using units of seconds for time and kilometers for distances. Using
these values will lead to significant roundoff errors.
Maximum Minimum
Positions 4.8 · 104 8.1 · 10−2
Velocities 1.2 · 104 7.6 · 10−2
Table 3.2: Maximal and minimal (absolute) values for the positions, velocities and
accelerations using units of sidereal years for time and hundred thousand kilometers
(105 km) for distances. Using these values significantly reduces roundoff errors and
enables meaningful values of tolerances.
the space curvature produced by unit rest mass, respectively. In general

relativity, β = γ = 1. Modifying these values allows exploring different
relativity theories.
3.2 Rescaling
Once the user has supplied the initial conditions, PAT2 rescales the problem
so as to avoid the unwanted effects of roundoff error and to ensure that all
variables are at the same scale. In fact, the N-body problem considered here
with its usual units of time (seconds) and distance (kilometers) is prone to
a lot of numerical noise. For example, when computing the first term of
(2.1) which is proportional to r13 with r = O(108 ) we come below machine
precision and therefore introduce unnecessary numerical errors. Of course,
this will depend on the machine on which the simulation is run and on the
language used. With a 64 bit machine (MacBook Pro, Intel Core i5) and a 64
bit double precision language (MATLAB R2014a), the machine precision is
≈ 2.2·10−16 . As a reminder, machine precision is the smallest value , such
that 1 + and 1 have different floating point representations. Computing
with values smaller than causes a lot of unwanted and avoidable errors.
Note that a way of overcoming this issue could be to use another language,
such as FORTRAN that uses 128 bit representation. However, we will chose
to stay with MATLAB and rescale the problem.
The units chosen for the rescaling are sidereal years for time and hundreds
of thousand kilometers (105 ) for distances. This way, we avoid unnecessary
roundoff errors. Furthermore, by doing this we obtain positions, velocities
22
TRAJECTORIES TOOL
(and accelerations) that are of the same scale, which is a desirable property
when assigning tolerances to the numerical method. In fact, the variables
in this problem are of significantly different scales, as shown in Tables 3.1
and 3.2.
3.3 Numerical Integration Scheme

Now, with the problem rescaled, PAT2 calls the MATLAB function ode113.
This is a built-in function aimed at solving non-stiff differential equations.
It is a variable step size and variable order method, with orders ranging
from 1 to 13. ode113 was designed for problems with stringent error tol-
erances and for solving computationally intensive problems. It is based on
a variable order Adams-Bashforth-Moulton (or ABM) solver. This solver is
a predictor-corrector method based with PECE (Predict-Evaluate-Correct-
Evaluate) implementation. ABM methods are, as are all numerical inte-
gration methods, based on the fundamental theorem of calculus. Given a
problem
y 0 (t) = f (t, y(t)),
and points t0 , t1 , . . . , tk , tk+1 , we have
Z tk+1
y(tk+1 ) = y(tk ) + f (τ, y(τ ))dτ, k ∈ N.
tk
ABM methods use a number of previous points (tk−s , fk−s ),...,(tk−1 , fk−1 ),
(tk , fk ) to construct the Lagrange polynomial approximation p(t, y(t)) of
f (t, y(t)) passing through these points. Then
Z tk+1
y(tk+1 ) ≈ y(tk ) + p(τ, y(τ ))dτ, k ∈ N.
tk
can easily be integrated over the interval [tk , tk+1 ] as p it is a polynomial.

This method, known as the Adams-Bashfourth method (AB), is explicit
since all the pairs (tk−s , fk−s ),...,(tk , fk ) are known at time tk .
(AB) produces the predictor formula
Z tk+1
pk+1 = yk + p(τ, y(τ ))dτ
tk
that can be used for a rough estimate (tk+1 , pk+1 ) of (tk+1 , y(tk+1 )).
We can now do the same thing and construct a Lagrange polynomial approx-
imation of f , this time using the unknown point (tk+1 , fk+1 ). This yields
an implicit method known as the Adams-Moulton method (AM) which pro-
duces the corrector formula. Since it is an implicit method, it requires several
iterations. The value of the corrector pk+1 is used as the initial guess for
yk+1 . The corrector formula should give a more accurate approximation to
23
TRAJECTORIES TOOL
Solver
ode45 ode113
10−4 12’692 6’779
Tolerance
10−5 16’226 8’667
10−6 21’224 10’830
10−7 34’544 13’179
Table 3.3: Number of function evaluations needed to obtain a certain tolerance

level. Here, we consider the barycenter model used to propagate the positions of
the asteroid Apophis.
y(tk+1 ) than the predictor would because it is of higher order. Therefore,

by comparing the two, predictor and corrector, we can determine if the step
size should be increased or decreased.
This is a multistep method because it requires the solutions at several

preceding time points in order to compute the next solution. In contrast, a
one-step method only needs the previous solution in order to compute the
next. This is the case for Runge-Kutta methods, for example. Because ABM
is a multistep method, it requires the knowledge of several steps before it
can begin. These first steps typically are taken using one-step methods such
as Runge-Kutta.
The method presented here is not symplectic. The reason for not chos-
ing such a solver has been explained previously. Here, we remind the main
reason of this choice. As a reminder, the ultimate goal will be to have a
propagator that can be used in deflection missions and therefore some per-
turbing forces will be included in the equations of motions. These perturbing
forces will only be applied during a certain timespan, shorter than the in-
tegration timespan. The entire concept of deflection missions is based on
modifying the energy of a potentially hazardous asteroid in order to change
its trajectory, and hence we are not looking to conserve the total energy of
the system.
In order to evaluate ode113’s quality, we compare its results to those

obtained using ode45 (which is based on a Runge-Kutta method). The idea
is that the error is controlled in the numerical schemes through the choice of
tolerances. Therefore, the results obtained with ode113 or ode45 (or another
method) for a given tolerance should be almost the same. As a reminder, a
tolerance of 10−3 corresponds to 0.1% error. Similarly, a tolerance of 10−6
24
TRAJECTORIES TOOL
Error [km] Marg. improv. [km] Runtime [min]

10−3 ? ? ?
10−4 310.0709 - 8.5
Tolerance
10−5 273.1395 36.9314 10.5

10−6 272.4092 0.7303 11.3
10−7 272.4018 0.0074 15.5
10−8 272.4020 −2 · 10−4 20.5
Table 3.4: Numerical error (in absolute value, compared to HORIZONS) in the last
time step. Here, we consider the asteroid Apophis during the timespan 2014-2024.
Each error (computed against HORIZONS) is then compared to that obtained
with the previously assigned tolerance (one order of magnitude less), yielding the
marginal improvement of a smaller tolerance. The star symbol ? indicates that the
solver crashes given the specified tolerance.
corresponds to an error of 0.0001%. Tolerances have no unit.

For stable systems, if we control the error in each time step, we effectively
control the global error. In fact, for stable systems, the global error is
proportional to the per-unit-step error and therefore we can control the
global error by controlling the per-unit-step error. This however is not true
with unstable systems. In fact, when simulating an analytically unstable
system, it can no longer be assumed that the global integration error is
proportional to the per-unit-step integration error, since integration error
can accumulate excessively across multiple steps.
The difference in performance between two solvers will therefore be in how
efficient it is, i.e. how many times it needs to evaluate the equations of
motion. This in turn reflects how much runtime will be needed. Table
3.3 shows the total number of function evaluations for ode113 and ode45
given certain tolerances. We first notice that for both ode45 and ode113,
a smaller tolerance requires a higher number of function evaluations, as
expected. Moreover, the number of function evaluations, given a certain
tolerance, is always higher for ode45. This means that ode113 is more
efficient in obtaining similar results. Another way to look at it, is that
for a given number of function evaluations, ode113 will produce a solution
corresponding to one obtained with ode45 with a tolerance of at least two
to three orders of magnitude bigger.
Furthermore, we want to check that by using a smaller tolerance, we

obtain better numerical results. Table 3.4 shows the numerical errors ob-
tained in the last time step using ode113 with different tolerance levels.
25
TRAJECTORIES TOOL
The star symbol (?) that appreas in the first line of the table indicates that
the numerical scheme crashes given the tolerance of 10−3 . What happens is
that given a large tolerance, the solver produces numerical errors that are
relatively large. As the problem is unstable, these errors increase drastically
and therefore it is no longer possible for the solver to achieve desired accu-
racy without going below the allowed minimum step size. Therefore and the
scheme stops. Besides this phenomenon, we observe that when going from
a tolerance of 10−4 to 10−5 , i.e. increasing the accuracy, there is almost
a 40 kilometer improvement on the error in the last step. We notice that
this improvement dampens as tolerance is decreased more, suggesting that
the problem with this formulation and solver has a lower bound on error.
Note that there are different error metrics that can be considered. Here,
we chose to study the error in a given point (the last point). However, we
might be interested in comparing the "worst" or highest error that occurs
during the 10 year timespan. For the study of deflection, this metric makes
sense because of the ultimate goal to keep the error under a 1000 kilometer
threshold. Another error metric could be the integral of the error over the
entire integration. Then, by dividing by the timespan, we obtain an average
of the error. We will not explore this metric further.
3.4 Output
Once the integration is complete, PAT2 brings the results back to their
initial scales of seconds for time and kilometers for distances. If a timespan
is specified as [t0 , t1 , . . . , tn ], then the output of the numerical integration
will be the state vector evaluated in these n points. Otherwise, if only an
initial and final time t0 and tf are specified, then the solver will return a
structure which can be evaluated at any time between t0 and tf . In Appendix
A, we will briefly discuss polynomial interpolation that is used in order to
obtain solutions at any point in the timespan and the errors that occur as
a consequence.
Simple flags allow the user to make plots for the outputted trajectories, the
numerical errors and the energy.
26
Section 4
Benchmark comparisons
Figure 4.1: Artist’s conception of a catastrophic asteroid impact with the early
Earth.
† Image taken from http://solarsystem.nasa.gov/news/display.cfm?News_ID=23777
In this section, we will explore different possibilities regarding available

asteroid propagators and compare the results that they produce in order to
choose the solution that is most appropriate for our study.
4.1 Possible Candidates

As stated previously, there exist some ready to use propagators. Some ex-
amples are JPL’s HORIZONS, NASA’s GMAT (General Mission Analysis
27
SECTION 4. BENCHMARK COMPARISONS
Initialize
Init. state vector Timespan
model
Numerical
integration
ephemeris
Planetary
Output
Asteroid
state
vector
Figure 4.2: Flow chart describing the inner workings of HORIZONS (and many
other available propagators). The pink ellipses represent the input data, the green
rectangles the processes, the orange cylinder the required database and the purple
ellipses the output data. The arrows show the direction of flow
Tool) or AGI’s STK (Systems Tool Kit).

HORIZONS is an online solar system ephemeris computation service that
provides highly accurate ephemerides for solar system objects such as plan-
ets, asteroids, comets and satellites. Uncertainties in major planet ephemerides
in HORIZONS range from 10 centimeters to over 100 kilometers.The avail-
able timespan in HORIZONS differs from body to body, but in general
is quite large. Using its web-interface, one can specify the body whose
ephemeris is of interest as well as the coordinate origin, the time span of
the ephemeris and the units to be used. Unfortunately, the source code for
HORIZONS is proprietary and therefore it is impossible to see or modify
it. Hence it is not possible to model defleciton mission with this tool. The
flow chart illustrating the inner workings of HORIZONS is shown in Figure
4.2, where the arrows show the direction of flow between the pink ellipses
represent the input data, the green rectangles are the processes, the orange
cylinder the required database and the purple ellipses the output data.
GMAT is a more general tool that can be used for planning any space
mission. It is intended for mission optimization and mission analysis. By
modelling the asteroid whose trajectory we want to study as a satellite, we
can use GMAT to propagate its trajectory given initial positions and ve-
locities. GMAT uses data from JPL’s ephemerides for all the other bodies
28
considered in the model, i.e. it only propagates the position of the "satel-
lite". This makes for a fast-running propagator. GMAT is an open-source
system, allowing for a personalization of the software. GMAT has a user
interface that allows the user to select the desired coordinate system, initial
conditions, time span, force model and integrator.
STK is a tool used to model complex systems such as aircrafts, missiles

and satellites. It allows the user to visualize dynamic datasets in four di-
mensions (space and time). It has a component called "Astrogator" that
allows analysis in deep space. However, STK was built for missions close to
Earth. As a consequence, its deep space capabilities are far from being its
strong suit. As with HORIZONS, STK does not share its source code so no
user modifications are possible.
As real data (positions and velocities) for solar system bodies is hard to
obtain, we will be using HORIZONS as a reference and comparing all other
results from other propagators to the data generated by this system. Hence
HORIZONS data will, from now on, be considered as the truth.
In order to get an idea of the quality of each propagator, including PAT2 , we
will compare the propagated trajectories of a given asteroid, Apophis, over
the 10 year time span (Januray 1st 2014 to Januray 1st 2024). Note that the
close approach in 2029 is not included in this timespan. We mention this
fact because close approaches make for bigger numerical errors. In fact, a
close approach means a fast increase in acceleration which can be missed by
the numerical integration. However, if we are aware that a close approach
will happen during the numerical integration timespan, a tighter tolerance
could help in avoiding bad results. Note that this is not a property exclu-
sive to PAT2 . In fact, any numerical results obtained with an integrator not
specifically designed and built with fast changes in acceleration in mind (see
Appendix A) will suffer during close approaches. It is therefore important
to be aware of this issue and chose the appropriate tolerances.
Again, the trajectories will all be compared to those of HORIZONS

and the initial positions and velocities needed will be taken from the same
database. Note that HORIZONS is updated regularly and so initial condi-
tions queried from this system at different times may be different. As the
N-body problem is chaotic, these small changes in initial conditions may
have a big impact on the propagated trajectories.
Realistically, only GMAT would be an actual propagator candidate for the
study of deflection missions as it is the only one of the studied solutions
whose code is open source.
29
Runtime Error (last step) [km]

HORIZONS instantaneous -
GMAT a few seconds ≈ 5.4 · 106
PAT2 8.5 min. ≈ 310
Table 4.1: Benchmark comparisons using the 10 year timespan 2014-2024.
4.2 Benchmarking Results

The results obtained for the benchmarking problem using GMAT are shown
in Table 4.1 (section 5 provide more information on the definition of numer-
ical error). The results for STK are not shown here as obtaining meaningful
results seemed to be infeasible. In fact, it seems that STK uses an Earth-
centered coordinate system even when specified otherwise. Again, STK was
not built for deep space mission and therefore its use for analysis beyond
Earth is quite tricky. As this is not a viable solution for asteroid deflection
missions anyway, we will not be pursuing it any further in this context. A
similar remark can be made for GMAT. We believe that it is possible to
obtain better results using this tool. However, we were not able to do so
within a reasonable time frame and amount of work. In fact, in [19], the
authors claim to obtain 1500 kilometer errors with GMAT (and 8000 kilo-
meter errors within STK). However, with no details on the initial conditions
used, and the precise trajectories to which the output is compared to, it is
hard to reproduce similar results.
Again, if there were a close approach during the simulated timespan, the
errors would be bigger. As the ultimate goal is to have as small an error
as possible, we will not be using GMAT. Instead, PAT2 will be used for the
case studies in Section 5.
30
Section 5
Case Studies
Figure 5.1: Artist’s conception shows how families of asteroids are created. Over
the history of our solar system, catastrophic collisions between asteroids located in
the belt between Mars and Jupiter have formed families of objects on similar orbits
around the sun.
† Image taken from http://photojournal.jpl.nasa.gov/catalog/PIA17016
Now that the propagator PAT2 has been chosen, we will take a closer look
at the results obtained using it.
To do so, we will be taking a closer look at five different asteroids. These are
Apophis, Icarus, 2007 FT3, 2009 VZ39 and 2008 FF5 whose orbits are shown
in Figure 5.2. The last three of these have been picked from NASA’s Sentry
table [25] which lists the bodies with potential future Earth impact events.
Note that not all of the bodies on this list have diameters big enough for
them to enter Earth’s atmosphere and therefore they are not all hazardous
to us. The asteroids mentioned above are chosen because of the differences
31
SECTION 5. CASE STUDIES
Figure 5.2: Orbits of the asteroids considered for the case studies. These asteroids
are: Apophis, Icarus, 2007 FT3, 2009 VZ39 and 2008 FF5.
in three important orbital parameters: eccentricity, inclination with respect

to the ecliptic plane and size (which itself is related to mass). In fact, we
want to evaluate how PAT2 behaves with respect to different values in these
parameters in order to assess if it can be used on any asteroid and if some
types of asteroid parameters yield better numerical results.
The timespan studied here will be from January 1st 2014 to Januray 1st
2024, corresponding to 3652 Earth days.
Note that given the physical model and variables used in PAT2 , we are able
to study the trajectory of any of the objects considered i.e. the Sun. the 8
planets, Pluto and the Big 3. However, we will focus on that of the asteroids
because of the end goal of this thesis.
In this section as well, HORIZONS data will be considered as the truth.
Hence, from now on, the real positions and velocities will refer to those
queried from HORIZONS.
We suspect that bodies on orbits with small eccentricities will have bet-
ter results than those with high eccentricities, the key idea being that a
high eccentricity makes for big and fast changes in acceleration along the
trajectories which can be missed by the numerical integrator. We might also
suspect that high inclinations influence the accuracy of the solution.
32
As far as numerical errors are concerned, if t0 and tf are the initial

and final integration times respectively, rApo (tk ) is Apophis’ true position
(i.e. taken from HORIZONS) at time tk and r̃Apo (tk ) is the corresponding
numerical solution obtained with PAT2 , then the numerical error at time tk
(at output) is defined as
eApo (tk ) = kr̃Apo (tk ) − rApo (tk )k ∈ R ∀ t0 ≤ tk ≤ tf .
We will also use the component-wise error defined as
eApo (tk ) = r̃Apo (tk ) − rApo (tk ) ∈ R3 ∀ t0 ≤ tk ≤ tf .
Equivalent definitions will be used for the other asteroids studied here.
It is important to note that when using an integration scheme with error

control, the error (in each component, not the one defined above) at each
time step is computed. The idea is that for stable systems, controlling
the local error (error at each time step), we control the global error. In
fact, the local error is proportional to the global error for stable systems.
However, this is no longer true for unstable systems. This means that if
we have a stable system, we can check that the global error is consistent
with the prescribed tolerance. This cannot be done with unstable systems.
Nevertheless, error control is still important for unstable problems as we
know that small errors can explode.
5.1 Apophis
As stated previously, Apophis was put in the spotlight because of initial
observations and computations that indicated a high probability (around
2.7%) that it would collide with Earth. This lead to a big number of scientific
publications focusing on this asteroid. For this reason, Apophis will be the
object of our first case study.
Diam. Mass Orb. period Semi-maj. axis Inclin. Eccentr.

325 m 3.9 · 1010 kg 0.89 years 0.922 AU 3.3 deg 0.191
Table 5.1: Some of asteroid Apophis’ relevant physical data ([23]).
Some of Apophis’ relevant physical data is given in Table 5.1. Note

that its diameter of 325 meters would allow Apophis to survive entering
Earth’s atmosphere. According to Table 1.1, if a collision with Earth were to
occur, we could expect a large sub-global event causing around 500 thousand
33
Figure 5.3: Orbit of Apophis during the timespan 2014-2024 (data taken from
HORIZONS).
fatalities. Furthermore, with a small inclination of 3.3 degrees, Apophis’

orbit lies almost in the plane of the ecliptic along with the rest of the planets,
including Earth.
Apophis’ orbit is shown in Figure 5.3. We see that this orbit lies almost
in the ecliptic plane (which coincides with the xy plane).
Apophis’ trajectories are propagated using PAT2 over a 10 year period

(Januray 1st 2014-Januray 1st 2024) using an initial state vector extracted
from HORIZONS. The barycentric model is used here, i.e. for bodies with
moons (except Earth because its moon is explicitly included in the model)
we will consider the barycenters of the planet-moon systems, rather than
the body centers.
In Table 5.2, we compare the numerical errors in the last time step given
successively decreasing tolerances along with the marginal improvements
and associated runtimes. As explained in Section 3.3, the star symbol ?
indicates that the solver crashes given the tolerance of 10−3 due to the in-
herent instability of the N-body problem. Therefore, it is required to use a
tolerance smaller than 10−3 in order to solve this problem.
We expect to see a decreasing numerical error and an increasing runtime as
the tolerance is reduced. This is a typical tradeoff between accuracy and
runtime. In Table 5.2, we can see that by choosing a tolerance of 10−5 in-
34
Error (last step) [km] Marg. improv. [km] Runtime [min]

10−3 ? ? ?
10−4
Tolerance
310.0709 - 8.5
10−5 273.1395 36.9314 10.5
10−6 272.4092 0.7303 11.3
10−7 272.4018 0.0074 15.5
10−8 272.4020 −2 · 10−4 20.5
time step. Here, we consider the asteroid Apophis during the timespan 2014-2024.
marginal improvement of a smaller tolerance. The star symbol ? indicates that the
solver crashes given the specified tolerance
Figure 5.4: Numerical error in the position of Apophis during the timespan 2014-
2024. Numerical results are compared to HORIZONS data. The numerical error
is defined as eApo = kr̃Apo − rApo k where rApo is the true position and r̃Apo is the
numerical position of the asteroid. When rApo − r̃Apo is negative, we inverse eApo
so as to see the oscillatory behavior.
35
Figure 5.5: Numerical error in the position of Apophis during the timespan 2014-
2024 for each component x, y and z. Numerical results are compared to HORIZONS
data.
stead of 10−4 , we obtain a result that is approximately 37 kilometers closer

to the real solution in the last time step. Furthermore, we see that the run-
time is increased from 8.5 minutes to 10.5 minutes. We also notice that the
marginal improvement when going from a tolerance of 10−5 to 10−6 is very
small. However, the increase in runtime is not. This is a very poor tradeoff.
Given the results in this table, the tolerance 10−5 is preferred for Apophis
and we will therefore use this tolerance in the study of its trajectories.
The numerical error in the position of Apophis as a function of time is

shown in Figure 5.4, with time on the x axis and the error on the y axis. As
a reminder, the numerical error is defined as eApo = kr̃Apo −rApo k where rApo
is the true position and r̃Apo is the numerical position of the asteroid. When
r̃Apo − rApo is negative, we change the sign of eApo so as to see the oscillatory
behavior of the error. If we do not do this, all values of eApo are positive.
This oscillatory behavior suggests that there is something lacking in the
physical model - some additional gravitational term that is not included.
Interestingly, other studies and comparisons of asteroid propagators yield
similar oscillatory errors (e.g. [19]), suggesting that this is not an error
in PAT2 ’s code but rather in the physical description of the problem or
the equations of motion. None of the papers that present oscillatory errors
comment on this phenomenon.
A few more remarks about the error can be made here. First, the period of
36

10−4 1.505143·103 - 0.4
10−5 1.499333·103
Tolerance
5.8 0.5
10−6 1.498301·103 1.03 0.5
10−7 1.498274·103 2.7·10−2 0.7
10−8 1.498274·103 ≈0 0.9
Table 5.3: Numerical error in the position of Apophis using the Newtonian equa-
tions of motion (2.5).
oscillation is equal to the orbital period of the asteroid (0.89 Earth years).
Furthermore, the error seems to grow initially and then dampens after 10
years. A discussion with JPL’s Paul Chodas about this result concluded that
just as gravity pulls bodies together, it can also pull bodies apart. Therefore
it is possible for the errors to compensate for each other as it seems to be the
case here. However, we note that during most of the 10 years, the amplitude
of the error remains somewhat stable.
Figure 5.5 illustrates the numerical error in the position of Apophis, this
time decomposed into the x, y and z components. We observe that the
contribution from the z component is significantly smaller than that for the
two others. This is not a surprise as the orbit in question lies almost entirely
in the ecliptic plane (xy plane). Moreover, we notice that the x and y error
contributions have similar shapes. However, the error in the x component
drifts upward (i.e. increases). The maximal error over the 10 year timespan
is less than 500 kilometers which is small enough to study keyhole passages
and deflection missions.
For the sake of comparison, and in order to evaluate the impact of the
relativistic terms in equation (2.1), we propagate Apophis’ trajectory using
the Newtonian equations of motion (2.5). The numerical errors for different
tolerances are shown in Table 5.3 along with the marginal improvements
and associated runtimes. By comparing these results to those of Table 5.2,
we notice that the relativistic equations yield more precise trajectories, as
desired. For example, given a tolerance of 10−4 , the Newtonian equations of
motion yield an error in the last time step of 1.5·103 kilometers whereas the
relativistic version yields an error of about 310 kilometers, which is about 5
times smaller. However, the Newtonian equations, which are explicit, make
for significantly shorter runtimes.
An important concept in problems such as predicting orbits for celestial
37
Figure 5.6: Normalized simulated total energy of each body in the model as well
as that for the solar system as a whole. The normalization allows us to understand
how the total energy varies around its mean value.
38
Figure 5.7: Ratio of the simulated potential over kinetic energy for each body in
the model as well as the solar system as a whole.
39
Figure 5.8: Normalized total energy of each body in the model as well as that for
the solar system as a whole using data from HORIZONS. The normalization allows
us to understand how the total energy varies around its mean value.
40
Figure 5.9: Ratio of the potential over kinetic energy for each body in the model
as well as the solar system as a whole using HORIZONS data.
41
objects is that of energy. Figure 5.6 shows the normalized total (simulated)
energy of each body in the model as well as that for the entire solar sys-
tem, as a function of time. We decided to show values normalized by the
mean so as to illustrate the variation around the mean rather than the ac-
tual value. This is a simple post-processing exercise that can be done using
the outputted numerical positions and velocities. Ideally, we would like to
see a regular oscillation around a given value, indicating that the numerical
mechanical energy is approximately conserved. We check this because it is
well known that when solving the Kepler problem using the Euler Explicit
scheme, the trajectories spiral out (increase) while they spiral in (decrease)
when using the Euler Implicit scheme. These results were the reason for
developing the Euler Symplectic method, for which the trajectories neither
increase nor decrease.
For bodies such as Saturn, Uranus, Neptune and Pluto whose orbital peri-
ods are greater than 10 years (29 years for Saturn, 84 for Uranus, 165 for
Neptune and 247 for Pluto), a longer integration timespan could yield more
information. In fact, even if Uranus’ energy seems to be decreasing over a
10 year timespan, it may very well increase later on, once it has covered a
greater portion of its orbit. Hence if we are interested in studying the energy
of the outer solar system specifically, we could use a longer time span.
Figure 5.7 shows the ratio of potential over kinetic energy for each body
and for the solar system as a whole. For a perfectly circular orbit, we would
expect this ratio to be constant in time. However, due to eccentricities of
the orbits, there is an exchange between potential and kinetic energy which
translates to an oscillatory ratio, as seen in the figure.
For the sake of comparison, Figure 5.8 shows the normalized energy
obtained from HORIZONS data. The same timespan is used, i.e. Januray
1st 2014 to January 1st 2024. The desired oscillatory behavior is present for
the outer solar system bodies. However, it is not the case for the inner solar
system bodies, which also happen to be the smaller bodies. As stated in
Appendix A, there is more to the energy in the solar system than mechanical
energy. We also notice that the last four bodies (the Big3 and the Moon)
have relatively high variations, compared to the others. We should also note
that we are computing the energies for every day. By using a finer mesh, we
would see a more precise depiction of the evolution of the energy.
As stated previously, given the way PAT2 is designed, we have the pos-
sibility to study the positions (and velocities) of all of the bodies considered
in the model, not only the asteroid. Figures 5.10, 5.11, 5.12 and 5.13 il-
lustrate the numerical errors in the positions of the Sun and Jupiter both
in norm and component-wise. These two bodies are always very important
when modelling the solar system as they are extremely massive and there-
fore have a great influence on the other celestial bodies. The same graphs
42
Figure 5.10: Numerical error in the position of the Sun during the timespan 2014-
2024. Numerical results are compared to HORIZONS data.
Figure 5.11: Numerical error in the position of the Sun during the timespan 2014-
data.
43
Figure 5.12: Numerical error in the position of Jupiter during the timespan 2014-
2024. Numerical results are compared to HORIZONS data. Since Jupiter’s orbital
period is greater than 10 years, we might get more information by using a longer
timespan.
Figure 5.13: Numerical error in the position of Jupiter during the timespan 2014-
data. Since Jupiter’s orbital period is greater than 10 years, we might get more
information by using a longer timespan.
44
Figure 5.14: Numerical error in the position of Mercury during the timespan
2014-2024. Numerical results are compared to HORIZONS data.
Figure 5.15: Numerical error in the position of Mercury during the timespan 2014-
data.
45
Figure 5.16: Numerical error in the position of Venus during the timespan 2014-
Figure 5.17: Numerical error in the position of Venus during the timespan 2014-
data.
46
can be made for any of the 15 bodies considered. We notice that in both
cases, the x component of the error contributes more to the overall error
compared to the y and z components and this remark can be made for all
of the bodies in the model. As we are only studying a 10 year timespan, it
is interesting to do a similar analysis for bodies with shorter orbital periods.
We choose Mercury and Venus because they have the shortest orbital peri-
ods (88 days and 225 days, respectively) and therefore they orbit the Sun
a greater number of times compared to all the other bodies. Figures 5.14,
5.15, 5.16 and 5.17, give the numerical errors of Mercury and Venus both in
norm and component-wise. Again, for both these bodies, the x component
of the error dominates. We also notice a modulation of Mercury’s error,
which can again be associated with the fact that gravity not only pulls but
also pushes and therefore some errors might compensate for each other. Fi-
nally, Figure 5.15 is a great illustration of the instability of the problem, i.e.
how a small error in the state vector (after about 2500 days) can result in a
much greater error for future times.
5.2 Icarus
We will now consider the asteroid Icarus for our second case study. This
asteroid is chosen because it has a high eccentricity (0.8369) and a high
inclination with respect to the ecliptic plane (22.68 degrees). One of the
most important properties of a celestial body is its eccentricity. In terms
of numerical integration, a high eccentricity means that there will be fast
changes in the acceleration during each one of its orbits. If the integrator
does not catch these changes, the results can be very wrong. As a reminder,
an orbit with eccentricity e = 0 is circular. When 0 < e < 1 we have
an ellipse, when e = 1 a parabola and when e > 1 we have a hyperbola.
Furthermore, Icarus’ orbit crosses that of the Earth, reinforcing the interest
of its study.
A similar analysis to that for Apophis will be done for Icarus using PAT2 .
Diam. Mass Orb. period Semi-maj. axis Inclin. Eccentric.

1.4 km 1012 kg 1.12 years 1.078 AU 22.86 deg 0.8369
Table 5.4: Some of asteroid Icarus’ relevant physical data ([23], [24]).
Some of Icarus’ relevant physical data is given in Table 5.4.
successively decreasing tolerances. Marginal improvements and runtimes are
also shown. Given the results in this table, the tolerance 10−9 is preferred.
47

10−6 9.3369 · 103 - 13.1
10−7 8.5117 · 103
Tolerance
704.2484 15.7
10−8 269.4207 555.7565 19.7
10−9 241.2383 28.1824 23.2
10−10 231.1155 10.1228 27.4
time step. Here, we consider the asteroid Icarus during the timespan 2014-2024.
marginal improvement of a smaller tolerance.
Figure 5.18: Orbit of Icarus during the timespan 2014-2024 (data taken from
HORIZONS).
48
Figure 5.19: Numerical error in the position of Icarus during the timespan 2014-
We will therefore use a tolerance of 10−9 in the study of Icarus’ trajectories.

For the sake of comparison, we note that for Apophis, we chose a tolerance
of 10−6 and in fact, Apophis’ orbit is less eccentric than that of Icarus (0.191
for Apophis vs. 0.8369 for Icarus).
The numerical error in the position of Icarus is illustrated in Figure 5.19.

A few remarks can be made here. First, the error, as for Apophis, is peri-
odic. In fact, approximately every 408 days, there is a peak in the numerical
error and 408 days is the orbital period of Icarus. Second, the maximal error
adter 10 years is about 600 kilometers, which is below the target value of
1000 kilometers.
Figure 5.20 shows that the biggest errors (in black) happen when the as-
teroid is at perihelion, i.e. where it is closest to the Sun. This is because
at perihelion, the asteroid as a much greater acceleration than at any other
point on the orbit. Therefore, when approaching perihelion, the acceleration
increases quickly. Therefore the high eccentricity of this orbit is responsible
for the peeks in the numerical error. The points in blue in Figure 5.20 are
the points where the numerical error passes through the y = 0 line, i.e. when
the error is zero. We can see that this happens at aphelion, i.e. when the
asteroid is farthest from the Sun, i.e. when it is moving the slowest.
Moreover, the numerical error increases with time, as is usually the case in
numerical integration.
Figure 5.21 shows the x, y and z contributions to the numerical error. We
49
Figure 5.20: Points on Icarus’ orbit with biggest (in black) and smallest (in blue)
errors. The Sun is located at the right (invisible) focus of the ellipse, making the
left-most point of the ellipse the perihelion and the right-most point the aphelion.
Figure 5.21: Numerical error in the position of Icarus during the timespan 2014-
data.
50
Figure 5.22: Orbit of 2007 FT3 during the timespan January 1st 2014 - January
1st 2024 (data taken from HORIZONS).
can see that this time, the error in the z component contributes to the total
error, to the contrary of Apophis. This is to be expected as Icarus’ orbit is
inclined with respect to the ecliptic plane.
The numerical errors for the other bodies of the solar system are similar
to those obtained when studying the asteroid Apophis.
5.3 2007 FT3

We will now consider the asteroid 2007 FT3 for our third case study. This
asteroid, whose orbit is shown in Figure 5.22 is chosen because it has a high
inclination (almost 27 degrees) with respect to the ecliptic plane and a rel-
atively small eccentricity (0.307). Hence it shares similarities with Apophis
(w.r.t. eccentricity) and with Icarus (w.r.t. inclination). A similar analysis
to that for Apophis will be done for 2007 FT3 using PAT2 .
Diam. Mass Orb. period Semi-maj. axis Inclin. Eccen.

340 m 5.5 · 1010 kg 1.2 years 1.128 AU 26.83 deg 0.307
Table 5.6: Some of asteroid 2007 FT3’s relevant physical data ([25], [26]).
51

10−4
Tolerance
106.0741 - 8.2
10−5 90.1473 15.9268 9
10−6 90.6667 -0.5194 11
time step. Here, we consider the asteroid 2007 FT3 during the timespan 2014-2024.
Some of 2007 FT3’s relevant physical data is given in Table 5.6.
successively decreasing tolerances. Given the marginal improvements and
runtimes shown in this table, the tolerance 10−5 is preferred. Again, for the
sake of comparison, we note that for Apophis we chose a tolerance of 10−6
(small eccentricity) and 10−9 for Icarus (high eccentricity). Since 2007 FT3
has a relatively small eccentricity, this choice of tolerance corresponds to
our earlier hypothesis that small eccentricities will yield smaller numerical
errors. We will therefore use a tolerance of 10−5 in the study of 2007 FT3’s
trajectories.
Figure 5.23 shows the numerical error in the position of asteroids 2007
FT3 during the timespan Januray 1st 2014 to Januray 1st 2024. Similarly
to previous cases, the error is oscillatory with a period equal to the orbital
period of the body. The maximal error during the 10 year timespan is less
than 600 kilometers, which is acceptable for deflection mission analysis.
Figure 5.24 shows the decomposition of the error in the 3 component x, y
and z. This figure shows that all three of this components contribute to
the overall error, as expected because of the high inclination of the orbit.
Finally, we notice that the error increases with time, which is not surprising.
5.4 2009 VZ39

For our fourth case study, we consider the asteroid 2009 VZ39 whose orbit
is shown in Figure 5.25. This asteroid is chosen because it has both a small
inclination and a relatively small eccentricity. It is also the asteroid with
the highest number of potential impacts according to NASA’s Sentry list
[25] at this time. However, as this body’s diameter is well beneath the 140
meter threshold discussed previously, it would not survive entering Earth’s
52
Figure 5.23: Numerical error in the position of 2007 FT3 during the timespan
Figure 5.24: Numerical error in the position of 2007 FT3 during the timespan
2014-2024 for each component x, y and z. Numerical results are compared to
HORIZONS data.
53
Figure 5.25: Orbit of 2009 VZ39 during the timespan 2014-2024 (data taken from
HORIZONS).
atmosphere and therefore is not dangerous. Nevertheless, it is of academic

interest.
A similar analysis to that for Apophis will be done for 2009 VZ39 using
PAT2 .
Diam. Mass Orb. period Semi-maj. axis Inclin. Eccentr.

15 m 6.2 · 106 kg 1.81 years 1.4827 AU 2.52 deg 0.3824
Table 5.8: Some of asteroid 2009 VZ39’s relevant physical data. It is the aster-
oid with the highest number (815) of potential impacts with Earth, according to
NASA’s Sentry list ([25], [26]).
Some of 2009 VZ39’s relevant physical data is given in Table 5.8.
runtimes shown in this table, the tolerance 10−4 is preferred. Again, this
confirms that asteroids with low eccentricities will have smaller numerical
errors compared to asteroids with high eccentricities. We will therefore use
a tolerance of 10−4 in the study of 2009 VZ39’s trajectories.
Figure 5.26 shows the numerical error in the position of asteroid 2009
54

10−4
Tolerance
223.2047 - 6.8
10−5 226.9954 -3.7907 12
10−6 221.1609 5.8345 11
time step. Here, we consider the asteroid 2009 VZ39 during the timespan 2014-2024.
Figure 5.26: Numerical error in the position of 2009 VZ39 during the timespan
55
Figure 5.27: Numerical error in the position of 2009 VZ39 during the timespan
HORIZONS data.
Body Number of obs. used Solution number

Apophis 4077 JPL 190
Icarus 856 JPL 76
2007 FT3 14 JPL 6
2009 VZ39 8 JPL 8
2008 FF5 52 JPL 14
Table 5.10: Number of observations used in the statistical fit of an orbit (see
Appendix A) and solution number for each of the five asteroids considered. The
solution number represents the number of time that an asteroid’s orbit has been
calculated (this is done every time new observations are made or when there is a
change in the code).
VZ39 during the timespan January 1st 2014 to January 1st 2024. We notice
that this error looks more like that of Icarus rather than that of Apophis
or 2007 FT3. Figure 5.27 shows the decomposition of the error into its x,
y and z components. As for Apophis, 2009 VZ39 lies almost in the ecliptic
and therefore the z component of the error is very small.
The maximal error during this 10 year period is about 2000 kilometers,
which is not only above the desired maximal value of 1000 kilometers, but
is also the highest error out of the asteroids considered so far. This is a
56
surprising result. In fact, 2009 VZ39 has a moderate eccentricity and a very
small inclination and therefore we would expect better numerical results.
One reason for this higher than expected error could be the fact that there
high uncertainty about this body and its orbit, as illustrated in Table 5.10.
In this table we see that only a small number of observations is used in the
statistical data fit of the orbit. This is most likely indicates that there are
not many observations or measurement of this asteroid and therefore some
physical parameters might have high uncertainties. Moreover, the solution
number shown in this table represents the number of times a solution has
been computed, either because new observations have been made or there is
a change in the code (observation processing algorithm, for example). This
number should also give an idea as to how well the asteroid’s orbit is known
and again we see that this number is small compared to that for the other
bodies considered.
5.5 2008 FF5

For our last case study, we consider the asteroid 2008 FF5. This asteroid is
chosen because it has a very high eccentricity (0.965). Additionally, it has
a very big semi-major axis (almost 2.3 AU). Again, as this body’s diameter
is well beneath the 140 meter threshold discussed previously, it would not
survive entering Earth’s atmosphere and therefore is not dangerous. Never-
theless, it has interesting physical properties and therefore is of interest to
us.
A similar analysis to that for Apophis will be done for 2008 FF5 using PAT2 .
Diam. Mass Orb. period Semi-maj. axis Inclin. Eccentr..

81 m 7.1 · 108 kg 3.45 years 2.283 AU 2.6 deg 0.965
Table 5.11: Some of asteroid 2008 FF5’s relevant physical data ([26],[25]).
Some of 2008 FF5’s relevant physical data is given in Table 5.11.
runtimes, the tolerances 10−9 and 10−10 are the most interesting. We will
show the results of the numerical solution using 10−9 . We notice, once again,
that the tolerance needed in order to obtain reasonable errors is much smaller
than for less eccentric asteroids. However, it is important to note that
obtaining reasonable results is possible, but smaller tolerances are needed
and smaller tolerances translate to longer runtimes.
57
Figure 5.28: Orbit of 2008 FF5 during the timespan 2014-2024 (data taken from
HORIZONS).

10−6 3.8069·105 - 13.4
10−7 4.2613·104 3.38077·105
Tolerance
14.5
10−8 2.4135·103 4.0200·104 18.3
10−9 291.4986 2.1220·103 23.6
10−10 41.5419 249.9567 28.6
Table 5.12: Numerical error (in absolute value, compared to HORIZONS) in the
last time step. Here, we consider the asteroid 2008 FF5 during the timespan 2014-
2024. Each error (computed against HORIZONS) is then compared to that ob-
tained with the previously assigned tolerance (one order of magnitude less), yielding
the marginal improvement of a smaller tolerance.
58
Figure 5.29: Numerical error in the position of 2008 FF5 during the timespan
Figure 5.30: Points on 2008 FF5’s orbit with biggest (in black) errors. The Sun
is located at the right (invisible) focus of the ellipse, making the left-most point of
the ellipse the perihelion and the right-most point the aphelion.
59
Figure 5.31: Numerical error in the position of 2008 FF5 during the timespan
HORIZONS data.
Figure 5.29 shows the numerical error in asteroid 2008 FF5’s position.
This is the only one of the asteroids that we consider that does not have an
oscillatory error. This is also the body with the longest orbital period (3.45
years compared to a maximal orbital period of 1.8 years for the four others).
Figure 5.30 shows the point of maximal error (black point), which happens
at day 1636 (about four and a half years into the timespan). In this figure,
the Sun is located on the right-hand side of the ellipse and therefore the
maximal error happens at perihelion, similarly to Icarus as shown in Figure
5.20. Figure 5.31 shows the contribution of each component x, y and z to
the overall numerical error. As in previous cases, the error in z is extremely
small because this body almost lies in the ecliptic.
60
Section 6
Conclusions and Future

Work
Figure 6.1: The center of our galaxy, the Milky Way.

† Image taken from http://www.nasa.gov/multimedia/imagegallery/image_feature_649.html
The results obtained in Section 5 suggest that the propagator PAT2 can
be used in the study of asteroid deflection missions.
Table 6.1 shows the eccentricity and inclination of each of the five aster-
oids considered as well as the tolerances chosen for each body, the maximal
errors after 10 years and the runtimes. We notice that the maximal error
for the 10 year timespan of Januray 1st 2014 to Januray 1st 2024 ranges
from 489 kilometers (for Apophis) to 2.5 · 103 kilometers (for Icarus). Out
of the five asteroids, three have maximal errors under the threshold of 1000
61
SECTION 6. CONCLUSIONS AND FUTURE WORK
Figure 6.2: Tolerance versus runtime tradeoff.
Figure 6.3: Tolerance versus numerical error tradeoff (compared to HORIZONS).
62
Apophis Icarus ’07 FT3 ’09 VZ39 ’08 FF5

Eccentricity 0.191 0.8369 0.307 0.3824 0.965
Inclination [deg] 3.3 22.86 26.83 2.52 2.625
Tolerance 10−6 10−9 10−5 10−4 10−9
Max Error [km] 489.1 2.5·103 599.3 2·103 843.9
Runtime [min] 11.3 28.5 9 6.8 23.6
Table 6.1: Summary of the results obtained in Section 5. Max Error refers to the
maximal error over the 10 year timespan. These errors assume HORIZONS as the
truth. The values of the errors are rounded for simplicity.
kilometers. Moreover, the tolerances chosen for the asteroids with high
eccentricities (Icarus and 2008 FF5) are significantly smaller (3 orders of
magnitude) than those for the other asteroids. This is not surprising as
asteroids with eccentric orbits will have fast changes in acceleration. How-
ever, we also notice the accuracy of the solution does not only depend on
eccentricity. In fact, 2009 VZ39 has a relatively low eccentricity and has a
relatively big error. This can be attributed to the lack of knowledge that we
have about this asteroid and the uncertainty in the orbit in HORIZONS. In
fact, 2009 VZ39 has a small number of observations and therefore it’s input
data (initial state and mass) contains a lot of uncertainty, which can explain
the inaccuracy of the numerical solution and the orbit in HORIZONS has
only been computed a few times which means that it is not very precise.
Finally, a high inclination with respect to the ecliptic plane does not seem
to have an influence on the quality of the output.
Hence, when using PAT2 for orbit estimation, we can expect to obtain the
least accurate solutions when the asteroid in question has a high eccentric-
ity and a small number of observations or measurements. Nevertheless, in
other cases, the maximal error after 10 years should be smaller than 1000
kilometers.
Finally, we notice that the runtimes vary from less than 7 minutes (for 2009
VZ39, for which a tolerance of 10−4 is chosen) to 28.5 minutes (for Icarus,
for which a tolerance of 10−9 is chosen). This is a good illustration of the
impact of tolerance on runtime.
Figures 6.2 and 6.3 illustrate the relationships between tolerance and
runtime, and tolerance and error, respectively. As expected, runtimes in-
crease with decreasing tolerances and errors decrease with decreased toler-
ances. The rate of decrease in error depends on the orbit properties. These
two figures can be used when deciding what tolerance to use for a given
asteroid with a given orbit. Furthermore, they provide information on what
63
kind of errors can be expected.
As far as limitations go, we see here is that PAT2 can only be used on
planets and minor planets or asteroids but not active comets. The reason
is that active comets require additional orbital parameters to account for
the nongravitational perturbations caused by outgassing of the cometary
nucleus, which is not modeled here. Other than that, PAT2 is flexible and
can introduce more or less objects, as desired, with minimal effort.
The most important limitation, as mentioned previously, is related to the
limitation of the input data. The quality of the results for any N-body
problem will heavily depend on the input parameters such as the initial
state vector and the mass of each body. Since the problem is unstable,
even small errors in these parameters can (and will) lead to very bad re-
sults. For asteroids that are hard to observe, i.e. for which we have very
few measurements, we cannot expect good numerical results. Furthermore,
for applications where an extremely precise asteroid trajectory is necessary
(for rendez-vous missions, for example), PAT2 may not be the best solution.
However, for many applications, including deflection mission analysis, the
obtained accuracy is acceptable. Finally, when studying close approaches,
it is necessary to use smaller tolerances compared to those for periods with
no close approaches. The errors should be expected to be higher and the
runtimes longer.
Besides deflection mission analysis, PAT2 can be used in different con-

texts related to solar system dynamics.
In fact, using relativistic equations of motion and treating all bodies in the
system in the same way, it is possible to study each one of their positions and
velocities. Therefore, PAT2 can be used to generate planetary ephemerides
as well as asteroid ephemerides. Furthermore, the equations of motion can
be easily modified to include perturbing forces, i.e. to model all types of
deflection missions or even what if scenarios. The tool can also be used in
computing initial impact probabilities.
By allowing the PPN parameters β and γ to be modified (they are set to 1
in the general relativity theory), we have the possibility to study different
relativity theories. This is outside the scope of this project. However, it may
be interesting for physicists interested in alternative theories of relativity.
Future work related to the topics discussed in this thesis could include
further benchmarking against NASA’s GMAT. Studying periods of close
approach could yield important information as to the quality of each prop-
agator. Furthermore, some more digging into why the x component of the
numerical error dominates should be done. In theory, for asteroids that lie
(almost) in the ecliptic plane, the errors in x and y are expected to be similar
as the ecliptic plane corresponds to the xy plane. A deeper analysis of the
64
error for asteroid 2008 FF5 might also be interesting as it’s error was very
different compared to all the others presented.
Moreover, further work could include implementing the numerical integra-
tor in FORTRAN, a 128 bit language, from within MATLAB. This should
reduce the effect of roundoff errors and could potentially be beneficial to the
overall accuracy of the numerical solution.
Another path that could be explored in this context is that of develop-
ping a numerical integrator specifically built for orbit prediction, similarly
to JPL’s DIVA (see A). This would consist in finding a way to accurately
capture fast changes in accelerations which would then produce very high
quality ephemerides. We should note that DIVA was a multi-year endeavour
at JPL. Alternatively, licenses to DIVA can be purchased.
65
Appendices
66
Appendix A
JPL Visit Report
Context
Thanks to Professor Olivier de Weck and Doctor Paul Chodas, I had the
opportunity to spend a day at NASA’s Jet Propulsion Laboratory (JPL) in
Pasadena, California.
Dr. Paul Chodas is a senior scientist at JPL. He is the asteroid dynamics
and impact probabilities guru. He is the principal architect for orbit deter-
mination and ephemeris software used by the JPL NEO Program Office and
a co-developer of the Sentry impact monitoring system.
The goal of my visit was to discuss the work that I have done in the
context of my Master thesis as well as the HORIZONS system developed
by JPL. As HORIZONS was, and still is, a big project, there are a number
of people working on its different aspects, such as the mathematical model,
the numerical integrator, the database, the updating of data, etc. I also had
the opportunity to have one-on-one discussions about each of these topics
with their respective experts.
Expectations
Going into this incredible visit, I did not know what to expect. I had
prepared a PowerPoint presentation describing the work I had done until
then as well as some questions that I would have liked to get the answers to.
These questions were about the equations of motion that are used to generate
asteroid ephemerides as well as the numerical integrator and reference frame.
Furthermore, by reading some documentation about HORIZONS, I learned
that asteroids and planets are handled differently at JPL. This means that
the ephemerides for planets are generated separately from those for the
asteroids and these ephemerides are also stored separately. I therefore also
wanted to know why this distinction was made and what is the difference
67
APPENDIX A. JPL VISIT REPORT
in how those ephemerides are generated. Furthermore, as I had had the

experience of obtaining different positions for a given asteroid at a given
time with two separate queries from HORIZONS, I wanted to know more
about how, when and why the ephemerides are updated.
Outcome of the visit

The presentation that I gave went well. I was able to inform Paul on the
work I had done so far and the results that I had obtained. His first reaction
was to say that my results weren’t bad at all (I showed him Figures 5.4, 5.5,
5.19 and 5.21 from the case studies in Section 5). I described the model that
I have implemented as well as the equations of motion and the numerical
integrator.
Paul and I discussed Hamiltonian systems which are not implemented

in either the planetary or asteroid ephemerides. Paul remarked that when
considering the solar system as a whole, mechanical energy is not conserved.
There are a lot of other sources of energy that need to be taken into account
for a meaningful description of the exchange of energy.
The reason that the planets and asteroids are considered separately is
purely historical. In fact, asteroids have only started to be discovered re-
cently, long after trajectories of planets were being estimated. Therefore it
was easier to create a second database for asteroids with their own format
rather than trying to integrate these bodies with the planets.
Now, the way that the asteroid ephemerides are generated happens in a
few steps. First, the planetary ephemerides need to be known and therefore
previously computed. This is done using a numerical solver called DIVA.
This is a variable step size, variable order Adams method developed at JPL
over several years. It is specifically built to be able to handle fast changes
in accelerations, as is necessary when studying the elliptical orbits of plan-
ets. In fact, an object on an eccentric orbit will be greatly accelerated and
perihelion. Once positions are estimated, they are fit to data using the least
squares method. The bodies included in the model are the Sun and the
9 classical planets. Note that if DIVA had been developed later in time,
Pluto would most likely not have been included in the model. However, it
was considered a planet at that time and for this reason it was included in
the model. The equations of motion for the planets are the same as those
presented in equation (2.1), i.e. including the relativistic corrections asso-
ciated with the general theory of relativity. Interestingly, they use orbital
elements instead of Cartesian coordinates to describe the planets and aster-
oids’ ephemeris. The reason for this is that these orbital elements have a
68
geometrical meaning. A lot of information can be obtained simply by look-

ing at these, which is not possible with Cartesian coordinates. However, the
orbital elements are immediately converted into Cartesian coordinates for
the numerical integration.
Once the trajectories are estimated and fit to data, they are stored in the
form of polynomials so that it is possible to obtain the positions at any point
in time simply by evaluating the polynomial. Appendix B briefly explains
the type of interpolation used.
The way that the asteroid ephemerides are generated is entirely different.
In fact, all the planets’ trajectories are prescribed and only the trajectory
of the asteroid is computed (1 body problem). Because the ephemerides
are stored as polynomials, it is easy to obtain the necessary information
about any planet at any time. However, before this can be done for a
newly discovered object, a handful of astrometric observations need to be
made. Then a preliminary orbit is computed using a 2-body problem (i.e.
no perturbations from any of the planets is included). These are of course
very approximate, but this is ok. It allows the people who are tracking this
object to know where to look in the sky. Then, when observations have been
made over several days (thanks to the initial orbit estimation), a definitive
orbit can be computed, this time with planetary perturbations included.
The advantage of this method is that there is no (or minimal) error in
the positions of the planets and therefore the trajectories of the asteroid
cannot be any more accurate. The only errors are related to measurements
and numerical integration, contrarily to the method presented in Section 2
where the error in the positions of the planets will increase the errors in the
positions of the asteroid. Another advantage of this method is the reduced
computation time. In fact, the problem shrinks from a 15 body problem
to a 1 body problem and therefore there should be 15 times less variables
and function evaluations. Moreover, the equations of motion governing the
Earth-Moon system are extremely complicated; the tidal forces for example
have a big impact on the orbits. Therefore using precomputed orbits for the
Earth and Moon should produce better results.
Conclusions
There are two major conclusions to be made here. First of all, the way to
get the most accurate asteroid ephemerides is to use a previously generated
planetary ephemerides and solve a 1 body problem. This cuts down on run-
time tremendously.
The reason for the planet-asteroid separation is purely historic. However,
after a lot of trial and error, it was found that this separation is convenient
for asteroid ephemeris generating.
69
Second, a numerical integrator necessary for solving an N-body problem

(for the planets) accurately, is highly complex. Ordinary differential equa-
tion solvers are not equipped with handling fast changes in accelerations
that happen with elliptical orbits. There needs to be a special handling of
errors for these fast changes, especially because the system is unstable.
It is important to keep in mind the context of this thesis and the goal of
the propagator. Even if we do not have access to an extremely high qual-
ity numerical integrator as they do at JPL, using the method of prescribed
planetary positions, we should still obtain as estimate of the asteroid’s tra-
jectory that is accurate enough for the study of deflection missions.
Finally, a short remark can be made about the study of keyhole passages.
In fact, the definition of a keyhole is simply a region in space such that if an
asteroid were to pass through the keyhole, the asteroid would collide with
Earth. However, when studying deflection missions, the only important
thing to know is whether or not the asteroid will collide with Earth and
therefore any simulation should be run until it can be determined if this will
happen or not (i.e. we do not stop the integration at the keyhole but at the
impact). In fact, the only way we can know if a point in space is part of a
keyhole is to start an asteroid at that point and see if it then collides with
the Earth. It would be interesting to determine how much ∆V is needed to
avoid a particular asteroid hitting Earth. This ∆V will grow a lot in the
keyhole region. Note that in these situations, we consider all bodies as point
masses and therefore a point-mass asteroid is either in or out of a keyhole,
but cannot be partly in and partly out as would be the case in reality.
70
Appendix B
Chebyshev Polynomial
Interpolation
In order to minimize the interpolation error (in the norm of the maximum,
i.e. the maximum error), we compute the roots of the Chebyshev polyno-
mials and use them as the nodes for a Lagrangian interpolation polynomial
(or any other interpolation polynomial of the same degree). For a general
interval [a, b] and a function f that is interpolated by a polynomial Pn of
degree n, by using
b+a b−a (2i − 1)π

xi = + cos
2 2 2(n + 1)
for i = 1, . . . , n+1, the interpolation error that occurs is bounded as follows:

n+1
b−a Bn+1

|f (x) − Pn (x)| ≤ n
,
2 2 (n + 1)!
with Bn+1 = maxx∈[a,b] |f (n+1) (x)|. More details can be found in [27].
71
Appendix C
How Asteroids Get Their

Names
Reading Section 5, some may wonder how these asteroids get their names.
The answer is far from obvious.
First, we should note that asteroids are being discovered on a daily basis
and will keep on being discovered for decades and probably more. These
asteroids are the smaller ones, that do not pose much of a threat to Earth.
The larger asteroids are quite well known. In fact, we have information on
all of these potentially dangerous asteroids.
Once an asteroid is discovered, it is given a name. The first four digits of

this name correspond to the year that the asteroid was discovered. The first
letter after these four digits corresponds to half years. Let me explain. A
corresponds to the first half of January, B to the second half of January, C
to the first half of February, and so on. The rest of the name comes from the
fact that a great number of asteroids are being discovered and therefore we
need a second letter and some more numbers to distinguish all the asteroids
discovered in a given half year.
For example, 2007 FT3 was discovered in the year 2007. More precisely, in
the second half of March (letter F).
72
Bibliography
[1] NASA, 2006 Near-Earth Object Survey and Deflection Study. Technical
report, DRAFT: Pre-Decisional Material, 2006
[2] J. Laskar, M.Gastineau Existence of collisional trajectories of Mercury,

Mars and Venus with the Earth. Nature, Vol 459, 2009
[3] Q. Wang, The existence of global solution of the N-body problem. Chin.
Astron. Astrophys. 10, 1991
[4] Victor Brumberg, On derivation of EIH (Einstein–Infeld–Hoffman)

equations of motion from the linearized metric of general relativity the-
ory. Springer Science+Business Media B.V. 2007
[5] E. Myles Standish, James G. Williams, Orbital Ephemerides of the Sun,

Moon, and Planets. Chapter 8
[6] Alan Pitz, Christopher Teubert, Bong Wie Earth-Asteroid Im-

pact Probability Computations of Disrupted Asteroid Fragments using
GMAT/STK/CODES. AAS 11-408
[7] Ernst Hairer, Christian Lubich, Gerhard Wanner, Geometric Numerical

Integration Illustrated by the Stoermer–Verlet Method. Cambridge Uni-
versity Press, 2003, pp. 399–450
[8] Ernst Hairer, Christian Lubich, Gerhard Wanner, Geometric Numerical

Integration: Structure-Preserving Algorithms for Ordinary Diff erential
Equations. 2004 Springer seies in computational mathematics 2nd Edi-
tion,
[9] Estabrook, F. B, Derivation of Relativistic Lagrangian for n-Body

Equations Containing Relativity Parameters β and γ. JPL D-18611,
JPL Interoffice Memo (internal document) Jet Propulsion Laboratory,
Pasadena, California, June 14, 1971
[10] François E. Cellier, Ernesto Kofman, Continuous System Simulation.

Springer US, 2006
73
BIBLIOGRAPHY
[11] E. V. Pitjeva, High-Precision Ephemerides of Planets—EPM and De-

termination of Some Astronomical Constants. Solar System Research,
Vol. 39, No. 3, 2005, pp. 176–186. Translated from Astronomicheskii
Vestnik, Vol. 39, No. 3, 2005, pp. 202–213.
[12] L. F. Shampine, Limiting Precision in Differential Equation Solvers.

Mathematics of computation, Vol. 28, No. 125, January, 1974
[13] John H. Mathews, Kurtis K. Fink, Predictor-Corrector Methods. Nu-

merical Methods Using Matlab, 4th Edition, 2004, Section 9.6
[14] L. F. Shampine, M. K. Reichelt, The MATLAB ode suite
[15] Florin Diacu, Solution of the n-body problem. The Mathematical Intel-
ligencer, Vol. 18, No. 3, 1996
[16] Cellier, François E., Kofman, Ernesto, Basic Principles of Numerical

Integration. Continuous System Simulation, Chapter 2, Springer, 2006
[17] V. A. Brumberg Essential General Relativity Theory. Institut de

Mecanique Celeste et de Calcul des Ephemerides, Astronomy and Celes-
tial Mechanics Course, Lecture 3
[18] Floating-Point Arithmetic in MATLAB. http://people.uib.no/

fra001/radu/num2_bergen/float.pdf
[19] Alan Pitz, Chritopher Teubert, Bong Wie Earth-impact probability com-
putation of disrupted asteroid fragments using GMAT/STK/CODES.
AAS 11-408
[20] http://earthguide.ucsd.edu/virtualmuseum/ita/08_1.shtml
[21] https://solarsystem.nasa.gov/planets/profile.cfm?Object=
Dwa_Ceres
[22] http://neo.jpl.nasa.gov/apophis
[23] http://ssd.jpl.nasa.gov/?horizons/sbdb.cgi#top
[24] http://nssdc.gsfc.nasa.gov/planetary/factsheet/
asteroidfact.html
[25] http://neo.jpl.nasa.gov/risk/
[26] http://ssd.jpl.nasa.gov/sbdb.cgi
[27] http://www.math.wsu.edu/faculty/genz/448/lessons/l303.pdf
74

Thesis Patricia Egger

Uploaded by

Copyright:

Available Formats

Thesis Patricia Egger

Uploaded by

Document Information

Original Description:

Copyright

Available Formats

Share this document

Share or Embed Document

Sharing Options

Did you find this document useful?

Is this content inappropriate?

Copyright:

Available Formats

Thesis Patricia Egger

Uploaded by

Copyright:

Available Formats

Ecole Polytechnique Fédérale de

Efficient computational approaches

Patricia Egger Under the supervision of:

July 14, 2014

Thanks so much to all of you! I will never forget this experience.

1 Introduction and Motivation 1

1.1 Typical impact types, frequency and conesquences . . . . . . . . . 2

3.1 Positions, velocities and accelerations in units of seconds and

4.1 Benchmark comparisons . . . . . . . . . . . . . . . . . . . . . 30

5.1 Some of asteroid Apophis’ relevant physical data . . . . . . . 33

6.1 Summary of results . . . . . . . . . . . . . . . . . . . . . . . . 63

1.2 Differences between orbits of a typical NEA and PHA . . . . 3

2.2 Illustration of the ecliptic plane . . . . . . . . . . . . . . . . . 7

3.2 Flow chart describing the inner workings of PAT2 . . . . . . . 21

4.2 Flow chart describing the inner workings of HORIZONS . . . 28

5.2 Orbits of the asteroids considered for the case studies . . . . 32

5.27 Component-wise error in the position of 2009 VZ39 . . . . . . 56

6.2 Tolerance versus runtime tradeoff. . . . . . . . . . . . . . . . . . 62

Aphelion: Defined as the point on a planet’s orbit where it is at

Asteroid: Also known as a minor planet, this is a large rocky body in

AU (Astronomical Unit): Is the mean distance between the Earth

Barycenter: Of a body, is its center of mass. Since the Sun makes

Deflection mission: A technological solution for deviating the tra-

∆V : Delta V: Change is velocity (denoted V ) needed to modify a

Ecliptic, or ecliptic plane: Defined as the imaginary plane contain-

Ephemeris: Gives the positions of naturally occurring astronomical

Ephemerides: Plural of the ephemeris.

Hamiltonian: Mathematical formalism to describe the equations of

ICRF (International Celestial Reference Frame): With its ori-

NEO - respectively NEA (Near-Earth Object - resp. Near-

Orbital eccentricity: Of a celestial object, is a real and positive

Perihelion: Defined as the point on a planet’s orbit where it is closest

PHO - respectively PHA (Potentially Hazardous Object -

PN formalism: Post-Newtonian formalism - In situations where

PPN formalism (Parametrized Post-Newtonian formalism):

Propagator: Simulates the trajectory of one or more objects over

Introduction and Motivation

Event type Diameter Fatalities Impact interval

Table 1.1: Typical impact types, frequency and conesquences. [1]

The remaining part of this thesis is structured as follows. Section 2

Physical and Mathematical

Figure 2.1: Illustration of the N-body problem.

This section will provide an overview of the physical and mathematical

increasing N, we obtain more precise trajectories.

2.1 The N-body problem

Even though the N-body problem is unstable, a convergent expansion

"The convergence of my power series is admittedly unsatisfactory, and so

It is noteworthy that in his proof, Wang considers bodies as point masses

transformation in order to prove the existence of a solution to the N-body

2.1.1 Equations of motion

For each body i, the point mass acceleration is given by:

is the state vector of the N bodies considered, then

M (y)y 0 = f (t, y) (2.3)

In the context of numerical integration, explicit problems are preferred.

For the sake of comparison, we also introduce the Newtonian equations of

Additional forces may be included in the physical model in order to

• Solar radiation pressure or SRP:

• Perturbations from the less massive asteroids in the belt, modeled by

• Perturbations from the 21 largest Trans-Neptunian Objects (or TNOs):