Effect of Coordinate Switching on Simulation Accuracy of
Translunar Trajectories
Except where reference is made to the work of others, the work described in this
thesis is my own or was done in collaboration with my advisory committee. This
thesis does not include proprietary or classi ed information.
Mana P. Vautier
Certi cate of Approval:
John E. Cochran, Jr.
Department Head and Professor
Aerospace Engineering
Andrew J. Sinclair, Chair
Assistant Professor
Aerospace Engineering
David A. Cicci
Professor
Aerospace Engineering
Joe F. Pittman
Interim Dean
Graduate School
Effect of Coordinate Switching on Simulation Accuracy of
Translunar Trajectories
Mana P. Vautier
A Thesis
Submitted to
the Graduate Faculty of
Auburn University
in Partial Ful llment of the
Requirements for the
Degree of
Master of Science
Auburn, Alabama
09 August, 2008
Effect of Coordinate Switching on Simulation Accuracy of
Translunar Trajectories
Mana P. Vautier
Permission is granted to Auburn University to make copies of this thesis at its
discretion, upon the request of individuals or institutions and at
their expense. The author reserves all publication rights.
Signature of Author
Date of Graduation
iii
Vita
Mana Vautier, son of Ian Philip and Hinekaitangi Awhi (Skipwith) Vautier, was
born 5 February, 1980, in Auckland, New Zealand. He graduated from Saint Kentigern
College in 1997 and attended the University of Auckland for one year. In 2003
he transferred to Brigham Young University in Provo, Utah and graduated with a
Bachelor of Science degree in Physics-Astronomy. He began his graduate studies
in Aerospace Engineering at Auburn University in August 2006. In August 2004
he married Annette (Hughes) Vautier, daughter of Eugene Marion and Constance
(Hu man) Hughes. They have two sons, Michael and Benjamin.
Figure 1: My Family.
iv
Thesis Abstract
Effect of Coordinate Switching on Simulation Accuracy of
Translunar Trajectories
Mana P. Vautier
Master of Science, 09 August, 2008
(B.S., Brigham Young University, 2006)
70 Typed Pages
Directed by Andrew J. Sinclair
This thesis focuses on the e ect of round-o error in the accurate simulation of
translunar trajectories. The three-body dynamics can be posed in either an Earth-
centered (EC) or Moon-centered (MC) frame. In this study, multiple translunar
trajectories were simulated to determine if there is an optimal switch point from an EC
to MC frame that minimizes round-o error. A high delity baseline simulation was
rst created, and the entire trajectory was propagated in EC. Comparison trajectories
were then simulated at lower precision using both EC and MC frames. The trajectory
was propagated rst in EC and then switched to MC at a preselected switch point.
By testing a range of switch points, it was determined that switching to the MC frame
during the rst 10% of the trajectory led to signi cantly higher round-o errors. For
any later switch points, there was little sensitivity in the round-o error. There was,
however, some correlation to switch points located at the radii of well established
gravitational spheres of the Moon. The importance of round-o error in simulation of
translunar trajectories is also demonstrated by calculation of the Lyapunov exponents
along the trajectory.
v
Acknowledgments
I would like to acknowledge all my professors who helped me get to where I am
today. I am especially grateful to my advisor, Dr. Andrew Sinclair, for his constant
encouragement, guidance and support. He went above and beyond to help me succeed
in my graduate studies, while still having fun with the BBQ?s, the Biscuits and the
Wartburg Warts (who came up with that name anyway?). I also want to acknowledge
Odyssey Space Research, L.L.C. who took me on as an intern during the summer of
2007, during which time the idea for my research was rst inspired. I must also
mention my fellow graduate students with whom I shared the last two years. Their
friendship, humour and timely help made even the toughest assignments bearable.
Much deserved gratitude also goes to my parents and family for their unending love
and support throughout my entire career as a full-time student (all 20 years worth!).
My parents instilled in me the belief that I could do anything if I worked hard and
truly put my mind to it, and my siblings never failed to help me feel brainy while
making sure my head never got too big. I want to express my love and gratitude to
Michael and Benjamin, two of the best sons a father could hope for. Their infectious
laughter, their love of life, and their joy in the simple things served as a constant
reminder as to what is most important in life. Finally, I want to acknowledge my
sweet angel wife, Annette. She is truly a light in my life, and I am forever grateful for
her continual patience, encouragement and love. She has stood by my side through
the thick and the thin, the good and the bad, the easy and the rough. I take much
joy in knowing that together we stand side-by-side as eternal companions through
this great journey called life.
vi
Style manual or journal used Journal of Approximation Theory (together with the
style known as \aums"). Bibliograpy follows van Leunen?s A Handbook for Scholars.
Computer software used The document preparation package TEX (speci cally
LATEX) together with the departmental style- le aums.sty.
vii
Table of Contents
List of Figures x
List of Tables xi
1 Introduction 1
2 System Model 3
3 Gravitational Spheres 10
3.1 Sphere of Gravitation . . . . . . . . . . . . . . . . . . . . . . . . . . . 11
3.2 Sphere of In uence . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12
3.3 Hill Sphere . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15
3.4 Kislik Sphere . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 17
4 Chaos 20
4.1 Historical Perspective . . . . . . . . . . . . . . . . . . . . . . . . . . . 20
4.2 Chaos Theory and Dynamic Systems . . . . . . . . . . . . . . . . . . 22
4.3 Lyapunov Exponents . . . . . . . . . . . . . . . . . . . . . . . . . . . 24
4.4 State Transition Matrix . . . . . . . . . . . . . . . . . . . . . . . . . 29
5 Simulation Model 31
6 Results 38
7 Discussion 46
7.1 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 46
7.2 Future Work . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 47
Bibliography 48
Appendices 50
A MATLAB Code 51
A.1 Double Precision Simulation . . . . . . . . . . . . . . . . . . . . . . . 51
A.2 Single Precision Simulation . . . . . . . . . . . . . . . . . . . . . . . . 52
A.3 Sub-Routines . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 54
viii
B Historical Review 56
B.1 Claudius Ptolemaeus . . . . . . . . . . . . . . . . . . . . . . . . . . . 56
B.2 Nicolaus Copernicus . . . . . . . . . . . . . . . . . . . . . . . . . . . 56
B.3 Galileo Galilei . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 56
B.4 Johannes Kepler . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 56
B.5 Isaac Newton . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 57
B.6 Joseph Louis Lagrange . . . . . . . . . . . . . . . . . . . . . . . . . . 57
B.7 Pierre-Simon Laplace . . . . . . . . . . . . . . . . . . . . . . . . . . . 57
B.8 Henri Poincar e . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 58
B.9 Aleksandr Lyapunov . . . . . . . . . . . . . . . . . . . . . . . . . . . 58
Index 59
ix
List of Figures
1 My Family . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . iv
2.1 Position vectors . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8
3.1 Sphere of In uence size variance . . . . . . . . . . . . . . . . . . . . . 14
3.2 Zero velocity curves . . . . . . . . . . . . . . . . . . . . . . . . . . . . 16
4.1 Evolution of Lorenz?s experiment . . . . . . . . . . . . . . . . . . . . 21
4.2 Bifurcation diagram for the Logistic map . . . . . . . . . . . . . . . . 23
4.3 The Lorenz Attractor . . . . . . . . . . . . . . . . . . . . . . . . . . . 24
4.4 Evolution of a hypersphere . . . . . . . . . . . . . . . . . . . . . . . . 27
5.1 System model position vectors . . . . . . . . . . . . . . . . . . . . . . 31
5.2 Translunar trajectory . . . . . . . . . . . . . . . . . . . . . . . . . . . 33
5.3 Schematic of initial conditions . . . . . . . . . . . . . . . . . . . . . . 34
6.1 Acceleration ratios . . . . . . . . . . . . . . . . . . . . . . . . . . . . 38
6.2 RSS position error . . . . . . . . . . . . . . . . . . . . . . . . . . . . 39
6.3 RSS error for three switch points . . . . . . . . . . . . . . . . . . . . 40
6.4 RMS values . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 41
6.5 Magni ed RMS values . . . . . . . . . . . . . . . . . . . . . . . . . . 42
6.6 Representation of simulation results . . . . . . . . . . . . . . . . . . . 43
6.7 Lyapunov Exponents . . . . . . . . . . . . . . . . . . . . . . . . . . . 44
6.8 Spacecraft range . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 45
x
List of Tables
3.1 Radii of Gravitational Spheres . . . . . . . . . . . . . . . . . . . . . . 19
5.1 Physical Parameters of the Simulation Model . . . . . . . . . . . . . . 35
6.1 Gravitational Spheres and Switch Point Times . . . . . . . . . . . . . 43
xi
Chapter 1
Introduction
NASA?s Project Constellation is committed to the long-term human and robotic
exploration of the Moon, Mars, and beyond. Although the spacecraft developed will
still have the capability of delivering payload to Low Earth Orbit (LEO) , there is a
renewed emphasis on leaving the gravitational in uence of the Earth, and traveling to
the Moon. Compared to the era of the Apollo missions, numerical simulations have
grown in capability and importance. As a result, simulation accuracy has also im-
proved signi cantly, and numerical results can now be calculated with much greater
precision. Until recently, these more advanced tools have primarily been used in sim-
ulating LEO?s. The new focus on lunar missions, however, motivates the development
of precise translunar simulations.
The issue of model choices a ects the numerical error in the simulation of translu-
nar trajectories. Numerical errors typically include several sources such as discretiza-
tion error, round-o error, and iterative error. This work presents an investigation
of the role of coordinate choice in reducing round-o error to improve simulation
accuracy. The importance of round-o error in the simulation accuracy will also be
investigated by demonstrating the chaotic nature of the trajectories.
Chapter 2 discusses the three-body problem and representations of the spacecraft
trajectory. Chapter 3 presents an overview of four de nitions for a gravitational
sphere. Chapter 4 introduces the concept of chaos theory and the sensitivity of a
trajectory to initial conditions. The speci c model used for the simulations in this
study are presented in Chapter 5, and results are given in Chapter 6. Finally, a
1
discussion of possible future work is found in Chapter 7 along with some concluding
remarks.
2
Chapter 2
System Model
Published in 1687, Sir Isaac Newton?s Philosophi Naturalis Principia Mathe-
matica or more simply, the Principia, is considered by many historians to be one of
the supreme achievements of the human mind [1]. In it, he introduced the three laws
of motion which help describe most non-relativistic dynamical systems.
By letting r, v and a respectively denote the position, velocity and acceleration
of a body of mass m from a xed origin so that
v = drdt and a = dvdt = d
2r
dt2; (2.1)
Newton?s rst and second laws can be summarized by the mathematical expression
F = d(mv)dt (2.2)
where F is the impressed force acting on the body. Assuming that the body?s mass
is constant, equation (2.2) can then be generalized for the case where n multiple
external forces are acting on the system producing
nX
i=1
Fi = mid
2ri
dt2 : (2.3)
Also the third law, again in vector notation, can be expressed mathematically by
F1 = F2: (2.4)
3
These three laws of motion enunciated by Newton laid the foundations of the eld of
dynamics.
Another signi cant scienti c contribution made by Newton was his formulation
of the law of universal gravitation. This law can be stated quite simply that every
particle in the universe is simultaneously attracted to every other particle in the
universe with a gravitational force which is directly proportional to the product of the
masses and inversely proportional to the square of the distance between them. Thus,
for any two bodies, we can express the magnitude of this force by the relationship
F = Gm1m2r2 (2.5)
where F is the magnitude of the gravitational force of attraction between masses m1
and m2 which are separated by a distance r. The proportionality factor, G, is known
as the constant of universal gravitation and has a value of (6:6741 0:0008) 10 11
m3=(kg s2) [2]. Although this law holds true speci cally for point mass particles,
it also holds for spherical bodies whose masses are evenly distributed. Since most
planets and stars can be approximated as such, Newton?s law of universal gravitation
has formed the basis of celestial mechanics, as well as the more modern eld of
astrodynamics.
One of the most fundamental problems of astrodynamics is to accurately describe
the motion of n gravitationally interacting massive particles. Obviously the most
simplistic scenario occurs when n = 2, and is known as the two-body problem. This
problem states, ?Given the position and velocity of two massive particles at any time,
calculate their position and velocity for any other time.? The signi cance of this
problem lies in the two main facts that a) it is the only problem in gravitational
4
dynamics for which a general analytical solution is obtainable, and b) it can be used
as an approximation for a wide variety of multi-body orbital problems.
A more general n-body problem can be stated by considering a system of n point
masses m1;m2;m3;:::;mn. An analytical description of the interactions and resulting
motion of these bodies can be provided using the four previously mentioned laws
published by Newton. Assume a suitable inertial coordinate system in which the
positions of the n masses are given by
ri = xiix +yiiy +ziiz (i = 1;2;3;:::;n)
where ix, iy, and iz represent the respective unit vectors along thex, y andz coordinate
axes. Applying Newton?s law of universal gravitation, the magnitude of the force of
attraction exerted on mi by mj is
Fi = Gmimjr2
ij
(i = 1;2;3:::n j = 1;2;3;:::;n j6= i) (2.6)
where
rij =jrj rij= p(rj ri) (rj ri).
The direction of the force can be conveniently expressed by multiplying equation
(2.6) by the unit vector directed along the line of centers of the two masses mi and
mj to obtain
Fi = Gmimjr2
ij
rij
jrijj (2.7)
where
rij = rj ri = rji.
5
Then, by summing all such gravitational forces acting on the ith body
Fi = Gmim1r3
i1
ri1 +Gmim2r3
i2
ri2 +::::+Gmimnr3
in
rin: (2.8)
Since the ith body cannot exert a force upon itself, equation (2.8) does not include
the term Gmimir3
ii
rii and can be simpli ed using the summation notation such that
Fi = G
nX
j=1
mimj
r3ij rij (i=1,2,3,...,n and j6= i): (2.9)
Now applying Newton?s second law of motion
mid
2ri
dt2 = G
nX
j=1
mimj
r3ij rij (i=1,2,3,...,n and j6= i): (2.10)
Finally, canceling mi from both sides gives the second order, non-linear, vector dif-
ferential equation of motion for the ith body of the system of n particles all acting
under their mutual gravitational attractions.
d2ri
dt2 = G
nX
j=1
mj
r3ij rij (i=1,2,3,...,n and j6= i): (2.11)
By setting n = 2, it is a relatively simple matter to reduce equation (2.11) to an
expression for only two bodies. When coupled with a set of initial conditions for the
positions and velocities at time t0, the motion of the two bodies is fully described,
and the positions and velocities at any future time t may be calculated. This forms
the solution to the famous two-body problem stated and solved by Newton.
For the case where n > 2 however, a solution is not so easily obtained. By
the end of the nineteenth century, it was considered an important and challenging
enough question that for his 60th birthday, Oscar II, King of Sweden, announced that
a prize would be given to anyone who could nd the solution to the problem. The
6
speci c announcement was: \Given a system of arbitrarily many mass points that
attract each according to Newton?s law, under the assumption that no two points
ever collide, try to nd a representation of the coordinates of each point as a series in
a variable that is some known function of time and for all of whose values the series
converges uniformly" [3].
In case the problem could not be solved, any other important contribution to
celestial mechanics would then be considered to be prizeworthy. The prize was nally
awarded to Henri Poincar e, even though he did not solve the original problem. One
of the judges, the distinguished Karl Weierstrass, said, \This work cannot indeed be
considered as furnishing the complete solution of the question proposed, but that it
is nevertheless of such importance that its publication will inaugurate a new era in
the history of celestial mechanics." The seminal work produced by Poincar e would
later lead to the development of chaos theory as discussed in Chapter 4.
It will be instructive here to develop the di erential equations of relative motion
for the three-body problem. We have from equation (2.11) the relation
d2ri
dt2 = G
nX
j=1
mj
r3ij rij (i=1,2,3...n and j6= i) (2.12)
Now let the reference body be mass m1. Its equation of motion with respect to an
inertial (unaccelerated) coordinate frame is then
d2r1
dt2 = G
m2
r312 r12 +G
m3
r313 r13; (2.13)
while the equation of motion for m2 in the same coordinate axes is
d2r2
dt2 = G
m1
r321 r21 +G
m3
r323 r23: (2.14)
7
Figure 2.1: System model position vectors.
By subtracting equation (2.13) from (2.14), the motion of m2 relative to m1 is
obtained
d2r2
dt2
d2r1
dt2 = G
m
1
r321 r21
m2
r312 r12 +
m3
r323 r23
m3
r313 r13
: (2.15)
Now de ne the following variables
r r2 r1 r3 r1 d r2 r3 (2.16)
as shown in Fig. 2.1 where r1, r2 and r3 are position vectors from an inertial reference
point, O, to the corresponding masses m1, m2 and m3. By factoring appropriate terms
and noting that rij = rj ri, equation (2.15) can be simpli ed to
d2r
dt2 +
G(m1 +m2)
r3 r = Gm3
1
d3d +
1
3
: (2.17)
Following a similar process, the motion of m2 relative to m3 can shown to be
8
d2d
dt2 +
G(m2 +m3)
d3 d = Gm1
1
r3r
1
3
: (2.18)
In both instances, the second term on the left hand side corresponds to the cen-
tral acceleration due to the primary body, while the term on the right hand side is
associated with the disturbing force from the perturbing body.
It is also worthwhile to mention that there exists many other external forces
that may act upon any given body, other than the gravitational forces described
above. These forces, known as perturbations, may include, but are not limited to,
atmospheric drag, thrusting, solar radiation pressure, outgassing and charged or un-
charged particles. One of the most signi cant perturbations is due to non-sphericity
of the bodies. As noted previously, Newton?s law of gravitation applies only if the
mass of a body is evenly distributed in homogenous spherical shells. For this reason,
the expressions derived above are approximations, and the inclusion of additional
perturbing forces could more accurately model the true trajectory.
9
Chapter 3
Gravitational Spheres
There are many applications in which the n-body problem may be broken down
into a series of two-body problems, particularly for the case when the mass of the
primary body is much greater than the mass of all other bodies considered in the
system. Implementing this method provides an appropriate rst approximation to
the determination of precise spacecraft orbits.
In many such cases, the true trajectory of the spacecraft is represented by an
approximate one composed of a sequence of conic sections, constructed by deriving
solutions of individual two-body problems in a particular sequence throughout the
trajectory. For example, if a spacecraft were traveling from the Earth to Mars, the
journey can be modeled by rst treating the trajectory as a two-body system com-
posed of the Earth and the spacecraft. At some point, the trajectory is modeled
by switching to the spacecraft-sun two-body system, and nally as the spacecraft
approaches Mars, the solution to the spacecraft-Mars two-body system is used.
This particular method is known as the patched-conic approximation, and despite
the fact that it is only an approximation to true trajectory, it does a ord a convenient
means of investigating various initial and boundary conditions, while avoiding lengthy
computation times. The obvious question then, is at what point along the trajectory
should the reference frame be switched such that the terminal position and velocity
of one section is equivalent to the initial position and velocity of the next section.
The concept of a gravitational sphere introduces the notion that there is a certain
region of space associated with any given body. In celestial mechanics, there have been
10
many di erent approaches to constructing these gravitational spheres, the boundaries
of which can be examined as a possible location for an optimal switch point.
3.1 Sphere of Gravitation
The de nition of a gravitational sphere most directly related to the gravitational
forces between two bodies is that derived by Gurzadyan [4]. Consider two celestial
bodies of masses m1 and m3. It is evident that a spacecraft of mass m2 will be
in uenced by gravitational forces from both such bodies. If it is assumed that the
spacecraft is stationary and located on the vector line , then there always exists a
distance d from m3 at which the gravitational accelerations exerted by both bodies
on the spacecraft will be equal to each other. Mathematically this relationship can
be expressed by
F23 = F12.
Then from equation (2.5)
Gm2m3d2 = Gm1m2r2 : (3.1)
Now setting d = Rg1, and substituting in r = Rg1
m3( Rg1)2 = m1R2g1 (3.2)
Rg1
Rg1 = (m1=m3)
1=2
Rg1 = Rg1(m1=m3)1=2
= Rg1 1 + (m1=m3)1=2
11
which can easily be solved for Rg1, producing the simple expression for the radius of
the sphere of gravitation
Rg1 = 1 + (m
1=m3)1=2
: (3.3)
Rauschenbakh et al. [5] refer to the sphere of gravitation as the sphere of at-
traction and make the further simplifying assumption that m1 >>m3. By virtue of
this assumption, the radius of the gravitational sphere lies at a small enough distance
from m3 such that r . Under this simplifying condition, and setting d = Rg2,
equation (3.1) is reduced to
Rg2 =
rm
3
m1: (3.4)
Chebotarev [6], who de nes his sphere of gravitation identically to Rauschenbakh
et al., makes the interesting observation that of all the satellites of the major planets,
only the Earth?s moon is at all times beyond the limits of the sphere of gravitation
of the planet. This suggests that the Moon should be attracted more strongly by the
Sun than by the Earth. Indeed, if the Moon were to suddenly stop its motion, it
would fall onto the Sun, and not the Earth around which it orbits. Of course such a
scenario is not feasible, and herein lies the limitation to the signi cance of the concept
of the sphere of gravitation, since neither do bodies at absolute rest exist nor can they
be suddenly stopped.
3.2 Sphere of In uence
Pierre-Simon Laplace, an 18th century mathematician and astronomer, consid-
ered the accurate computation of the motion of a comet which approaches very near
12
to a disturbing planet. He hypothesized that \if the planet be Jupiter, its attraction
upon the comet may exceed that of the sun, and this action can entirely change the
elements of its orbit. This singular case, which appears to have taken place with the
rst comet of the year 1770, deserves particular attention" [7].
By introducing the four acceleration vectors C1, P1, C3 and P3 to respectively
represent the central and perturbing accelerations of mass m2 with respect either m1
or m3, equations (2.17) and (2.18) can be rewritten as
d2r
dt2 + C1 = P3 (3.5)
and
d2d
dt2 + C3 = P1: (3.6)
According to Laplace, the advantage of equation (3.5) over (3.6), or vice-versa,
is dependant upon the ratio of the disturbing acceleration to the central acceleration.
The idea is to select the appropriate central body so that the patched-conic approx-
imation can be applied in a reference frame in which the perturbing acceleration is
less, thereby maximizing the accuracy of the computation. The boundary de ned by
the equality of these two ratios
jP3j
jC1j =
jP1j
jC3j (3.7)
is the almost spherical surface rst described by Laplace as the sphere of in uence.
A detailed derivation given by Chobotov [8] and Battin [9] yields for the radius Ri of
the sphere of in uence, the following expression
13
Ri =
m
3
m1
2=5 1
1 + 3 cos2
1
10
(3.8)
where is the angle between the vectors d and . It can be readily seen from varying
from 0 to 90 that the surface of the sphere of in uence is a attened spheroid
with a radius variation of approximately 15% as shown in Fig. 3.1.
?8 ?6 ?4 ?2 0 2 4 6 8
x 104
?8
?6
?4
?2
0
2
4
6
8x 10
4
x (km)
y (km)
Actual surface Standard definition Laplace?s definition
Moon
Figure 3.1: Variation in the size of the Sphere of In uence.
When Laplace rst derived the radius of this sphere, he used a value of = 0 .
Ri =
m2
3
2m21
1=5
: (3.9)
14
However, it is generally agreed upon in the literature that the sphere of in uence of
the planet is de ned by the value of the maximum radius of the sphere which occurs
when = 90 .
Ri =
m
3
m1
2=5
: (3.10)
Of all the di erent de nitions of a gravitational sphere, the sphere of in uence, some-
times referred to as the sphere of activity [6] or sphere of action [8], is the most
commonly used.
3.3 Hill Sphere
The concept of the Hill sphere, named after 19th century American astronomer
George William Hill, was developed upon the principles of the restricted three-body
problem. This is a special case of the three-body problem where one of the masses is
in nitesimal, such as a spacecraft, so that its motion is a ected by the gravitational
in uences of the other two bodies, but conversely it has no appreciable a ect on either
of them. Fundamental to the solution of this problem is the expression
_x2 + _y2 + _z2 = n2(x2 +y2) + 2Gm1r + 2Gm3d C (3.11)
which is known as Jacobi?s integral. Here C is a constant of integration, n is the mean
motion pG(m1 +m3)= 3, and the left hand side is the square of the magnitude of
the velocity of the in nitesimal mass relative to the rotating reference frame.
If the velocity of the in nitesimal body goes to zero, equation (3.11) reduces to
n2(x2 +y2) + 2Gm1r + 2Gm3d = C (3.12)
15
where C is uniquely xed by the initial conditions of position and velocity. For
certain values of C, there are inaccessible regions in the rotating frame, as shown by
the shaded area in Fig. 3.2. These inaccessible regions, whose boundaries are called
zero velocity curves, divide the accessible regions into three regions, two interior and
one exterior [10].
Figure 3.2: Zero velocity curves for a speci ed value of C.
Since a particle?s velocity cannot be imaginary, it holds true that those accessible
regions, known historically as Hill?s regions, are governed by the inequality 2U > C
where
U = n
2
2 (x
2 +y2) + Gm1
r +
Gm3
d : (3.13)
16
Thus equation (3.12) is important in that it de nes for a speci c value of C
the boundaries of the regions within which the particle is con ned inde nitely. As
the value of C decreases, the area of the accessible region grows continually larger,
while the forbidden regions shrink gradually smaller until they vanish altogether. The
distance to the positions along the rotating x-axis at which the zero velocity curves
of the two interior regions rst open to each other and to the exterior region can be
approximated as
Rh =
m
3
3m1
1=3
(3.14)
These points constitute two of the three co-linear Lagrange points, and the given
approximation de nes the radius of the Hill sphere [11]. For a more comprehensive
development of the Jacobi integral and the associated zero velocity curves, the reader
is referred to Roy [12] and Battin [9].
3.4 Kislik Sphere
A nal concept for the development of a gravitational sphere was presented more
recently by Kislik in 1964 [13]. His approach is also based on the principles of the
restricted three-body problem and invokes the use of Jacobi?s integral. Using the
example of a Sun-planet system, he examined the point at which the trajectory of
an arti cial celestial body changes from rectilinear motion directed radially outward
from the planet, to elliptical motion about the sun.
If m1 and m2 are the masses of the Sun and planet respectively, the constant of
Jacobi?s integral at initial time t = t0 can be expressed in unitless quantities as
C0 = r20 + 2 r
20
+ 21 r
10
V2 r220 (3.15)
17
where r0, r10 and r20 are the initial distances of the spacecraft from the Sun-planet
barycenter, the center of the Sun and the center of the planet respectively, V is
the magnitude of the absolute velocity of the spacecraft, and is the mass ratio
m2=(m1 + m2). In the limit for r0 ! (1 ), r10 ! 1 and r20 ! 0, equation (3.15)
reduces to
C0 = 2 4 + 3 h where h = V2 r220: (3.16)
The constant of Jacobi?s integral which describes the initial part of the elliptical
trajectory at some later time t, has the similar form
C = r2 + 2 r
2
+ 21 r
1
V2 (3.17)
where V is the magnitude of the relative velocity of the spacecraft. There exists
a point along the trajectory where C = C C0 has a minimum value, and the
boundary of the gravitational sphere de ned by Kislik passes through this point. The
expression he derived for the radius Rk, of this sphere is given by
Rk = 1:15
m
3
m1
1=3
: (3.18)
Along with Kislik?s original paper, a more detailed discussion of the Kislik sphere can
be found in Rauschenbakh et al. [5].
A simple expression comparable to the Hill sphere and Kislik sphere can be made
by equating the mean motion of the in nitesimal mass m2 about the secondary body
m3, with the mean motion of m3 about the primary mass m1.
p
G(m2 +m3)=R3 =
p
G(m3 +m1)= 3: (3.19)
18
By making the simplifying assumption that m1 >>m3 >>m2, m2 can be neglected
from the left hand side of equation (3.19) and m3 can be neglected from the right
hand side. It then becomes a matter of basic algebra to obtain the expression
R =
m
3
m1
1=3
: (3.20)
In general, the de nition of a sphere of gravitational in uence of a particular
body can be viewed as the region of space within which it is possible to neglect
the gravitational in uence of third bodies and assume without great error that any
motion within it occurs in a Keplerian orbit [14]. A summarized list of the expressions
for each gravitational sphere is provided in Table (3.1) along with the values of the
radii of the various gravitational spheres for the Earth-Moon system and the Sun-
Earth system. Although there are many di erent de nitions for exactly where the
boundary of a gravitational sphere lies, the concept is certainly convenient and widely
used, particularly in planning translunar and interplanetary trajectories.
Table 3.1: Radii of various Gravitational Spheres for the Earth-Moon
(E-M) and Sun-Earth (S-E) system
Name of Sphere Expression E-M (km) S-E (km)
Sphere of Gravitation 1 Rg1 = 1+(m
1=m3)1=2
38,375 258,813
Sphere of Gravitation 2 Rg2 =
q
m3
m1 42,631 259,262
Sphere of In uence Ri =
m3
m1
2=5
66,182 924,648
Hill Sphere Rh =
m3
3m1
1=3
61,523 1,496,560
Kislik Sphere Rk = 1:15
m3
m1
1=3
102,042 2,482,175
19
Chapter 4
Chaos
4.1 Historical Perspective
In the late 1950s, a meteorologist at MIT named Edward Lorenz acquired a
Royal McBee LGP-30 computer in order to run his weather prediction simulations
which incorporated a system of twelve di erential equations to model a miniature
atmosphere. The computer was the size of a refrigerator, contained 16kB of internal
memory and it could calculate at the rate of 60 multiplications per second and print
at the speed of six lines of numbers per minute [15]. Such capabilities were extremely
impressive for the time.
In 1961 after noticing a particularly interesting sequence, Lorenz decided to re-
run the simulation and view it in more detail. However, rather than start the simu-
lation over from the beginning, he simply entered in initial values obtained from the
printout of the simulation results. Returning a short time later, he discovered that
the sequence he had observed earlier, had evolved along a completely di erent trajec-
tory. This behavior can be seen in Fig. 4.1 which illustrates a numerical simulation
of one variable in the Lorenz system. The curves represent the propagation of initial
conditions di ering by only 0.0001. At rst they appear to coincide, but soon chaotic
dynamics leads to independent, widely divergent trajectories.
After some initial investigation, Lorenz realized that the discrepancies lay in
the fact that calculations on his computer were performed with a 6-digit precision,
whereas the numbers he entered from the printout only contained three signi cant
20
Figure 4.1: Evolution of Lorenz?s experiment (Reproduced from [16]).
digits. This observation suggested that there exists, for certain types of systems, a
strong sensitivity to initial conditions. Thus, an arbitrarily small perturbation of the
trajectory at an initial time t0, may lead to signi cantly di erent behavior at some
future time.
Sensitivity to initial conditions is a characteristic of chaotic motion, and is popu-
larly known as the butter y e ect. The name comes from the title of a paper presented
by Edward Lorenz in 1972 [17] to the American Association for the Advancement of
Science in Washington, D.C. entitled \Predictability: Does the Flap of a Butter-
y?s Wings in Brazil set o a Tornado in Texas?" The paper addressed the issue of
whether the apping of a single butter y?s wing produces a tiny change in the state of
the atmosphere. Over a period of time, the tiny change can evolve into a signi cant
departure [16].
Although many others, including Hadamard, Poincar e, Birkho , Kolmogorov,
Cartwright, Littlewood and Smale, had initial insights to chaos theory prior to Lorenz,
much of it was under the auspices of ergodic theory. As such, chaos was not truly
formalized until the latter part of the 20th century, and despite the fact that Lorenz?s
1963 publication [18] was largely ignored for ten years, it is now seen as a prescient
beginning to the study of chaos.
21
4.2 Chaos Theory and Dynamic Systems
A dynamic system is a set of states combined with a rule for computing the
current values of the states from previous values. Typically, they are categorized
as either linear or nonlinear. The only possible behaviors of a linear system are (i)
unbounded growth, (ii) decay to equilibrium, or (iii) peridoic oscillation where both
(ii) and (iii) represent asymptotically periodic motion. Nonlinear systems, on the
other hand, can exhibit a broader class of behaviors including a separate type of
motion known as chaos.
Chaos is concerned with the irregular behavior of solutions to deterministic equa-
tions of motion and has received much attention from mathematicians and physicists
over recent years. The equations must be nonlinear to generate chaotic solutions, but
apart from that can be remarkably simple [19]. According to the Poincar e-Bendixson
Theorem, chaotic behavior can only arise in continuous dynamical systems whose
phase space has three or more dimensions. However the theorem does not apply
to discrete dynamical systems, where chaotic behavior can arise in two or even one
dimensional systems.
As mentioned previously, a characteristic feature of chaotic motion is that any
tiny errors in current knowledge makes accurate prediction into the future very dif-
cult to obtain. It is not however the same as stochastic motion, which occurs in
the case where even exact knowledge of the current state does not allow prediction of
future behavior. The presence of chaotic motion in models of nature can be signi -
cant, since in any physical system it is almost impossible to obtain perfectly precise
knowledge of current conditions.
The rule mentioned above is sometimes referred to as a map if it is used to denote
some function f whose domain space is equal to the range space. A particular one
22
dimensional map, called the Logistic map, is a convenient example of how a very
simple, nonlinear dynamical equation can exhibit complex and chaotic motion. It is
based on the logistic growth model equation
f(x) = ax(1 x) (4.1)
rst created by Pierre Fran cois Verhulst in 1845 to describe population growth and
decay. For the equation to possess non-trivial dynamical behavior, the growth rate
a, and the state variable x, must respectively remain on the intervals 1 4 or x exceeds unity, the system diverges towards 1 [20]. Fig. 4.2,
known as a bifurcation diagram, shows the long term possible values of x plotted
against various values of a starting at a = 2:4.
Figure 4.2: A bifurcation diagram for the Logistic map (reproduced from [21]).
Some key features to observe from this diagram include the period doubling at
approximately a = 3, the start of the period-four trajectory at a 3:45 and the
23
beginning of chaotic motion at around a = 3:57. It is also interesting to note the
windows that appear, particularly the prominent three-period window in the region
of a = 3:84.
One of the most famous maps is the 3-dimensional Lorenz attractor which shows
how the state of a dynamical system can evolve over time in a complex, non-repeating
pattern. It is governed by the ordinary di erential equations (ODE?s) used in Lorenz?s
simpli ed model of the atmosphere, and from certain perspectives, has the coinciden-
tal appearance of a butter y (see Fig. 4.3). Since the dynamics on it are chaotic, it
is referred to as a strange attractor.
Figure 4.3: The Lorenz Attractor (reproduced from [22]).
4.3 Lyapunov Exponents
Various methods have been devised over the years to make quantitative state-
ments about the exponential divergence of in nitesimally close trajectories in chaotic
24
systems. One such method, named after the Russian mathematician Aleksandr Lya-
punov, incorporates the use of special constants called Lyapunov exponents. A closely
related value is the Lyapunov number, and both can be viewed as the averaged rate
of divergence (or convergence) of nearby points along the trajectory [23].
More formally, let f be a smooth map of the real line R. The Lyapunov number
L(x1) of the trajectory fx1;x2;x3;:::g is de ned as [15]
L(x1) = limn!1(jf0(x1)j:::jf0(xn)j) 1n (4.2)
if this limit exists. The Lyapunov exponent (x1) is de ned as
(x1) = limn!1 1n[lnjf0(x1)j+:::+ lnjf0(xn)j] (4.3)
if this limit exists. Notice that = lnL, and therefore exists if and only if L exists
and is nonzero. With this de nition of the Lyapunov exponent rmly set, it is now
possible to formally de ne chaotic motion. Once again, let f be a smooth map of
the real line R. The trajectory fx1;x2;x3;:::g generated by f is chaotic if all the
following conditions are satis ed:
1) it is bounded,
2) it is not asymptotically periodic,
3) the Lyapunov exponent (x1) is greater than zero.
A simple example of calculating L(x1) and (x1) can be illustrated by examining
the following system known as a tent map
25
T(x) =
8
>><
>>:
2x; if x 1=2
2(1 x); if x 1=2.
(4.4)
Di erentiating with respect to x yields
T0(x) =
8
>><
>>:
2; if x 1=2
2; if x 1=2.
(4.5)
Then for any trajectory,
L(x1) = limn!1(2n) 1n = 2 (4.6)
regardless of initial condition. Also, (x1) = ln(2) = 0:6931.
It is relatively easy to obtain an analytical solution for the Lyapunov exponent
of the tent map. However, for most dynamical systems, especially continuous sys-
tems, the solution is not so easy to obtain analytically and must therefore be found
numerically. One such numerical method is known as Wolf?s algorithm. Consider
two trajectories initially separated by d0. After n iterations they are separated by dn.
This separation can be expressed as
dn d0e n; (4.7)
which can then be solved for as follows
26
dn
d0 e
n
ln
d
n
d0
n
lim
n!1
1
nln
d
n
d0
:
Now consider a higher dimensional system on Rm where m is the number of
elements in the state vector. For an initial condition x0, there exists a neighboring
initial condition, namely x0 + x, in any arbitrary direction which belongs to a
hypersphere about x0. The rate of separation between the two initial conditions over
time can be di erent for di erent orientations of the initial separation vector, r, thus
transforming the hypersphere into a hyperellipse as illustrated in Fig. 4.4. A system
will have as many Lyapunov exponents as the dimension of the state space, as can
be seen in the following de nition.
Figure 4.4: Evolution of an initial in nitesimal hypersphere.
27
Let f be a smooth map on Rm, and for k = 1;:::;m, let rnk be the length of the
kth longest orthogonal axis of the ellipsoid. If r = j xj, the kth Lyapunov number
of x0 is de ned by
Lk = limn!1limr!0
rn
k
r
1
n
(4.8)
if this limit exists. The kth Lyapunov exponent is k = lnLk. Also by de nition,
L1 L2 ::: Lm and 1 2 ::: m.
The largest such exponent is of most interest because it determines the pre-
dictability of a dynamical system. It measures the average rate of exponential diver-
gence of nearby trajectories of the chaotic mapping along the direction of greatest
expansion. A positive exponent will suggest that trajectories tend to separate rather
than to remain parallel or close to one another, indicative of chaotic motion [24].
For various values of , the trajectory of a dynamical system exhibits the following
behaviors.
< 0: The trajectory attracts to a stable xed point or stable periodic orbit. Negative
Lyapunov exponents are characteristic of systems that exhibit asymptotic stability.
The more negative the exponent, the greater the stability, and the system heads to-
wards its equilibrium point as quickly as possible.
= 0: A Lyapunov exponent of zero indicates that the system is in some sort of
steady state mode. Such systems exhibit Lyapunov stability.
> 0: The trajectory is unstable and chaotic. Nearby points, no matter how close,
will diverge to any arbitrary separation. The greater the value of , the more chaotic
the motion, thus the stochastic model is characterized by the case where !1.
28
4.4 State Transition Matrix
Up to this point, discrete systems have been discussed. Now, the chaotic motion
of continuous systems will be considered. As illustrated in equation (4.8), the value
of the Lyapunov number can be calculated by having an explicit knowledge of the
ratio of the axes of the hyperellipsoid to the radius of the initial hypersphere. One
method for evaluating this ratio invokes the use of the state transition matrix, which
is a matrix of partial derivatives which maps deviations in the state vector from one
time to another [25].
Consider a dynamical system governed by the rst-order di erential equation
and a given initial condition
_x = f(x;t) x(t0) = x0 (4.9)
and de ne the state transition matrix as
@x(t)@x
0
(4.10)
where x(t) is the state vector at some future time t. It is easily observed that eval-
uating at the initial time t0 yields the identity matrix. Then by taking the time
derivative of , the following sequence of expressions is obtained
29
_ = d
dt
@x(t)
@x0
= @@x
0
dx
dt
= @@x
0
(f(x;t))
= @f@x @x@x
0
_ = A where A = @f
@x (4.11)
By using the di erential equations which describe the system, A can be calculated
and used to obtain a solution for . Then, if x0 represents a small variation from
the initial condition, the deviation from the nominal trajectory at some later time,
t, can be expressed as x(t) = x0. The magnitude of this deviation can then be
calculated as j xj=p x0T T x0. Finally, by dividing through by the magnitude
of the variation in the initial condition
j xj
j x0j =
p x
0T T x0p
x0T x0 =
pD
k x0T x0p
x0T x0 =
p
Dk (4.12)
it can be seen that the square roots of the m eigenvalues, Dk, of the symmetric ma-
trix T give the ratio of the axes of the hyperellipse to the radius of the initial
hypersphere. This value of course can then be used directly in equation (4.8) to com-
pute the associated Lyapunov exponents, remembering that the exponent of greatest
interest is the largest exponent 1.
30
Chapter 5
Simulation Model
The motion of a spacecraft in cislunar space is governed primarily by the grav-
itational elds of the Earth and the Moon. The e ects of solar gravity and the per-
turbations arising from the nonspherical shape of the attracting bodies are important
in a nal analysis but are neglected in obtaining the trajectories for this simulation.
In this study the positions of the spacecraft, Earth, and Moon are described
by the relative position vectors as shown in Fig. 5.1, where r is the position of the
spacecraft relative to the Earth, d is the position of the spacecraft relative to the
Moon, and is the position of the Moon relative to the Earth. For simplifying
purposes, the dynamical model used is assumed to be coplanar. The mass of the
Earth is me, the mass of the Moon is mm, and the mass of the spacecraft is ms.
Figure 5.1: Position vectors between the vehicle, Earth, and Moon.
The translunar trajectories studied in this work can be considered solutions to
the three-body problem. As described in Chapter 2, the relative motion of the ve-
hicle can be described with respect to either the Earth or the Moon. This gives
31
rise to alternative sets of equations of motion in either an Earth-centered (EC) or
Moon-centered (MC) reference frame. The total acceleration in either representation
contains a central-body acceleration and a perturbative term. By the condition that
ms << mm << me, it is a reasonable simpli cation to neglect ms. Then by in-
troducing the standard gravitational parameters, e=Gme and m=Gmm, equations
(2.17) and (2.18) can be re-written in a slightly simpli ed form. Thus, the equation
of motion in the EC frame becomes
d2r
dt2 +
e
r3 r = m
d
d3 +
3
(5.1)
where the Earth is the central body and the Moon causes a perturbative acceleration.
The equation of motion in the MC frame is similar, however, the Moon is now treated
as the central body and the Earth causes the perturbative acceleration.
d2d
dt2 +
m
d3 d = e
r
r3
3
: (5.2)
Since the eccentricity of the Moon?s orbit is very small, e = 0:0549, the trajectory
of the Moon was modeled as a circular orbit using the following expression.
=
2
64cos(!mt)
sin(!mt)
3
75: (5.3)
Here, the value of the Moon?s angular velocity, !m = 2 =Tm, was obtained using
Kepler?s third law of planetary motion
T
m
2
2
=
3
G(me +mm) (5.4)
32
0 10 20 30 40 50 60
0
10
20
30
40
50
60
Figure 5.2: Diagram of the translunar trajectory.
to solve for Tm, and a value of = 384;400 km was used as the semi-major axis of
the Moon?s orbit.
A nominal trip time of 3.5 days was chosen for the trajectory, and an integration
time step of 20 secs was used for all simulations. The vehicle trajectory in the EC and
MC reference frames were simulated by numerical integration of equations (5.1) and
(5.2) respectively. The Moon position was updated every half time step to maintain
consistency with the the fourth-order Runge-Kutta (RK4) integration scheme (see
Appendix A). This RK4 scheme was used to perform all integrations in this study,
and simulation accuracy was improved signi cantly by updating the ephemerides at
each intermediate time step.
33
Figure 5.3: Plot of the Earth, low Earth orbit, point of v, and initial trajectory.
The starting location of the Moon was set at a distance directed out from the
Earth along the positive x-axis as shown in Fig. 5.2. Altitudes of h = 359,750 m to
360,250 m in 50 m increments were chosen for the initial low Earth orbit (LEO), and
a v of 3102.13 m/s was applied at varying angles from = 36:970 to 36:906
in increments of 0:008 from the negative y-axis. A magni ed view of the initial
part of the trajectory is illustrated in Fig. 5.3, and for convenience, table 5.1 lists the
physical parameters used in this simulation.
There are two main types of error invoked in this particular simulation, namely
discretization and round-o errors. Whereas iterative error may be signi cant in
other simulations, it was not the case in this study. In the equations of motion,
calculation of the perturbative acceleration can include larger amounts of round-o
error than the central-body acceleration. This is because, for certain con gurations
34
Table 5.1: Physical Parameters of the Simulation Model
Parameter Value Units
t0 0 secs
t 20 secs
tf 3.5 days
m 4.90266 x 1012 m3=s2
e 3.98600436 x 1014 m3=s2
re 6,372.797 km
v 3,102.13 m/s
-36.970 to -36.906 degrees
h 359.750 to 360.250 km
384,400 km
of the bodies, calculation of the perturbative acceleration can require subtraction of
similar numbers. Considering the EC representation, when the spacecraft is close to
the Earth, summing d=d3 and = 3 can require di erencing similar numbers. When
the spacecraft is close to the Moon, calculating d = r requires di erencing similar
numbers. One method of circumventing the former di culty is by de ning two new
variables, q and f
q = r (r 2 ) (5.5)
f = 3q + 3q
2 +q3
1 + (1 +q)3=2: (5.6)
Now equation (5.1) can be expressed using f as follows [9]
d2r
dt2 =
e
r3 r
m
jr j3(r +f ): (5.7)
In order to investigate the e ect of coordinate choice on round-o error, a MAT-
LAB simulation was created that allowed switching at various points between the
35
EC and MC models. This simulation was implemented with single precision compu-
tations. To estimate the amount of roundo error, the results were compared to a
\truth" simulation which was computed using double precision and the entire trajec-
tory was propagated in EC. The root-mean-square over the entire trajectory of the
norm of the position errors between the single and double precision calculations was
chosen for the error estimate.
In MATLAB, a oating-point number handled in double precision format uses
64 bits (8 bytes) of memory storage based on IEEE Standard 754. Double-precision
variables accurately represent values to approximately fteen decimal places. The
lower 32-bit single-precision variables represent data to about seven decimal places.
Chapter 4 introduced the idea that in modern dynamic-system theory, the propa-
gation of small initial errors in numerical simulations can be related to chaotic motion.
The di erences in roundo errors between the EC and MC frames can be related to
the Lyapunov exponents of the two representations. The loss of precision due to tak-
ing the di erence of similar numbers is in some sense equivalent to high sensitivity to
small perturbations, a characteristic of chaotic motion.
For the EC reference frame, the Lyapunov exponents were calculated as follows.
Since x = [r _r]T, then by equation (4.9)
_x =
2
64_r
r
3
75 = f(x;t): (5.8)
Then from equation (4.11),
A = @f@x =
2
64@_r@r @_r@_r
@ r
@r
@ r
@_r
3
75 =
2
640 I
G 0
3
75 (5.9)
where
36
G =
2
64@ rx@rx @ rx@ry
@ ry
@rx
@ ry
@ry
3
75: (5.10)
If rx and ry are given respectively as
rx = 1rx(r2
x +r2y)3=2
m
"
rx x
[(rx x)2 + (ry y)2]3=2 +
x
( 2x + 2y)3=2
#
(5.11)
ry = 1rx(r2
x +r2y)3=2
m
"
ry y
[(rx x)2 + (ry y)2]3=2 +
y
( 2x + 2y)3=2
#
(5.12)
then
G11 = ejjrjj3 + 3 er
2
x
jjrjj5
m
jjr jj3 +
3 m(rx x)2
jjr jj5 (5.13)
G12 = 3 erxryjjrjj5 + 3 m(rx x)(ry y)jjr jj5 (5.14)
G21 = G12 (5.15)
G21 = ejjrjj3 + 3 er
2
y
jjrjj5
m
jjr jj3 +
3 m(ry y)2
jjr jj5 : (5.16)
Using A to integrate , and nding the eigenvalues, D, of the matrix T , Lyapunov
exponents were computed as
= lnpD (5.17)
37
Chapter 6
Results
0 0.5 1 1.5 2 2.5 3 3.510
?8
10?6
10?4
10?2
100
102
104
106
Time (days)
Acceleration Ratio
|Pm|
|Ce| =
|Pe|
|Cm|
ratioEC
ratioMC
Figure 6.1: Ratio of perturbative to central-body accelerations.
Figure 6.1 presents the ratios of perturbative to central accelerations for EC
and MC frames. The ratios clearly illustrate that the gravitational in uence of the
Earth dominates the trajectory at rst, but as the spacecraft approaches the Moon,
the gravitational pull of the Moon gradually becomes more dominant. Upon close
examination, it is observed that the two curves intersect at t 2:88 days which
corresponds to a distance from the Moon of approximately 59,400 km.
To examine the e ect of coordinate choice on round-o error, the root-sum-
squared (RSS) position errors, , between single and double precision simulations
38
0 0.5 1 1.5 2 2.5 3 3.50
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
Time (days)
Position Error,
?, (km)
SP = 0.83tf = 2.905 days
Figure 6.2: Plot of the RSS position error.
were calculated at each time step. These values, which represent the magnitude
of the position di erences, were obtained for the entire trajectory and plotted with
respect to time. Fig. 6.2 illustrates how the position error grows over time for the
case where h = 360 km, = 36:930 and where the reference frame is switched
from EC to MC at t = 2:905 days. It is particularly noticeable that there is a distinct
change in the error behavior at the switch point time, proving that a careful selection
of the time at which to switch reference frames can a ect the round-o error and
improve the simulation accuracy.
Figure 6.3 shows sample results for the round-o error between single and dou-
ble precision calculations using several di erent switch points. These results clearly
indicate that switching coordinates either very early or very late in the trajectory
introduces large round-o errors. As mentioned previously, these errors can be miti-
gated by selecting a more optimal point at which to switch from EC to MC. For this
39
study, the switch point time was expressed as a percentage of the nal time tf. By
plotting the RSS errors for ninety nine switch points, it was observed that switching
at t = 2:905 days produces a particularly low position error. This corresponds to a
distance from the Moon of approximately 57,200 km.
0 0.5 1 1.5 2 2.5 3 3.50
0.5
1
1.5
2
2.5
3
3.5
4
Time (days)
Position Error,
?, (km)
SP = 0.83tf = 2.905 days
SP = 0.05tf = 0.175 days
SP = 0.95tf = 3.325 days
Figure 6.3: RSS position error for three switch point times.
To fully characterize this behavior a range of switch points was investigated, and
each switch point was tested for a family of initial conditions. The average error of
the plot shown in Fig. 6.2 can be quanti ed by nding the root-mean-square (RMS)
error, rms, over the entire trajectory as determined from the following expression
rms =
r
20 + 21 + 22+;:::;+ 2n
n =
vu
ut1
n
nX
i=0
2i (6.1)
where n = tf= t. The RMS errors for ninety nine switch points were then calculated
for 121 di erent sets of initial conditions. The initial conditions investigated were
40
taken from a combination of the various values for h and as described in Chapter
5. Figure 6.4 shows the maximum, minimum, and average RMS values for all initial
conditions at each switch point. Although the RMS values begin high, they quickly
approach a steady value for later switch point times. Note that the maximum RMS
value is less than 20 km.
0 20 40 60 80 1000
2
4
6
8
10
12
14
16
18
20
Switch Point (Percentage of tf)
RMS Error,
?, (km)
Figure 6.4: Plot of the Maximum, Minimum and Average RMS values
A magni ed plot of the average RMS values is shown in Fig. 6.5. Although
the lowest average RMS values are obtained be switching reference frames early in
the trajectory, certain initial conditions can still induce much higher round-o errors.
Therefore, the two points of main interest in the plot appear at switch point times of
approximately 67.5% and 85% of tf. The rst point corresponds to a time of 2.3625
days and a distance from the Moon of 102,300 km. Switching from EC to MC at
this time produces a sharp reduction in the average RMS value indicative of greater
simulation accuracy. The second later point occurs at 2.975 days and a distance from
41
0 20 40 60 80 1002.6
2.65
2.7
2.75
2.8
2.85
2.9
2.95
3
Switch Point (Percentage of tf)
RMS Error,
?, (km)
Figure 6.5: Magni ed Plot of the Average RMS values
the Moon of 51,200 km. Figure 6.5 shows a very slight increase in RMS error at this
point.
It was also necessary to determine if the results in this study were representative
of a wide range of initial conditions. To investigate this issue, a switch point time of
1.75 days (0.5tf) was chosen. The average value of the RMS error at that particular
switch point was then calculated for varying numbers of simulation trials from 1 to
121, where a simulation trial is de ned as a simulation with a speci c set of initial
conditions. Figure 6.6 shows the results of this investigation. It is evident that by
examining less than twenty sets of initial conditions, there is a lot of variance in
the results. However, after approximately eighty simulation trials, the average value
has converged to a fairly steady value. It can therefore safely be assumed that an
investigation of anymore than eighty sets of initial conditions will essentially produce
a negligible change in the results.
42
0 20 40 60 80 100 120 1401
1.5
2
2.5
3
3.5
Number of Simulation Trials
Average RMS value (km)
Figure 6.6: Plot to show if simulation results are representative of many initial con-
ditions.
Finally, Fig. 6.7 displays the Lyapunov exponents for the translunar trajectory
over one hour intervals. Since the values of the exponents are greater than zero
throughout the entire trajectory, it can be asserted that the motion is always chaotic.
It is also evident, however, that the values of are greater while the spacecraft is
in close proximity to the Earth or to the Moon. This suggests that the trajectory is
more chaotic during the initial and nal stages of the journey.
Table 6.1: Radii of various Gravitational Spheres and Corresponding Switch
Point Times
Sphere Radius (km) Switch Point (Days)
Rg1 38,375 3.12
Rg2 42,631 3.07
Rh 61,523 2.85
Ri 66,182 2.80
Rk 102,042 2.36
43
0 0.5 1 1.5 2 2.5 3 3.58
8.1
8.2
8.3
8.4
8.5
8.6
8.7
8.8
8.9
Lyapunov Exponent,
?
time (days)
earth centered
moon centered
Figure 6.7: Plot of the Lyapunov Exponents for a translunar trajectory.
A summary of the radii of the ve gravitational spheres is shown in Table 6.1
along with the corresponding switch point times. Note that the spacecraft enters the
Kislik sphere after approximately 2.36 days. This corresponds very well to the sharp
decline in the plot shown in Fig. 6.5. Both the Hill sphere and the sphere of in uence
boundaries occur before the slight rise in the RMS error occurs at 2.975 days. Also,
the range of the spacecraft from the Earth and from the Moon is shown in Fig. 6.8.
This plot indicates that the spacecraft is equidistant from the Moon and the Earth
after approximately 1.062 days.
44
0 0.5 1 1.5 2 2.5 3 3.50
0.5
1
1.5
2
2.5
3
3.5
4 x 10
5
Time (days)
Distance (km)
r
d
?
Figure 6.8: Plot of the Spacecraft Range from the Earth and Moon
45
Chapter 7
Discussion
7.1 Conclusion
In the simulation of translunar trajectories, the vehicle motion can be described
in either an EC or MC frame. These coordinate choices have signi cant impact on the
amount of round o error in the simulation. It was shown that the switching from
EC to MC coordinates at an appropriate point along the trajectory can minimize
precision loss caused by this round o error. The preferred switch point is correlated
with the historical concepts of the various gravitational spheres.
The accuracy was signi cantly improved (regardless of integration frame) if the
Moon?s ephemerides were updated at each intermediate integration time step. This is
due to the fourth-order Runge-Kutta integration scheme used as described in Chapter
5.
From the results presented in Chapter 6, the magnitude of the position di erence,
induced by round-o errors, ranged from less than 1 km to over 4 km depending of
the choice of switch time. Of course, part of the reason for these large magnitudes
was the use of single-precision calculations. On the scale of the solar system, these
error values might be considered to be very small. In terms of an accurate simulation
however, errors of more than a few meters could be viewed as signi cant. The results
obtained in this work indicated a large reduction in accuracy if coordinate switching
is used in the rst 10% of the trajectory. Fortunately, the results also indicated
46
that after this initial period, coordinate switching can be performed without large
sensitivity to the particular switch point.
Finally, it was distinctly shown that the translunar trajectories investigated in
this simulation exhibit chaotic motion. As a result, there is a very apparent sensitivity
of the trajectory to initial conditions.
7.2 Future Work
Future work may further investigate other sources of numerical error in translunar
trajectory simulations. By decreasing the time step of the integration routine and
increasing the amount of precision used in the simulation, the discretization and
round-o errors could be respectively reduced. Making both these revisions however
would lead to the obvious disadvantage of a signi cant increase in the simulation
runtime. Richardson extrapolation can be used to estimate discretization error, and
the method of nearby problems can be used to characterize the total numerical error.
It may also be insightful to obtain more data by examining a wider variety of initial
conditions. This might include varying v, as well as extending the range of h and 0
values used. A more detailed analysis of the Lyapunov exponents, and the similarities
of the simulated lunar trajectories to chaos theory could also be be undertaken.
47
Bibliography
[1] Turnball, H., The Great Mathematicians, Methuen & Co Ltd, London, 1951.
[2] Gillies, G. T. The Newtonian Gravitational Constant: Recent Measurements and
Related Studies, Rep. Prog. Phys. 60, 1997, pp. 151-225.
[3] Cvitanovi c, P. Artuso, R. Mainieri, R. Tanner, G., and Vattay, G., Chaos: Clas-
sical and Quantum, Niels Bohr Institute, Copenhagen, 2005, URL:
http://chaosbook.dk/postscript.shtml [cited 30 March 08].
[4] Gurzadyan, G. A., Theory of Interplanetary Flights, Gordon and Breach Pub-
lishers, The Netherlands, 1996.
[5] Rauschenbakh, B. V., Ovchinnikov, M. Yu., McKenna-Lawlor, S., Essential
Space ight Dynamics and Magnetospherics, Kluwer Academic Publishers, The
Netherlands, 2003.
[6] Chebotarev, G. A. Gravitational Spheres of the Major Planets, Moon and Sun,
Soviet Astronomy - AJ. Vol. 7, No. 5, 1964, pp. 618-622.
[7] Laplace, P., Celestial Mechanics, Chelsea, New York, 1966.
[8] Chobotov V., Orbital Mechanics, AIAA Education Series, Virginia, 2002.
[9] Battin, R., An Introduction to the Mathematics and Methods of Astrodynamics,
AIAA Education Series, Virginia, 1999.
[10] Ross, S. D., Scheeres, D. J., Multiple Gravity Assists, Capture, and Escape in
the Restricted Three-Body Problem, SIAM J. Applied Dynamical Systems Vol.
6, No. 3, 2007, pp. 576-596.
[11] Hamilton, D. P., Burns, J. A., Orbital Stability Zones about Asteroids, Icarus 92,
1991, pp. 118-131.
[12] Roy, A. E., Orbital Motion, Adam Hilger Ltd, Bristol, 1982.
[13] Kislik, M. D., Spheres of In uence of Large Planets and the Moon, Cosmic Re-
seach Vol. 2, Issue 6, 1964, pp. 853-858.
48
[14] Radzievskii, V. V., Gravitational Capture of Cosmic Dust by the Sun and Planets
and Evolution of the Circumterrestrial Cloud, Soviet Astronomy Vol. 11, No. 1,
1967, pp. 128-136.
[15] Alligood, K. T., Sauer, T. D., Yorke, J.A., Chaos: An Introduction to Dyanmical
Systems, Springer Verlag, New York, 1996.
[16] Stewart, I., The Mathematics of Chaos, Basil Blackwell Ltd, Oxford UK, 1989.
[17] Hilborn, R., Chaos and Non Linear Dyanmics, Oxford University Press, New
York, NY, 1994.
[18] Lorenz, E. N., Deterministic Nonperiodic Flow, Journal of the Atmostpheric
Sciences Vol. 20, 1963, pp. 130-141.
[19] Casdagli, M., Chaos and Deterministic versus Stochastic Non-Linear Modelling,
Journal of the Royal Statistical Society B, Vol. 54, No. 2, 1992, pp. 303-328.
[20] May, R. M., Simple Mathematical Models with Very Complicated Dynamics, Na-
ture, Vol. 261, 1976, pp 459-467.
[21] Nonlinear Dynamics and Chaos - A Toolbox for Complex Systems Research,
Hubler, A. W., URL: http://www.how-why.com/ph510/LogisticMap Bifurca
tionDiagram.JPG [cited 24 March 08].
[22] Wikimedia, URL: http://upload.wikimedia.org/wikipedia/commons/f/f4/
Lorenz attractor.svg [cited 25 March 08].
[23] Hedrih, K., Nonlinear Dynamics and Aleksandr Mikhailovich Lyapunov, Scien-
ti c Technical Review, Vol. 57, No. 1, 2007, pp. 3-7.
[24] Wolf, R. C., Local Lyapunov Exponents: Looking Closely at Chaos, Journal of
the Royal Statistical Society B, Vol. 54, No. 2, 1992, pp. 353-371.
[25] Tapley, B. D., Schutz, B. E., Born, G. H., Statistical Orbit Determination, Else-
vier Academic Press, London, 2004.
[26] Bate, R., Mueller D. , White J., Fundamentals of Astrodynamics, Dover Publi-
cations, New York, 1971.
[27] Chaos on the Web, Cross, M., URL: http://www.cmp.caltech.edu/ mcc/Chaos
Course/Outline.html [cited 30 March 08].
[28] Freedman, R. A., Kaufmann, W. J., Universe, W. H. Freeman, New York, 2002.
49
Appendices
50
Appendix A
MATLAB Code
A.1 Double Precision Simulation
%========================================================================================%
% PROPOGATING THE EC DOUBLE PRECISION STATES %
%========================================================================================%
clear all; close all; clc;
format long; format compact;
path(path,'nMatlab2nm files')
%================================ SETTING UP THE PROBLEM ================================%
% User Inputs
for alt = 359750:50:360250; % meters
for angle = 36.890: 0.008: 36.970; % degrees
v = 3102.13; % m/s
fprintf('alt = %6.0f , ang = % 6.3f nn', alt, angle);
t0 = 0; % secs
del t = 20; % secs
tf = 3.5*86400; % secs 86400 secs = 1 day
folder = char(['alt ' int2str(alt) ' ang ' int2str(abs(angle*1000))]);
% Constants
GMm = 4.90266e12; % m3/s2
GMe = 3.98600436e14; % m3/s2
radE = 6372797; % meters
EMdist = 384400000; % meters
omega = sqrt((GMe+GMm)/EMdist?3); % rad/sec
Tm = 2*pi/omega; % secs 27d 6h 49m 50.34879957310977s
alpha = angle*pi/180; % radians
h = radE + alt; % meters
Vorbit = sqrt(GMe/h); % m/s
V = Vorbit + v; % m/s
% Initial Conditions
rxi = h*sin(alpha); % meters
ryi = h*cos(alpha); % meters
vxi = V*cos(alpha); % m/s
vyi = V*sin(alpha); % m/s
x = [rxi;ryi;vxi;vyi];
%================================ PERFORM THE PROPOGATION ===============================%
% EARTH CENTERED FRAME
frame = 1;
51
xmt = EMdist*cos(omega*t0);
ymt = EMdist*sin(omega*t0);
fid = fopen([folder 'nDblECI ' int2str(del t) '.txt'], 'w');
f = fopen([folder 'nDblmoonECI ' int2str(del t) '.txt'], 'w');
facc = fopen([folder 'nAccECI ' int2str(del t) '.txt'], 'w');
fprintf(fid, '% 6.3f, % 6.3f, % 6.3f, % 6.3f, % 6.3f nn', tf/86400, del t, ...
alt/1000, v, angle);
fprintf(fid, '% 6.0f, % .16e, % .16e, % .16e, % .16e nn', t0, x);
fprintf(f, '% 6.0f, % .16f, % .16f nn', t0, xmt, ymt);
for t = t0+del t:del t:tf
% Use the moon position at the end of the previous interval as the
% position at the start of this interval.
xm = xmt;
ym = ymt;
% Calculate the moon position at the middle and end of this interval.
xmh = EMdist*cos(omega*(t + 0.5*del t));
ymh = EMdist*sin(omega*(t + 0.5*del t));
xmt = EMdist*cos(omega*(t + del t));
ymt = EMdist*sin(omega*(t + del t));
% Propogating the states forward in time
[x,accem] = integrate(x, del t, xm, ym, xmh, ymh, xmt, ymt, frame);
fprintf(fid, '% 6.0f, % .16e, % .16e, % .16e, % .16enn',t,x);
fprintf(f, '% 6.0f, % .16f, % .16f nn', t, xmt, ymt);
fprintf(facc, '% 6.0f, % .16f, % .16f, % .16f, % .16f, % .16f, % .16f,' ...
'% .16f, % .16f nn', t, accem);
end
fclose(fid);
fclose(f);
fclose(facc);
end
end
A.2 Single Precision Simulation
%========================================================================================%
% PROPOGATING THE EC > MC SINGLE PRECISION STATES %
%========================================================================================%
clear all; close all; clc;
format long; format compact;
path(path,'nMatlab2nm files')
% User Inputs
for alt = 359750:50:360250; % meters
fprintf('alt = % 4.4f nn', alt);
for angle = 36.890: 0.008: 36.970; % degrees
fprintf('angle = % 4.4f nn', angle);
v = 3102.13; % m/s
t0 = 0; % secs
del t = 20; % secs
tf = 3.5*86400; % secs 86400 secs = 1 day
for fSPvec = 0.01:0.01:0.99;
fprintf('fSP = % 4.4f nn', fSPvec);
fSP = SinglePrec(fSPvec, alt, angle, v, t0, del t, tf);
end
end
end
52
%================================ SETTING UP THE PROBLEM ================================%
function fSP = SinglePrec(fSP, alt, angle, v, t0, del t, tf)
% User Inputs (from SwPts.m)
alt = single(alt); % meters
angle = single(angle); % degrees (CCW from the negative y axis)
v = single( v); % m/s
t0 = single(t0); % secs
del t = single(del t); % secs
tf = single(tf); % secs 86400 secs = 1 day
fSP = single(fSP);
spt = fSP*100;
tSP = floor(fSP*tf/del t + 0.5) * del t;
folder = char(['alt ' int2str(alt) ' ang ' int2str(abs(angle*1000))]);
% Constants
GMm = single(4.90266e12); % m3/s2
GMe = single(3.98600436e14); % m3/s2
radE = single(6372797); % meters
EMdist = single(384400000); % meters
omega = single(sqrt((GMe+GMm)/EMdist?3)); % rad/sec
% Tm = 2*pi/omega = 27d 6h 49m 50.34879957310977s
alpha = single(angle*pi/180); % radians
h = single(radE + alt); % meters
Vorbit = single(sqrt(GMe/h)); % m/s
V = single(Vorbit + v); % m/s
% Initial Conditions
rxi = single(h*sin(alpha)); % meters
ryi = single( h*cos(alpha)); % meters
vxi = single(V*cos(alpha)); % m/s
vyi = single(V*sin(alpha)); % m/s
x = [rxi;ryi;vxi;vyi];
%====================== PERFORM THE PROPOGATION =========================================%
% EARTH CENTERED FRAME
frame = single(1);
xmt = EMdist*cos(omega*t0);
ymt = EMdist*sin(omega*t0);
fid = fopen([folder 'nSgl ' int2str(del t) ' ' int2str(spt) '.txt'], 'w');
f = fopen([folder 'nSglmoon.txt'], 'w');
fprintf(fid, '% 6.3f, % 6.3f, % 6.3f, % 6.3f, % 6.3f, % 6.3f, % 6.3f nn', fSP, tSP,...
tf/86400, del t, alt/1000, v, angle);
fprintf(fid, '% 6.0f, % .16e, % .16e, % .16e, % .16e nn', t0, x);
fprintf(f, '% 6.0f, % .16f, % .16f nn', t0, xmt, ymt);
for t = t0+del t:del t:tSP
% Use the moon position at the end of the previous interval as the
% position at the start of this interval.
xm = xmt;
ym = ymt;
% Calculate the moon position at the middle and end of this interval.
xmh = EMdist*cos(omega*(t + 0.5*del t));
ymh = EMdist*sin(omega*(t + 0.5*del t));
xmt = EMdist*cos(omega*(t + del t));
ymt = EMdist*sin(omega*(t + del t));
% Propogating the states forward in time
x = integrate(x, del t, xm, ym, xmh, ymh, xmt, ymt, frame);
fprintf(fid, '% 6.0f, % .16e, % .16e, % .16e, % .16enn',t,x);
fprintf(f, '% 6.0f, % .16f, % .16f nn', t, xmt, ymt);
end
53
% SWITCH POINT
% Convert state vector from EC to MC
moon = [xmt;ymt];
xmdot = EMdist*omega*sin(omega*(t+del t));
ymdot = EMdist*omega*cos(omega*(t+del t));
mndot = [xmdot;ymdot];
x = x [moon;mndot];
% MOON CENTERED FRAME
frame = single(2);
for t = tSP+del t:del t:tf
% Use the moon position at the end of the previous interval as the
% position at the start of this interval.
xm = xmt;
ym = ymt;
% Calculate the moon position at the middle and end of this interval.
xmh = EMdist*cos(omega*(t + 0.5*del t));
ymh = EMdist*sin(omega*(t + 0.5*del t));
xmt = EMdist*cos(omega*(t + del t));
ymt = EMdist*sin(omega*(t + del t));
% Compute moon state at end of interval
moon = [xmt;ymt];
mndot = EMdist*omega*[ sin(omega*(t + del t));cos(omega*(t + del t))];
% Propogating the states forward in time
x = integrate(x, del t, xm, ym, xmh, ymh, xmt, ymt, frame);
% Compute the EC state from the MC state
xeci = x+[moon;mndot];
fprintf(fid, '% 6.0f, % .16e, % .16e, % .16e, % .16enn',t,xeci);
fprintf(f, '% 6.0f, % .16f, % .16f nn', t, xmt, ymt);
end
fclose(fid);
fclose(f);
A.3 Sub-Routines
%========================================================================================%
% M FILES USED FOR TRAJECTORY PROPOGATION SUB ROUTINES %
%========================================================================================%
function [xnew,accem] = integrate(x, del t, xm, ym, xmh, ymh, xmt, ymt, frame)
% Perform a 4th Order Runge Kutta Integration
x0 = x;
x0dot = deriv(x0,xm,ym,frame);
x1 = x0 + 0.5*del t*x0dot;
x1dot = deriv(x1,xmh,ymh,frame);
x2 = x0 + 0.5*del t*x1dot;
x2dot = deriv(x2,xmh,ymh,frame);
x3 = x0 + del t*x2dot; % No 0.5 factor
x3dot = deriv(x3,xmt,ymt,frame);
x4 = x0 + del t/6.0*((x0dot+x3dot) + 2.0*(x1dot+x2dot));
xnew = x4;
% Calculate the seperate acceleration components
% of the vehicle due to the earth and moon
accem = acc(x,xm,ym);
54
%========================================================================================%
function xdot = deriv(x,xm,ym,frame)
GMm = 4.90266e12; % m3/s2
GMe = 3.98600436e14; % m3/s2
pos = x(1:2); rv = pos;
vel = x(3:4);
moon = [xm;ym]; rm = moon;
if frame == 1 % For Earth Centered Frame
acc = GMe*rv/nrm(rv)?3 GMm*((rv rm)/nrm(rv rm)?3 + rm/nrm(rm)?3);
elseif frame == 2 % For Moon Centered Frame
acc = GMm*rv/nrm(rv)?3 GMe*((rv+rm)/nrm(rv+rm)?3 rm/nrm(rm)?3);
end
xdot = [vel;acc];
%========================================================================================%
function accem = acc(x,xm,ym)
GMm = 4.90266e12; % m3/s2
GMe = 3.98600436e14; % m3/s2
pos = x(1:2); rv = pos;
moon = [xm;ym]; rm = moon;
d=rv moon;
acce EC = GMe*rv/nrm(rv)?3;
accm EC = GMm*((rv rm)/nrm(rv rm)?3 + rm/nrm(rm)?3);
% ratioEC = norm(accm EC)/norm(acce EC);
accm MC = GMm*d/nrm(d)?3;
acce MC = GMe*((pos)/nrm(pos)?3 rm/nrm(rm)?3);
% ratioMC = norm(acce MC)/norm(accm MC);
accem = [acce EC' accm EC' accm MC' acce MC'];
% ratio = [ratioEC ratioMC];
%========================================================================================%
function N = nrm(vec)
N = sqrt(sum(vec.?2));
55
Appendix B
Historical Review
B.1 Claudius Ptolemaeus
Ptolemaeus (circa 100 - 170 AD), better known as Ptolemy, was a Greek as-
tronomer, mathematician and geographer. He expanded upon the work of Hipparchus
and codi ed the geocentric model of the universe using epicycles and deferents. The
astronomical data in the library at Alexandria enabled Ptolemy to deduce the sizes
and rotation rates of the epicycles and deferents of the Sun, Moon and planets. Much
of his work is assembled in the 13 volume treatise known as the Almagest.
B.2 Nicolaus Copernicus
Copernicus (1473 - 1543), a Polish lawyer, physician, mathematician and as-
tronomer, was a strong proponent of the heliocentric model of the universe, although
he still considered all planets on circular orbits. His seminal work, De Revolutionibus
Orbium Coelestium, was not published until after his death, and remained on the
index of forbidden books from 1616 - 1758.
B.3 Galileo Galilei
Galileo (1564 - 1642), born in Pisa, Italy, was one of the rst people to use a
telescope to observe the heavens. He discovered craters on the Moon, sunspots on
the Sun, the rings of Saturn, the phases of Venus and four moons orbiting Jupiter.
His discoveries strongly suggested a heliocentric structure of the universe, which con-
tradicted the prevailing opinion of the time. As such, he spent the last eight years of
his life under house arrest for suspicion of heresy.
B.4 Johannes Kepler
Kepler (1571 - 1630) was a German mathematician and astronomer. He was a
supporter of the Copernican heliocentric model, but introduced the idea of elliptic
orbits of the planets. Using the observations of Tycho Brahe, which are credited as
56
being the most accurate of the time, Kepler was able to develop the modern laws of
planetary motion.
B.5 Isaac Newton
Newton (1642 - 1727) was born in Lincolnshire county, England and is consid-
ered by some people to have been the most in uential person in the history of science.
Among many other things, he is credited with inventing the re ecting telescope, for-
mulating an empirical law of cooling and shares credit with Gottfried Leibniz for the
development of calculus. Perhaps his must signi cant contribution to science however
include his three laws of motion and his law of universal gravitation. Although twen-
tieth century scientists found that Newton?s laws break down in certain situations
(on atomic and relativistic scales), Newtonian mechanics has become the cornerstone
of modern physical science.
B.6 Joseph Louis Lagrange
Lagrange (1736 - 1813), born in Torino, Italy, introduced key concepts and de-
veloped innovative methods in mathematical analysis, di erential equations, number
theory and classical and celestial mechanics. His 1788 treatise on analytical me-
chanics, M ecanique Analytique, o ered the most comprehensive treatment of classical
mechanics since Isaac Newton. His attempts to solve the three-body problem resulted
in the discovery of the Lagrangian points in 1772. He also predicted the existence
of (Trojan) asteroids at the triangular libration points of the Sun-Jupiter system 134
years before astronomical observations con rmed their existence.
B.7 Pierre-Simon Laplace
Laplace (1749 - 1827) was born in Normandy, France, and became a professor at
the Ecole Militaire in Paris at the age of 18. Some of his major contributions include
the formulation of Laplace?s equation, inventing the Laplace transform, introducing
the concept of the potential function and his work on the stability of the solar system.
He was also one of the rst scientists to postulate the existence of black holes, and his
seminal work M ecanique C eleste helped transform the study of classical mechanics
from a geometry based study to one based on calculus.
57
B.8 Henri Poincar e
Poincar e (1854 - 1912) was a French mathematician and physicist who established
the concept of non-integrable dynamical systems. One of Poincar e?s major works
on celestial mechanics includes Les M ethodes Nouvelles de la M ecanique C eleste, in
which he aimed to completely characterize all motions of mechanical systems. His
work helped lay the foundation of modern chaos theory, as well as the the eld of
topology. He is acknowledged as a co-discoverer, with Albert Einstein and Hendrik
Lorentz, of the special theory of relativity. There is now a mathematical institute in
Paris that is named after him.
B.9 Aleksandr Lyapunov
Lyapunov (1857 - 1918) received his early education from his father, Mikhail
Lyapunov, a well known Russian Astronomer. In 1876, he entered Saint Petersburg
State University where he published two papers on hydrostatics during his fourth
year. After successfully defending his doctoral thesis on The general problem of the
stability of motion in 1892, he was appointed as a professor at Kharkov University.
He remained there for almost a decade, and in 1901 was elected as a member of the
Russian Academy of Sciences. Much of his work focussed on the stability of systems,
and probability theory. Today his name is found in many areas of mathematics such
as Lyapunov fractal, Lyapunov function, Lyapunov equation, Lyapunov stability and
the Lyapunov exponent.
58
Index
Apollo, 1
Astrodynamics, 4
Bifurcation, 23
Butter y E ect, 21
Celestial Mechanics, 4, 7, 10
Chaos, 7, 21, 36, 43, 47
Deterministic, 22
Dynamics, 4, 22, 29
Eigenvalue, 30, 37
Ergodic Theory, 21
Error
Discretization, 1, 34, 47
Iterative, 1, 34
Round-o , 1, 34, 35, 38, 46
Gravitational Constant, 4
Hill Sphere, 15
Jacobi?s Integral, 15, 17
Kepler, Johannes, 32, 56
King Oscar II, 6
Lagrange Points, 17
Laplace, Pierre-Simon, 12, 57
Lorenz Attractor, 24
Lorenz, Edward, 20
Low Earth Orbit, 1, 34
Lyapunov
Aleksandr, 25
Exponent, 25, 28, 30, 36, 43
Lyapunov Exponent, 47
Lyapunov, Aleksandr, 58
Map, 22, 25, 28
Mean Motion, 18
Newton
Law of Gravitation, 4, 5
Law of Motion, 3, 6
Newton, Isaac, 3, 57
Patched-Conic Approximation, 10
Perturbation, 9, 21, 32, 35, 38
Poincar e, Henri, 7, 21, 58
Project Constellation, 1
Runge-Kutta, 33, 46
Sensitivity, 21, 47
Sphere
Hill, 44
Kislik, 17, 44
of In uence, 12, 44
State Transition Matrix, 29, 37
Stochastic, 22
Three-Body Problem, 7, 31
Two-Body Problem, 4, 6, 10
Verholst, Pierre Fran cois, 23
Weierstrass, Karl, 7
Zero Velocity Curves, 16
59