On the Stability of Quasi-Satellite Orbits in the Elliptic Restricted Three-Body Problem Application to the Mars-Phobos System
Francisco da Silva Pais Cabral
Dissertação para a obtenção de Grau de Mestre em Engenharia Aeroespacial
Júri Presidente:
Prof. Doutor Fernando José Parracho Lau
Orientador:
Prof. Doutor Paulo Jorge Soares Gil
Vogal:
Prof. Doutor João Manuel Gonçalves de Sousa Oliveira
Dezembro de 2011
Acknowledgments
The author wishes to acknowledge
His thesis coordinator, Prof. Paulo Gil, for, one, presenting this thesis opportunity in the author’s field of interest and, second, for the essential orientation provided to the author that made this very same thesis possible,
His professors, both in IST and TU Delft, for the acquired knowledge and transmitted passion in the most diverse fields,
His university, IST, for providing the means to pursue the author’s academic and professional goals,
His colleagues for their support and availability to discuss each others’ works,
His friends for keeping the author sane.
iii
Abstract In this thesis, the design of quasi-satellites orbits in the elliptic restricted three-body problem is addressed from a preliminary space mission design perspective. The stability of these orbits is studied by an analytical and a numerical approach and findings are applied in the study of the Mars-Phobos system. In the analytical approach, perturbation theories are applied to the solution of the unperturbed Hill’s equations, obtaining the first-order approximate averaged differential equations on the osculating elements. The stability of quasi-satellite orbits is analyzed from these equations and withdrawn conclusions are confirmed numerically. We also use the fast Lyapunov indicator, a chaos detection technique, to analyze the stability of the system. The study of fast Lyapunov indicator maps for scenarios of particular interest provides a better understanding on the characteristics of quasi-satellite orbits and their stability. Both approaches are proven to be powerful tools for space mission design.
Keywords: Quasi-Satellite Orbits, Elliptic Restricted Three-Body Problem, Stability, Perturbation Theory, Fast Lyapunov Indicator, Mars-Phobos System.
v
Resumo ´ ˆ corpos el´ıptico Nesta tese, o planeamento de quase-orbitas no contexto do problema restrito dos tres ˜ espaciais. A estabilidade destas e´ estudado de uma perspectiva do planeamento preliminar de missoes ´ ´ de uma aboradagem anal´ıtica e de uma abordagem numerica ´ orbitas e´ estudada atraves e as descober˜ aplicadas ao caso do sistema Marte-Fobos. Na abordagem anal´ıtica, teorias de perturbac¸ao ˜ tas sao ˜ aplicadas a` soluc¸ao ˜ das equac¸oes ˜ ˜ perturbadas, obtendo-se as equac¸oes ˜ sao de Hill nao diferenciais ´ ´ medias aproximadas de primeira ordem nos elementos osculadores. A estabilidade de quase-orbitas ´ destas equac¸oes ˜ e as conclusoes ˜ retiradas sao ˜ confirmadas numericamente. Use´ analisada atraves ´ o indicador rapido ´ ´ ˜ de caos, para analisar a amos tambem de Lyapunov, uma tecnica de detecc¸ao ´ ´ estabilidade do sistema. O estudo dos mapas do indicador rapido de Lyapunov para cenarios de par˜ das caracter´ısticas de quase-orbitas ´ ticular interesse vai-nos providenciar uma melhor compreensao e da sua estabilidade. Ambas as abordagens demonstram ser ferramentas poderosas no planeamento ˜ espaciais. de missoes
´ ˆ Corpos El´ıptico, Estabilidade, Teoria de Palavras Chave: Quase-Orbitas, Problema Restrito dos Tres ˜ Indicador Rapido ´ Perturbac¸ao, de Lyapunov, Sistema Marte-Fobos.
vii
Contents Acknowledgments . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
iii
Abstract . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
v
Resumo . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
vii
List of Tables . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
xi
List of Figures
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
xiii
List of Acronyms . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
xv
List of Symbols . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xix 1. Introduction
1
1.1. Motivation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1
1.2. Problem Statement . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
2
1.3. Quasi-Satellite Orbits
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
3
1.4. Stability . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
4
1.5. Bibliographic Review . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
4
1.5.1. Orbits In The Three-Body Problem . . . . . . . . . . . . . . . . . . . . . . . . . . .
5
1.5.2. Chaos Indicators . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
6
1.6. Thesis Overview . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
7
2. Dynamics
9
2.1. Classical Mechanics . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
9
2.1.1. Lagrangian Mechanics . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
10
2.1.2. Hamiltonian Mechanics . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
11
2.1.3. The Variational Equations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
12
2.2. The N-Body Problem . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
13
2.3. The Restricted Three-Body Problem . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
14
2.4. Modifications of the Restricted Three-Body Problem . . . . . . . . . . . . . . . . . . . . .
17
2.4.1. The Spatial Restricted Three-Body Problem . . . . . . . . . . . . . . . . . . . . . .
17
2.4.2. The Elliptic Restricted Three-Body Problem . . . . . . . . . . . . . . . . . . . . . .
17
2.4.3. Change Of Origin
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
19
2.5. Equations of Motion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
20
2.5.1. Hamilton’s Equations of Motion . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
20
2.5.2. An Invariant Relation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
20
ix
2.6. Chaos Indicators . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
22
2.6.1. Lyapunov Characteristic Exponents . . . . . . . . . . . . . . . . . . . . . . . . . .
23
2.6.2. Fast Lyapunov Indicator . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
24
3. QSO Solutions and Stability
27
3.1. Linearization of the Equations of Motion . . . . . . . . . . . . . . . . . . . . . . . . . . . .
27
3.2. Unperturbed Hill’s Equations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
29
3.2.1. Homogeneous Solution . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
30
3.2.2. General Solution . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
32
3.2.3. Stability Considerations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
33
3.2.4. Constants of Integration . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
35
3.2.5. Constant Transformation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
35
3.3. Influence of the Second Primary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
37
3.3.1. Region of Stability . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
37
3.3.2. Approximate Solutions in the Osculating Elements . . . . . . . . . . . . . . . . . .
41
3.4. Application to the Mars-Phobos System . . . . . . . . . . . . . . . . . . . . . . . . . . . .
47
4. Numerical Exploration of QSOs
49
4.1. Numerical Integration . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
49
4.2. Implementation Validation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
49
4.2.1. Validation Test . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
50
4.2.2. Implementation Language . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
50
4.3. Computational Parameters & Methodology . . . . . . . . . . . . . . . . . . . . . . . . . .
53
4.3.1. Integration Method and Time-Step . . . . . . . . . . . . . . . . . . . . . . . . . . .
53
4.3.2. Orbit Escape . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
56
4.4. FLI Maps . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
56
4.4.1. Planar QSOs . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
57
4.4.2. Three-Dimensional QSOs . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
65
4.4.3. Velocity Maps . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
73
5. Conclusions
77
A. Programming Code
79
Bibliography
92
x
List of Tables 1.1. Mars & Phobos Parameters . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
3
3.1. Amplitudes for Values of the Mean Motion Ratio . . . . . . . . . . . . . . . . . . . . . . . .
48
4.1. Validation Test Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
51
4.2. MatLab & C Performance Comparison I . . . . . . . . . . . . . . . . . . . . . . . . . . . .
51
4.3. MatLab & C Performance Comparison II . . . . . . . . . . . . . . . . . . . . . . . . . . . .
52
4.4. Reference Values . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
53
4.5. Parameter Computation With Chosen Time-Step . . . . . . . . . . . . . . . . . . . . . . .
55
4.6. Orbits - Mean Orbital Motion Ratio . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
60
xi
List of Figures 1.1. Hill Sphere and Region of Influence of Phobos . . . . . . . . . . . . . . . . . . . . . . . .
2
3.1. Analysis of the behavior of xnp (f ) and ynp (f ) . . . . . . . . . . . . . . . . . . . . . . . . .
33
3.2. Analysis of the behavior of xnp (f ) and ynp (f ) with C3 = −e C2
. . . . . . . . . . . . . . .
33
3.3. Parametric plot of the stable solutions . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
34
3.4. Osculating Elements . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
37
3.5. Orbital Resonance . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
48
4.1. FLI Behavior . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
52
4.2. Performance Comparison . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
54
4.3. Position and Velocity Error Analysis
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
55
4.4. FLI Map - x0 Vs. y˙ 0 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
58
4.5. FLI Map - y0 Vs. x˙ 0 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
59
4.6. Planar QSOs - Tangential Entry . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
59
4.7. 2:1 Mean Motion Orbit . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
60
4.8. FLI Map - y0 Vs. y˙ 0 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
62
4.9. FLI Map & QSO - y0 Vs. y˙ 0 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
63
4.10.FLI Map - Initial True Anomaly . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
64
4.11.FLI Map - Vertical Velocity . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
66
4.12.3D QSO . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
67
4.13.FLI Map - z0 Vs. y˙ 0 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
69
4.14.FLI Map - z0 Vs. z˙0 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
71
4.15.3D QSOs - Amplitude Ratio . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
72
4.16.3D QSOs - Large Amplitude Orbits . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
72
4.17.FLI Map - Velocities I . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
74
4.18.FLI Map - Velocities II . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
75
xiii
List of Acronyms 3BP APLE CI ER3BP FLI GALI LCE MEGNO mLCE ODE PCR3BP
Three-Body Problem Averaged Power Law Exponent Chaos Indicator Elliptic Restricted Three-Body Problem Fast Lyapunov Indicator Generalized Alignment Index Lyapunov Characteristic Exponent Mean Exponential Growth factor of Nearby Orbits maximum Lyapunov Characteristic Exponent Ordinary Differential Equation Planar Circular Restricted Three-Body Problem
PSS
Poincare´ Surface of Section
QSO
Quasi-Satellite Orbit
R3BP RKn RLI SALI SER3BP S/C
Restricted Three-Body Problem n-order Runge-Kutta Relative Lyapunov Indicator Smaller ALignment Index Spatial Elliptic Restricted Three-Body Problem Spacecraft
SSN
dynamical Spectra of Stretching Numbers
VOP
Variation Of Parameters
xv
List of Symbols
a
Orbit Semi-Major Axis
C
(Modified) Jacobi Integral
ci
Arbitrary Constants
Ci
Integration Constant
d
Distance
D
Norm of the Deviation Vector
e
Orbital Eccentricity
E
3BP Energy (as function of generalized coordinates)
E(θ, k)
Incomplete Integral of the Second Kind
E(k)
Complete Integral of the Second Kind
f
Orbit True Anomaly
F (θ, k)
Incomplete Integral of the First Kind
h
Time-Step
H
Hamiltonian
J2
Second Gravitational Moment
K(k)
Complete Integral of the First Kind
L
Lagrangian
L1,2
Lagrangian Points
m
Mass
xvii
n
Orbital Mean Motion
p
Semi-Latus Rectum
P
Orbital Period
pi
Generalized Momenta
P (λ)
Characteristic Polynomial
qi
Generalized Coordinates
r1
Distance To First Primary
r2
Distance To Second Primary
RH
Hill’s Sphere Radius
T
Kinetic Energy
tL
Lyapunov Time
u, v, w
Perturbation Coordinates
U
Potential Energy
V
Velocity
W
Wronskian Determinant
w(t)
Deviation Vector
x(t)
State Vector
x, y, z
Cartesian Coordinates
Y(t)
Fundamental Matrix
G
Gravitational Constant
α, φ, δx , δy , γ, ψ
Osculating Elements
∆V
Impulse
λi
Eigenvalues
µ
Mass Parameter
∇
Vector Differential Operator
ω
Orbital Angular Velocity
xviii
Ω
Amended Potential
Φ
Fundamental Matrix
ξ, η, ζ
Pulsating Coordinates
χ1
Maximum LCE
xix
A perfection of means, and confusion of aims, seems to be our main problem.
CHAPTER
(Albert Einstein 1879 - 1955)
1
Introduction
The objective of the present work is to study the stability of orbits around the Martian moon Phobos and to analyze the motion of such orbits. The problem is complex as Keplerian orbits are not possible due to the moon’s small mass which demands the account of Mars gravitational influence in a three-body problem. The problem is addressed with both analytical and numerical approaches to identify sufficiently stable orbits around the Martian moon. The consideration of Phobos’ orbital eccentricity increases the complexity of the problem as the system becomes non-autonomous, i.e., time-dependent.
1.1. Motivation The Martian moon Phobos is one of the most prominent candidates for future space exploration missions. The interests of a mission to this moon range from scientific to engineering purposes that make it an appealing research subject. One of the major challenges in the design of such a mission is the search for sufficiently stable orbits around the moon. A mission to Phobos would answer many unresolved questions regarding the origin of the Martian satellites, believed to be captured asteroids, as well as the formation process and origin of the Solar system. If Phobos was formed by a material within the Solar System, a sample collection from this moon would provide unprecedented information on its formation and evolution. This way, a sample return mission to the moon has been a discussed option for space exploration programs, being already part of the current Russian space program — a mission was scheduled to be launched on November, 2011 (Mukhin et al., 2000; Galimov, 2010). This mission failed to begin its interplanetary flight and is set to crash on Earth during the present month. Space exploration missions are complex, expensive, and time and mass-constrained, therefore, a space station in the Solar system to refuel the spacecraft should overcome the mass and cost related limitations. Phobos is a good candidate to host such resource space station, as the ∆V ’s required to reach Phobos are smaller than those to reach Earth’s own moon, although the travel times are much larger. This can be explained by the smaller required energy to break in order to orbit the Martian moon when compared to our moon. This resource station would open the door to new possibilities for space
1
exploration programs (Wiesel, 1993). The selection of a sufficiently stable orbit to circumnavigate Phobos is very important. Depending on the mission, the spacecraft must orbit the moon enough time to update ephemerides, perform scientific experiments, select a landing site, or prepare for a landing approach and the stability of such orbit will increase the probability of mission success (Gil and Schwartz, 2010). Some space missions, more or less successful, have already targeted Phobos. The first data on the moon was obtained by the early NASA missions to Mars, Mariner 7 & 9, and the Russian mission Mars 5. The first mission to exclusively target Phobos, however, only took place in 1988 — the Russian probes Phobos 1 & 2. This mission did not succeed as the control of both satellites was lost. The second probe, however, managed to send back to Earth a number of high-definition images before being lost. Since then, data about the Martian moons has only been collected by missions targeting Mars. However, there is a Russian sample-return mission, Phobos-Grunt, that is set to be launched November 2011 after successive delays. A successful mission to Phobos is yet to be accomplished, restating such a mission as an important milestone to be overcome (Marov et al., 2004).
1.2. Problem Statement The small mass of Phobos, when compared to Mars’ (about seven orders of magnitude larger), relinquishes any possibility of a Keplerian orbit as the region of influence of Phobos is below its own surface, and its Hill sphere, with the Lagrange points L1,2 on its surface, is just a few kilometers above its surface (see fig. 1.1) with radius (Hamilton and Burns, 1992) RH ≈ a(1 − e2 )
mP h 3(mM + mP h )
1/3 (1.1)
with a being the semimajor axis, e the orbital eccentricity of Phobos, and mP h and mM the mass of Phobos and Mars, respectively. Consequently, the problem should be treated as a three-body problem, hereafter 3BP, with Mars, Phobos, and the spacecraft (S/C) (Gil and Schwartz, 2010).
Figure 1.1. Sketch of the ellipsoid representing Phobos as viewed from the North Pole looking down on the equator. The minor axis is presented as a solid line. The interior and exterior dashed lines are, respectively, the region of influence and the Hill sphere containing the Lagrange points L1 and L2 (Gil and Schwartz, 2010).
The challenge goes further. The nonspherical gravitational potential of Mars and Phobos influences the stability of possible orbits about the moon, specially the latter at close distances to the moon. The
2
inclusion of Phobos’ orbital eccentricity also poses a challenge as the problem becomes time-dependent. The exploration of this time-dependence in moons orbiting with a small eccentricity, such as Phobos, is one of the main objectives of this work. The oblateness of Mars and Phobos’ irregular shape and rotation are, however, not considered — this work is developed under the point mass approximation for both bodies. The solar radiation pressure is also neglected. The relevant bulk and orbital parameters of both Mars and Phobos are presented in table 1.1. Mars
Phobos
Bulk Parameters Mass [kg] Equatorial Radius [km] 2
6.4185 × 1023
1.06 × 1016
3396.2
13.4 × 11.2 × 9.2
4.283 × 10
7.1 × 10−4
1.96045 × 10−3
0.100523
4
µ [km/s ] J2
Orbital Parameters 2.2792 × 108
9.3772 × 103
Sidereal Orbit Period [day]
686.980
0.31891
Sidereal Rotation Period [day]
1.02595
0.31891
Orbit Inclination [deg]
1.850
1.08
Orbit Eccentricity
0.0935
0.0151
Obliquity to orbit [deg]
25.19
Tidally locked
Semimajor Axis [km]
Table 1.1. Mars and Phobos bulk and orbital parameters. Phobos radius is represented as a set of three axes: major, median and minor. The data was collected from (NASA, 2010; Gil and Schwartz, 2010).
A specific type of orbits denominated quasi-satellite orbits, QSOs, can be found under certain initial conditions. The search of sufficiently stable QSOs by analytical and numerical approaches is the main subject of this work.
1.3. Quasi-Satellite Orbits Keplerian orbits about Phobos are not possible. As sketched in fig. 1.1, the Hill sphere, region for which a body is the main attracting force, is just above Phobos’ surface which prevents us from treating this as a two-body problem. Mars influence has to be considered. However, it is possible to find sufficiently stable orbits around Phobos in the setup of the 3BP. In the case one of the primaries has a much larger mass than the other, m1 >> m2 , a special type of orbits assumes particular interest — the so-called quasi-satellite orbits (QSOs). They are a special type of orbits as they are not closed periodic trajectories although they tend to occupy the same region in space.
3
QSOs are also known as distant retrograde orbits, DROs, or instant satellite orbits. In this thesis we chose to adopt the nomenclature quasi-satellite orbits. We would also like to emphasize that a number of articles refer to QSOs as orbits near the 1:1 resonance but here this concept is extended to all resonances. All sufficiently stable QSOs are resonant orbits or close to these. Resonant orbits are any orbit that forms with the second primary (in our case Phobos) a system that orbits the same primary (Mars) whose orbital mean motions are in the ratio of small whole numbers (Peale, 1976). As the orbital amplitude of a QSO increases, the orbit tends to a 1:1 resonance, often called quasi-synchronous orbits, whereas at close distance sufficiently stable 2:1, 3:2, and 4:3 resonances can be found (Wiesel, 1993).
1.4. Stability This thesis concerns about the stability of QSOs, but first it is necessary to define the concept of stability. In the literature there are many definitions of stability and we adopt the one that best fit to help us achieve our objective, i.e., find suitable orbits around Phobos for the purposes of mission design. One of the most used definitions states that an orbit is stable if the distance to an initially nearby orbit increases linearly with time, whereas it is chaotic if the distance to an initially nearby orbit increases exponentially with time (Meyer et al., 2009). This concept is known as exponential divergence. In this work we use chaos detection techniques to analyze the stability of QSOs that are based on this definition of stability. However, apart from the use of the chaos indicator, this definition does not suit our purposes as there are orbits that escape from Phobos and enter in a orbit around Mars that present a stable nature. The stability definition best fit to our purpose is one that can fulfill the orbit’s objective, i.e., orbit Phobos during enough time to perform any reconnaissance, or scientific activities required for the mission without colliding with the moon or escaping from orbit. This way, in our work, we define a sufficiently stable orbit as follows. A sufficiently stable orbit about Phobos is one that will orbit the moon for, at least, a period of 100 revolutions of the moon around Mars without colliding against it or get more than 1,000 km away from it. The number of 100 revolutions of Phobos about Mars, about a month, is chosen. Other studies in the literature use similar timespans (25 days in (Wiesel, 1993) and 30 days in (Gil and Schwartz, 2010)). The upper limit of 1,000 kms is chosen following (Gil and Schwartz, 2010).
1.5. Bibliographic Review A fair amount of research work has been performed and published in the topics concerning this thesis, namely, orbits in the 3BP and their stability, and chaos detection techniques.
4
1.5.1. Orbits In The Three-Body Problem The study of orbital stability in the three-body problem is one of the most researched subjects in celestial mechanics. This problem concerns the study of the stability of moons in a Sun-planet system or the stability of orbits in a planet-moon system. The systems Sun-Earth-Moon, Mars-Phobos-S/C, and Jupiter-Europa-S/C are amongst the most interesting cases for their scientific interest. One of the earliest contributions to this matter is due to George W. Hill and his studies on the lunar theory (Hill, 1878). He considered the limiting case where the mass parameter of the second primary, µ, tends to zero, and scaled his spatial dimensions by the factor µ1/3 . The equations of motion under this technique are known to this day as the Unperturbed Hill’s Equations. A variant of this method is employed in this work in the analytical approach to our problem. ´ Methodes Nouvelles de la Mecanique ` Another important contribution to this field was Poincare’s ´ ´ 1993)) and Sur le probleme ` ´ Celeste (english-translated version: (Poincare, des trois corps et les equations ´ sufaces of section (PSS) is a very used technique to analyze stability in de la dynamique. Poincare’s systems with two-degrees of freedom. Victor Szebehely’s Theory of Orbits (Szebehely, 1967) is one of the most important books on the 3BP. His book collects all the major contributions on the field. Szebehely’s own contribution is the most important to this thesis; his work on developing a coordinate system suitable for the Elliptic Restricted Three Body-Problem (ER3BP), the pulsating coordinates, is of great help as the time-dependent terms reduce to a factor on the problem’s potential. ´ Henon contributed to the research of the 3BP with his numerical study of periodic orbits and their stability in the planar circular restricted three-body problem (PCR3BP). At first he studied the case of ´ primaries with equal masses (Henon, 1965a,b), extending his search to the Hill case, µ → 0, afterwards ´ (Henon, 1969). He also expanded his study to the analysis of quasi-periodic orbits and their stability ´ (Henon, 1970). Finally, he performed vertical stability analysis for the case of primaries with equal ´ masses and for Hill’s case (Henon, 1973, 1974). Between 1974 and 1977, Benest published four articles about the existence of retrograde orbits in the PCR3BP and the effect of the mass parameter µ on their stability. Amongst another research topics, he approached the 3D stability of planar orbits and the effect of the mass parameter in this (Benest, 1974, 1975, 1976, 1977a). He also contributed to the numerical exploration of stable periodic and non-periodic orbits (Benest, 1977b). In (de Broeck, 1989) quasi-synchronous orbits — QSOs near the 1:1 resonance — are studied. This work was based in the simplified model of the planar circular restricted three-body problem and uses the unperturbed Hill’s equations as a starting point applying then a perturbation to obtain the first-order approximate equations. Theoretical conditions to obtain stability are derived and applied to the case of Mars-Phobos. These conditions are verified through numerical integration of the equations of motion. One of the methodologies used in our analytical approach is based on this work but our approach differs in the used model — we consider the eccentricity of Phobos and extend our analysis to the threedimensional case. The work developed in (Kogan, 1989) studies QSOs in the vicinity of the smaller primary separated
5
by distances greater than the radius of the Hill sphere in the context of the spatial circular restricted three-body problem. The long-term evolution of QSOs is described and stability conditions are investigated. This work also used the unperturbed Hill’s equations of motions as a starting point and described their solutions in terms of the osculating elements. The variation of these elements, under the perturbation of the second primary, is then derived through the method of variation of the arbitrary constants. The differential equations on the osculating elements are averaged and, from these, stability conditions are derived. This method is also adopted in our work but the eccentricity of the primaries’ orbits are considered. In (Wiesel, 1993) sufficiently stable QSOs are explored in Phobos and Deimos with a dynamics model that includes Mars oblateness, the moons’ orbital eccentricity, its irregular shape, and rotation. The cases of zero and nonzero eccentricity are studied and sufficiently stable resonant orbits are found to exist only in the latter. Sensitivity to poorly known system parameters is explored, and accuracy requirements for spacecraft insertion maneuvers are established. We emphasize the study of resonant orbits in this work which require a more complete dynamics model as the one used by Wiesel. The proximity to Phobos (or Deimos) of these orbits requires usually neglected perturbations as the moons’ irregular shapes and their moments of inertia to be considered. An analytical approach to QSOs in the elliptic restricted problem is presented in (Lidov and Vashkov’yak, 1993). This work obtained second-order approximate solutions but these, as stated by the authors, are complex. Their method is restricted to QSOs with small inclinations. The stability of quasi-synchronous orbits around Neptune and Uranus are explored numerically in (Wiegert et al., 2000) for timespans as long as 109 years. It is found that orbits can only exist for such long time periods if at low inclinations relative to their accompanying planet and over a restricted range of heliocentric eccentricities. In (Gil and Schwartz, 2010) a numerical study of QSOs from a preliminary mission design perspective is performed. Their stability and impact on the design of a mission to Phobos (eclipses, observation conditions, etc.) is addressed. QSOs with an inclination to Phobos’ orbital plane are also studied.
1.5.2. Chaos Indicators The so-called chaos indicators (CIs) are numerical techniques that distinguish regular from chaotic motion in a dynamical system. Their number has increased over the past years. They can be divided in two categories: techniques based on the evaluation of the deviation vector (vector containing an initial deviation from the initial position) and how it evolves, and techniques based on the analysis of the orbit itself. The first are, in most cases, derivatives of the Lyapunov Characteristic Exponents, LCEs. In chapter 2 we address their worth for our work. In 1892, Lyapunov published a work on the stability of motion. This text was originally published in Russian, meanwhile translated (Lyapunov, 1992), and still remains as one of the most important references in the subject. Several methods and techniques received the name of Lyapunov as a tribute, namely the LCEs, only developed during late 1970s following theoretical work performed in the end of
6
the 1960s (Skokos, 2010). ´ Henon & Heiles, in their pioneer work (Henon, M. and Heiles, 1964), launched the bases for the development of CIs in the subsequent years. One of the first and most used techniques was the graphical ´ Surface of Section (PSS) developed for the analysis of 2D maps and of 2D Hamilmethod of Poincare’s tonian systems. However, this method is not suitable for the analysis of systems with more degrees of freedom. This graphical technique was then extended for 3D maps and 3D Hamiltonian systems by Froeschle´ (Froeschle, 1970, 1972). By this time, the analysis of the nature of trajectories in a system was limited to systems with 3 or less degrees of freedom and, hence, there was a need for tools that could work in systems with more degrees of freedom. Theoretical background for the computation of the LCEs was developed in (Oseledec, 1968). This work provided the basis for the seminal papers by Benettin et al. (Benettin et al., 1980a,b), where a combination of the important theoretical and numerical results on LCEs can be found as well as a developed explicit method for the computation of all LCEs. Most of the CIs based on the analysis of the deviation vector derive from the LCEs. The most used are based on the evolution of the deviation vector, and one may list, based on the evaluation of the maximum LCE: the Fast Lyapunov Indicator (FLI) (Froeschle´ et al., 1997), the smaller alignment index, SALI, the Generalized Alignment Index (GALI) (Skokos et al., 2007), the Mean Exponential Growth factor ´ 2000), the Relative Lyapunov Indicator (RLI) (Sandor ´ of Nearby Orbits (MEGNO) (Cincotta and Simo, et al., 2000), and the averaged power law exponent, APLE; and based on the spectrum of LCEs: the stretching numbers, the helicity angles, the twist angles, and the study of the differences between such spectrum (Skokos, 2010). A survey of some of these methods and respective comparison can be found in (Maffione et al., 2011). In the category of indicators based on the analysis of a particular orbit, one may list the frequency map analysis of Laskar, the ”0–1” test, the method of the low-frequency spectral analysis, the ”patterns method”, the recurrence plots technique, and the information entropy index (Skokos, 2010). ´ its definition has evolved until more modern Since the first published work on the FLIs by Froeschle, definitions (for instance (Fouchard et al., 2002)). This CI has been used to study the order of periodic ´ 2001). orbits and distinguish resonant from non-resonant orbits (Lega and Froeschle, In (Villac and Aiello, 2005; Villac and Lara, 2005; Villac, 2008) the FLI maps, or stability maps, were introduced for the study of stability regions near the libration points in the restricted three-body problem. These maps are presented as design tools for preliminary trajectory design associated with planetary satellite orbiter missions and they are used to find stability regions in the Jupiter-Europa system. In (Villac and Lara, 2005) unstable periodic orbits are studied as a mean to perform thrust-free, dynamical transfers between stability regions.
1.6. Thesis Overview Mathematical and physical tools that are used in our work are reviewed in chapter 2. The 3BP is also reviewed. The Hamiltonian formulation of mechanics is introduced and applied in the context of
7
the Spatial Elliptic Restricted Three-Body Problem (SER3BP) in order to obtain the equations of motion. The variational equations in the Hamiltonian formulation are presented and the CIs of interest — the LCEs and the FLIs — are introduced at the end of the chapter. In chapter 3 we address the elliptic problem from the analytical point of view for the Hill’s limiting case, µ → 0. Perturbation theories, by two different methods, are applied to the solutions of the unperturbed equations leading to the derivation of stability conditions and approximate solutions. These are applied to the Mars-Phobos system. In the numerical approach in the fourth chapter we use the FLI to identify sufficiently stable QSOs. The study of these orbits is performed with the analysis of the so-called FLI maps where the variation of the FLI value over a set of initial conditions is addressed. This work ends with the conclusions where achievements and recommendations for future work are discussed.
8
Before we take to sea we walk on land, Before we create we must understand.
CHAPTER
(Joseph-Louis Lagrange 1736 - 1813)
2
Dynamics
In this chapter we review important mathematical and physical concepts relevant to this work. Classical mechanics are reviewed (from references (Sussman and Wisdom, 2000; Meyer et al., 2009)) as they will be required to solve the addressed problem. The Hamiltonian of our problem is derived and the equations of motions and the variational equations are obtained from the Hamiltonian. Furthermore, the theory behind CIs is addressed as well as the method for the computation of the CIs of interest.
2.1. Classical Mechanics This review of the classical mechanics is a survey from the most relevant concepts to this work from reference (Sussman and Wisdom, 2000). In the variational formulation the equations of motion are formulated in terms of the difference of the kinetic and potential energies. Neither of these energies depend on how the positions and velocities are specified; the difference is characteristic of the system as a whole. Thus, there is a liberty of choice regarding on how the system is chosen to be described, in contrast with the Newtonian formulation where there is a particle-by-particle inherent description. The variational formulation has numerous advantages over the Newtonian formulation. The equations of motion for those parameters that describe the state of the system are derived in the same way regardless of the choice of those parameters: the method of formulation does not depend on the choice of the coordinate system. If there are positional constraints among the particles of a system the Newtonian formulation requires that the forces maintaining these constraints to be considered, whereas in the variational formulation the constraints can be built into the coordinates. Furthermore, the variational formulation reveals the association of conserved laws with symmetries. Considering a mechanical system composed by point masses, extended bodies can be thought as composed by a large number of these particles with spatial relationships between them. Specifying the position of all the constituent particles of a system specifies the configuration of the system. The existence of constraints between parts of the system, such as those that determine the shape of an extended body, means that the constituent particles cannot assume all possible positions. The set of
9
all configurations of the system that can be assumed is called the configuration space of the system. The dimension of the configuration space is the smallest number of parameters that have to be given to completely specify a configuration. The dimension of the configuration space is also called the number of degrees of freedom of the system. Strictly speaking, these are not the same but for most systems they are identical. In order to talk about specific configurations a set of parameters is required to label the configurations. The parameters that are used to specify the configuration of the system are called the generalized coordinates. The number of coordinates does not need to be the same as the dimension of the configuration space, though there must be at least that many. More parameters than necessary may be defined, but then the parameters will be subject to constraints that restrict the system to possible configurations, i.e., to elements of the configuration space. Hence, sets of coordinates with the same dimension as the configuration space are easier to work with because there are no explicit constraints among the coordinates to deal with. The set of generalized coordinates are represented by the n-tuple q(t) = (q1 (t), . . . , qn (t)) and the set of the rate of change of the generalized coordinates, or generalized velocities, are repre˙ sented by the n-tuple q(t) = (q˙1 (t), . . . , q˙n (t)). Slightly different generalized coordinates are used in the Hamiltonian formulation.
2.1.1. Lagrangian Mechanics The Euler-Lagrange equations or just the Lagrange equations, use the principle of stationary action to compute the motions of mechanical systems, and to relate the variational and Newtonian formulation of mechanics. If L is a Lagrangian for a system that depends on time, coordinates, and velocities, and if q is a coordinate path then ∂L d ∂L − =0 dt ∂ q˙i ∂qi
(2.1)
with i = 1, . . . , n. These Lagrange equations form a system of second-order ordinary differential equations that must be satisfied. In order to use the Lagrange equations to compute the evolution of a mechanical system a suitable Lagrangian must be found. There is no general way to construct a Lagrangian for every system, but there is an important class of systems for which it is possible to identify Lagrangians in a straightforward way in terms of the kinetic energy, T , and the potential energy, U . The key idea is to construct a Lagrangian L such that Lagrange’s equations are Newton’s equations F~ = m~a. The reader may refer to (Sussman and Wisdom, 2000) for this deduction which leads to the Lagrangian ˙ = T (t, q, q) ˙ − U (t, q) L(t, q, q)
(2.2)
It is important to emphasize that Lagrangians are not in a one-to-one relationship with physical systems — many Lagrangians can be used to describe the same physical system. A quantity that is a function of the state of the system that is constant along a solution path is called a
10
conserved quantity, a constant of motion, or, following historical practice, an integral of motion. In fact, system symmetries are associated with integrals of motion. For instance, linear momentum is conserved in a system with translational symmetry; angular momentum is conserved if there is rotational symmetry; and energy is conserved if the system does not depend on the origin of time. A generalized coordinate that does not appear explicitly in the Lagrangian is called a cyclic coordinate. The generalized momentum component conjugate to any cyclic coordinate is a constant of motion. An example is a particle in a central force field where the Lagrangian, expressed in polar coordinates, does not depend on the polar coordinate θ and the angular momentum is conserved. Another important integral of motion is the energy. If the Lagrangian L does not depend explicitly on time, energy is conserved; systems with no explicit time-dependence are called autonomous systems.
2.1.2. Hamiltonian Mechanics The formulation of mechanics with generalized coordinates and momenta as dynamical state variables is called the Hamiltonian formulation. The Hamiltonian and the Lagrangian formulations of mechanics are equivalent, but each presents a different point of view. The Lagrangian formulation is especially useful in the initial formulation of a system; the Hamiltonian formulation is especially useful in understanding the evolution, particularly when there are symmetries and conserved quantities. For each continuous symmetry of a mechanical system there is a conserved quantity. If the generalized coordinates can be chosen to reflect a symmetry, the conjugate momentum is conserved; such conserved quantities allow the deduction of important properties of the motion. The Hamiltonian formulation is motivated by the desire to focus attention on the momenta. The momenta can be rewritten in terms of the coordinates and the velocities, so, locally, the velocities can be solved in terms of the coordinates and momenta. For a given mechanical system and given coordinates, the momenta and the velocities can be deduced from one another, thus, the dynamical state of the system can be represented in terms of the coordinates and momenta just as well as with the coordinates and the velocities. This formulation of the equations governing the evolution of the system has the advantage that if some of the momenta are conserved, the remaining equations are immediately simplified. The relation between the generalized velocities and the momenta, deduced with the so-called Legendre transformation, is pi = −
∂L . ∂ q˙i
(2.3)
The Hamiltonian formulation of dynamics provides much more than the stated goal of expressing the system derivative in terms of potentially conserved quantities. The Hamiltonian formulation is a convenient framework in which the possible motions may be placed and understood. It makes possible to see the range of stable resonance motions, and the range of states reached by chaotic trajectories, and discover other unsuspected possible motions. The Hamiltonian formulation leads to many additional insights.
11
The Hamilton’s equations are dqi ∂H = dt ∂pi
dpi ∂H =− dt ∂qi
and
(2.4)
where the momenta state variables are represented by the tuple p(t) = (p1 (t), . . . , pn (t)). The first equation in (2.4) is just a restatement of the relation between the velocities and the momenta whereas the second equation characterized the evolution of the dynamical system. The Hamiltonian is obtained from the Lagrangian by the Legendre Transformation H = pq˙ − L
(2.5)
which can be proven to be, for conservative systems, H = pq˙ − L = T + U
(2.6)
In an analogous way, conserved quantities can also be found in the Hamiltonian as they were in the Lagrangian but now with an advantage: if a coordinate does not appear in the Hamiltonian the dimension of the system of coupled equations that are remaining is reduced by two — the coordinate does not appear and the conjugate momentum is constant.
2.1.3. The Variational Equations The initial deviation vector contains the displacements of an initial state of a system. The stability of a system can be analyzed by the evolution of these initial displacements. The time evolution of orbits and deviation vectors is addressed in (Skokos, 2010) and reviewed here. The Hamilton’s equations of a system with N degrees of freedom can be written in matrix notation as " x˙ = f (x) =
∂H ∂p
∂H − ∂q
#T ≡ J∇x H(t, x)
(2.7)
where x(t) = (q1 (t), . . . , qn (t), p1 (t), . . . , pn (t)), and ∇x is the vector differential operator in respect to the state of the system x, where 0N J= −IN
IN 0N
(2.8)
with IN being the N × N identity matrix and 0N the N × N matrix with all its elements equal to zero. The solution of (2.7) is formally written as x(t) = Φt (x(0)).
(2.9)
where Φt is the fundamental matrix of solutions of (2.7) and maps the evolution of x(0) to x(t). Now, let us see how can we determine the evolution of the deviation vector. We denote by ∇x Φt the
12
linear mapping of the evolution of the deviation vector w(t) w(t) = ∇x Φt w(0)
(2.10)
where w(0) and w(t) are deviation vectors with respect to the orbit x(t) at times t = 0 and t > 0, respectively. An initial deviation vector w(0) = (δq1 (0), . . . , δqn (0), δp1 (0), . . . , δpn (0)) from an orbit x(t) evolves according to the so-called variational equations ˙ = w
∂f (x(t)) · w = J ∇2x H(t, x) · w =: A(t) · w ∂x
(2.11)
where ∇2z H(t, z) is the Hessian matrix of the Hamiltonian computed on the reference orbit x(t), i.e., ∇2x H(t, x)i,j
∂H = ∂xi ∂xj Φt (x(0))
(2.12)
i, j = 1, 2, . . . , 2N
Note that equation (2.11) represents a set of linear differential equations with respect to w, having time-dependent coefficients, since matrix A(t) depends on the particular reference orbit, which is a function of time t. The solution of (2.11) can be written as w(t) = Y(t) · w(0)
(2.13)
where Y(t) is denominated the fundamental matrix of solutions of (2.11), satisfying the following equation ∂f ˙ Y(t) = (x(t)) · Y(t) = A(t) · Y(t) , ∂x
(2.14)
Y(0) = I2N .
The analysis of the deviation vector is the basis for most CIs, taking advantage of the concept — introduced by Lyapunov — of exponential divergence, i.e., that initially nearby chaotic orbits separate roughly exponentially with time, whereas nearby regular orbits separate roughly linearly with time.
2.2. The N-Body Problem Following (Meyer et al., 2009), the N gravitating bodies problem is defined. Consider N point masses moving in a Newtonian reference system, R3 , with the only force acting on them being their mutual gravitational attraction. Let the i-th particle have position vector qi and mass mi > 0. ¨ i , is equal to Newton’s second law states that the mass times the acceleration of the i-th particle, mi q the sum of the forces acting on the particle. Additionally, Newton’s law of gravity says that the magnitude of the force on particle i due to particle j is proportional to the product of the masses and inversely 2
proportional to the square of the distance between them, Gmi mj / kqj − qi k (G is the gravitational constant). The direction of this force is along a unit vector from particle i to particle j, (qj −qi )/ kqj − qi k.
13
This way, the forces acting on particle i are
¨i = mi q
N X
Gmi mj (qj − qi ) kqj − qi k
j=1,i6=j
3
=−
∂U ∂qi
(2.15)
from which the potential energy U is derived X
U =−
1≤i0 ⇒ A A +√ A >0 6 1 2 3 2 2 α α ¯ ¯ 1 − e b2 − 4c > 0 2 2 3 µ 3 µ µ A1 A2 + √ A2 > 0 − 3 (A1 + A2 ) − 4 4− √ α ¯ α ¯6 ¯3 1 − e2 1 − e2 α
(3.51) Notice that the second condition is always satisfied since A1 , A2 > 0. The conditions are then reduced to two4
α>
µ(A1 + A2 ) 3 4 − √1−e 2
!1/3
2 2 3 µ µ 3 µ 4− √ − 3 (A1 + A2 ) − 4 A1 A2 + √ A2 > 0 α ¯ α ¯6 ¯3 1 − e2 1 − e2 α
(3.52)
For given values of e and µ a minimum value of the amplitude for which it is possible to find sufficiently stable orbits is defined by conditions (3.52). The maximum limit is defined by the validity of the approximation x, y, z > |p|. Let the unperturbed system x˙ i = fi (x, f )
(3.58)
xi = Xi (c1 , c2 , . . . , cn , f ) = Xi (c, f )
(3.59)
have general solution
with derivative
n
x˙ i =
∂Xi X ∂Xi + c˙j ∂f ∂cj j=1
(3.60)
Substituting in equation (3.57) we get n
∂Xi X ∂Xi + c˙j = fi (xi , f ) + pi (xi , f ) ∂f ∂cj j=1
(3.61)
but we know that Xi is the general solution of the unperturbed system defined in (3.58), i.e., X˙ i =
42
fi (xi , f ). This way, we have
n X ∂Xi j=1
∂cj
(3.62)
c˙j = pi (xi , f )
The variations of the osculating elements in equation (3.62) can be obtained in matrix form (Schaub and Junkins, 2003) −1
c˙ = [L]
∂R ∂c
T (3.63)
where [L] is the anti-symmetric matrix defined by the Lagrangian brackets Lij =
∂g ∂ci
T
∂h ∂cj
−
∂h ∂ci
T
∂g ∂cj
(3.64)
with Lij = −Lji , Lii = 0, and g and h defined by (3.53). R is the perturbation function written as a function of the osculating elements R=
µ 1 2 2 1 + e cos f (g1 + g2 + g32 )1/2
(3.65)
This way, the differential equations on ci are obtained
α˙ = µα
[3 + 2e cos f − cos (2(f + φ))] e sin f + 2 [2 sin (f + φ) + e sin (2f + φ)] δαx + 2 [cos (f + φ) + e cos (2f + φ)] 2r23 (1 + e cos f )(2 + e2 + 3e cos f )
δy α
4 + e [6 cos f + sin f sin (2(f + φ)) + e + e cos (2f )] + 2 [2 cos (f + φ) + e cos (2f + φ)] δαx + 2 [sin (f + φ) − e sin (2f + φ)] φ˙ = µ 2r23 (1 + e cos f )(2 + e2 + 3e cos f ) (1 + e cos f )(2 sin (f φ) + e sin φ) − e sin f δ˙x = µα r23 (2 + e2 + 3e cos f )
δx α
δ
− (1 + e cos f ) αy
(2 + e cos f )(cos (f + φ) + e cos φ) + (2 + e cos f ) δαx − e sin f δ˙y = µα r23 (1 + e cos f )(2 + e2 + 3e cos f ) γ˙ = µγ
δy α
δy α
cos (f + ψ) sin (f + ψ) r23 (1 + e cos f )
cos2 (f + ψ) ψ˙ = µ 3 r2 (1 + e cos f )
(3.66)
These equations are too complex to conclude on the stability of the system. Nevertheless, if we assume e 0
(3.70)
The solution of a non-homogeneous system of differential equations with constants coefficients, x˙ = Ax + u(f ), can be written as (Russell, 2007)
x(f ) = M (f )M (f0 )
−1
f
Z x0 +
M (f )M (s)−1 u(s) ds
(3.71)
v2 eλ2 f
(3.72)
f0
where, for a second-order system M (f ) = v1 eλ1 f
with λ1,2 and v1,2 being the eigenvalues and eigenvectors, respectively, of matrix A, and x0 the system’s initial conditions. The solution of the system is of the form
δ x0
¯˙ + ad φc + cos φ¯0 ab − φ¯˙ 2
!
√ cos ( abf ) +
r
a b
δx
=
δy
¯˙ + ad φc cos φ¯ ab − φ¯˙ 2 ! r ¯˙ √ bc + φd b ¯ = δy0 − sin φ0 cos ( abf ) − ˙ a ab − φ¯2 ˙ ¯ bc + φd + sin φ¯ ab − φ¯˙ 2
!
δy0
¯˙ bc + φd − sin φ¯0 ab − φ¯˙ 2
√ sin ( abf )
!
δx0
¯˙ + ad φc + cos φ¯0 ab − φ¯˙ 2
√ sin ( abf )
−
(3.73)
where c and d are the only factors depending on the orbital eccentricity of the second primary. The system (generalized in (3.70)) has characteristic polynomial P (λ) = λ2 + a b,
a, b > 0
(3.74)
and, as the eigenvalues λ1,2 are a pure imaginary conjugate pair, the system is stable. Solutions (3.73) are composed by two distinct periodic motions. The first has an amplitude composed by both terms √ that depend and not depend on the orbital eccentricity. This motion has frequency ω1 = a b (from the argument of the trigonometric functions) and period 2π 6π 2 α ¯3 α ¯3 P1 = √ = ≈ 37.1483 µ (5EK − 4E 2 − K 2 )1/2 µ ab
(3.75)
the other distinguishable motion has an amplitude that depends on the orbital eccentricity of the second primary and has the same frequency and period as φ¯ — ω = φ¯˙ and 2
P2 =
2π 2π 2 α ¯3 α ¯3 = ≈ 16.2992 E µ µ φ¯˙
(3.76)
45
The system (3.68) is also independent on the variable β¯ and its solution is β¯ = − arctan
r
p m−n 2 2 m − n (f + Cβ¯ ) tan m+n
(3.77)
with m = 3µE/(6π α ¯ 3 ) and n = 3µ(5E − 2K)/(6π α ¯ 3 ). The frequency and period of this oscillation are √ obtained in the same fashion, ωβ¯ = m2 − n2 and Pβ¯ = √
2π 6π 2 α ¯3 α ¯3 = = P1 ≈ 37.1483 2 2 1/2 µ µ (5EK − 4E − K ) m2 − n2
(3.78)
After substitution of β¯ the solution on γ¯ is obtained ¯ −1/2 γ¯ = Cγ¯ m − n cos (2β)
(3.79)
which oscillates two times faster than β¯ and, thus, Pγ¯ = Pβ¯ /2. The relation between the orbital mean motions of the QSO, nQSO , and of the second primary, n, in their orbits around the first primary can also be obtained. The QSO has orbital mean motion n = 1 + φ¯˙ QSO
and n is 1. µE µ nQSO = 1 + φ¯˙ = 1 + ≈ 1 + 0.385491 3 3 n πα ¯ α ¯
(3.80)
The second primary, located on the origin of the reference frame, acts as a restoring force on the the third-body. This force varies with the inverse of the distance between these two bodies. If the thirdbody gets too far from the origin, there is a chance that this restoring force will not be strong enough to maintain the third-body in orbit around the second primary. The inclination of the QSO influences the distance of the third body to the second primary as the motions in the z direction and on the x-y plane are independent. Consequently, the inclination of the QSO also influences the restoring capability of the second primary. A critical value for the ratio γ/α, for which this restoring capability vanishes, can be derived. The motion of the third-body in the x-y plane is nearly elliptic and the distance to the second primary is maximum when it crosses the y-axis in y = ±2α ± δy . The motion in the z-direction has also to be considered. The maximum distance between the second and third bodies is achieved when the third-body achieves the maximum height (in absolute value) z = ±γ in the same point that achieves the maximum distance to the second primary in the x-y plane, y = ±2α ± δy . This is the worst-case scenario for the analysis of the restoring capability of the second primary with the inclination of the QSO and it is defined by an angle β = ±π/2 (angle between the intersection line of the QSO plane with the primaries’ orbital plane and the positive x semi-axis). If we give up on the assumption that the quantity q = γ/α is small, the system composed by δ¯x and δ¯y
46
in (3.70) maintains the same form but now with
a=
4Kq − (q 4 + 3q 2 + 4)Eq µ ¯3 π(q 2 + 3)(q 2 + 4)3/2 α (3.81)
4(q 2 + 4)Eq − 4Kq µ b= ¯3 π(q 2 + 3)(q 2 + 4)1/2 α where the elliptic integrals Kq and Eq have module k =
p (q 2 + 3)/(q 2 + 4). The period of the main
motion in this case is 2π 2π 2 (q 2 + 3)(q 2 + 4) α ¯3 P =√ = 1/2 µ ab [(4Kq − (q 4 + 3q 2 + 4)Eq )(4(q 2 + 4)Eq − 4Kq )]
(3.82)
We now want to compute the critical ratio qc for which separation occurs. Analytically, the period is infinite for ejected orbits, hence, the quantity qc can be computed numerically by finding the root of the denominator in (3.82) 1/2 = 0 → qc = 0.961073 P |q=qc = ∞ → (4Kq − (qc4 + 3qc2 + 4)Eq )(4(qc2 + 4)Eq − 4Kq )
(3.83)
This conclusion is based on the averaged equations of motion where the peaks of periodic effects are neglected. Thus, it is expected that separation occurs before the value found for qc . The computed value is merely a statement that QSOs with q > qc will suffer orbit separation. The numerical exploration of the problem will tell us how accurate is this prediction.
3.4. Application to the Mars-Phobos System We now apply the results to the case of a QSO around Phobos. All the parameters computed here are for the case of initial true anomaly f0 = 0. Recall the numerical values of the relevant parameters: m1 = 6.4185 × 1023 kg m2 = 1.06 × 1016 kg m2 = 1.65148 × 10−8 µ= m + m 1 2 a = 9377.2 km e = 0.0151 P = 0.31891 days = 27553.8 s 2π n= = 2.288933 × 10−4 rad/s P
(3.84)
and the relation between the normalized amplitude, α, and the dimensional amplitude α∗ α=
α∗ (1 + e cos f ) a(1 − e2 )
(3.85)
47
which for f0 = 0 results in α∗ = α a(1 − e). We omit, for the sake of simplicity, the superscript
∗
as all
values computed in this section are dimensional. This way, the stability conditions for quasi-synchronous orbits are, from equations (3.52) α > 17.3774 km α > 29.4262 km
(3.86)
with the second condition prevailing over the first. The periods of the osculating elements are, as a function of the dimensional amplitude (in kilometers) α3 = 13.1961 α3 P1 = Pβ = 37.1483 µ a3 (1 − e)3 n α3 P2 = 16.2992 = 5.49413 α3 µ a3 (1 − e)3 n 37.1483 α3 = 6.59806 α3 Pγ = 2 µ a3 (1 − e)3 n
(3.87)
The relation between the orbital mean motions of the QSO and Phobos’ can also be studied through equation (3.80). In figure 3.5 it is presented how this ratio changes with the amplitude. Moreover, in table 3.1 the amplitudes of QSOs with ratios represented by two consecutive small integers are presented. From (Wiesel, 1993), it is known that sufficiently stable resonant orbits exist for these values of the ratio nQSO /n. Note that these non-synchronous orbits are not restricted by the conditions (3.86).
Figure 3.5. The change of the ratio nQSO /n with the amplitude is plotted. The values for which it is expected to find resonances 1:1, 2:1, 3:2, and 4:3 are marked.
nQSO /n
Amplitude [km]
2:1
17.117
3:2
21.5661
4:3
24.687
Table 3.1. Amplitudes for values of the ratio nQSO /n represented by a fraction of two small consecutive integers.
48
All truths are easy to understand once they are discovered; the point is to discover them.
CHAPTER
(Galileo Galilei 1564 - 1642)
4
Numerical Exploration of QSOs
In this chapter we integrate numerically the equations of motion and the variational equations to obtain the evolution of both the orbit and the deviation vector. The solutions obtained numerically are the basis for the study of the stability of the system with the FLI maps.
4.1. Numerical Integration Numerical integration is the numerical analysis of the solutions of ordinary differential equations (ODEs). This approach is particularly useful when the analytical solution of the set of ordinary differential equations is not possible, too complex, or just too time-consuming. As a numerical method, the technique fails to provide the same insight into a problem behavior as an analytical approach. Furthermore, it incorporates in the computed solution an error that depends on the specific numerical integration method used as well as on the computational parameters such as the time step. Runge-Kutta methods are a commonly used class of numerical integrators; they are regarded as the horsepower of celestial mechanics. In this work we are going to implement and test two variants of the Range-Kutta method: the commonly used explicit fourth-order Runge-Kutta, RK4, with four function evaluations, and an explicit eight-order Runge-Kutta method with thirteen function evaluations, adapted from (Prince and Dormand, 1981). This paper defines an adaptive method with RK7 and RK8 (seventh and eighth-order Runge-Kutta), however, we need the step size h to remain constant for the evaluation of the deviation vector evolution in the computation of the FLI value. This way, only the higher-order method is adapted to our problem.
4.2. Implementation Validation The integration methods RK4 and RK8 were first implemented in the programming language MatLab. However, the program proved to be fairly slow as the CPU times required for the computation of the FLI Maps were unaffordable. Consequently, the methods were implemented in C, a language that delivers faster CPU times when compared to MatLab. In this section, we validate both implementations with
49
the exact solution of a Keplerian orbit, i.e., a satellite orbiting a body of gravitational parameter µ (do not confuse with the mass parameter of the 3BP). After the validation of both implementations, a comparison between the performance of the RK8 method, implemented in MatLab and C, is also presented.
4.2.1. Validation Test The implementation of both Runge-Kutta methods, RK4 and RK8, are validated in C. For this purpose, we test a Keplerian orbit in the two-body problem. The equations of motion are d dt
x y z x˙ y˙ z˙
0
0 0 = −µx r3 0 0
0
0
1
0
0
0
0
0
0 µ − 3y r 0
0 0 µ − 3z r
0
0
0 0 0 1 0 0 0 0 0 0 0 0 0 1
x
y z x˙ y˙ z˙
(4.1)
with initial conditions X0 =
rp 0 0 0 Vp 0
(4.2)
and the parameters for an elliptic orbit
µ = GMM
2π P = n
µ n= 3 a
s 2 1 Vp = µ − rp a
s 2 1 Va = µ − ra a (4.3)
with rp = a(1−e), ra = a(1+e), and e = 0.5. Mars is used as the attracting body and the semi-major axis a is the semi-major axis of Phobos’ orbit. The equations are integrated over 100.5 periods with initial conditions x0 = rp and y˙ 0 = Vp . It is expected that the satellite meet the final conditions x0 = −ra and y˙ 0 = −Va . The final positions and velocities normalized with the apogee radius and velocity, respectively, are presented for both methods and for different time-steps in table 4.1. The results for both methods strongly suggest that the error decreases and the CPU time increases with the decrease of the size of the time-step. Furthermore, the implementations of both methods are validated as, with the right step-sizes, the results meet the final state of the system, known beforehand.
4.2.2. Implementation Language The implementations of the RK8 method in MatLab and C are tested and compared. The implementation in C is found to be much faster in the computation of the FLI. As an interpreting language MatLab does not possess the same computational speed as the more basic language C. The comparison of the two implementations is carried out with a sample orbit in our problem. The
50
RK4
RK8
Step
P/20
P/200
P/500
P/20
P/200
P/500
xf [ra ]
160.665957
-0.999564
-0.999999
-0.814849
-1.000000
-1.000000
yf [ra ]
-507.420687
0.017827
-0.000196
-0.395417
0.000000
0.000000
zf [ra ]
0.000000
0.000000
0.000000
0.000000
0.000000
0.000000
x˙ f [Va ]
0.704791
0.035159
0.000379
0.874789
0.000000
0.000000
y˙ f [Va ]
-2.222066
-0.999795
-1.000001
-0.802449
-1.000000
-1.000000
z˙f [Va ]
0.000000
0.000000
0.000000
0.000000
0.000000
0.000000
CPU [s]
0.000
0.05
0.11
0.02
0.22
0.55
Table 4.1. Validation test results for a Keplerian orbit using RK4 and RK8 with different time-steps.
state of the system, xf , the fundamental matrix of the deviation vector evolution at the final time, Y(ff ) = ∇x Φff (defined in (2.13)), the FLI, and the computational times are compared at the end of the integration period of 100 revolutions of Phobos about Mars. These results were obtained with a small time-step of π/125 and are presented in tables 4.2 and 4.3. MatLab
C
xf [km]
24.8497
24.8496
yf [km]
152.1574
152.1578
zf [km]
0.0000
0.0000
x˙ f [km/s]
-0.0157
-0.0157
y˙ f [km/s]
-0.0057
-0.0057
z˙f [km/s]
0.0000
0.0000
dmin [km]
74.6405
74.6138
dmax [km]
249.3773
249.3668
FLI
0.1213
0.1237
CPU [s]
86.7898
4.64
Table 4.2. Comparison of the final position, FLI, minimum and maximum distance and computational time of a sample orbit with y0 = −100 km and x˙ 0 = −20 m/s with MatLab and C implementations of RK8.
The comparison of the results obtained using the two different implementations indicate that there is a slight difference on the values of the fundamental matrix of the evolution of the deviation vector that results in a small difference in the FLI value. The final result of the FLI is not important, only its behavior for different orbits, i.e., its capacity to distinguish chaosticity. Further analysis has to be performed to assess if both implementations demonstrate the same capacity to characterize the stability of orbits despite the difference in the absolute value of the FLI. This assessment is performed by analyzing the FLI values for different orbits along the y-axis from -20 km to -300 for both implementations. The results are shown on figure 4.1 where a time-step of
51
Y(ff ) [×1059 ]
Language
MatLab
C
0.3835
0.0104
0
0.0387
0.3597
0
-3.5181
-0.0955
0
-0.3551
-3.2998
0
0
0
0.0000
0
0
0.0000
3.6375
0.0987
0
0.3671
3.4117
0
-0.2385
-0.0065
0
-0.0241
-0.2237
0
0
0
0.0000
0
0
0.000
0.3631
0.0121
0
0.0387
0.3402
0
-3.4673
-0.1156
0
-0.3695
-3.2491
0
0
0
0.0000
0
0
0.0000
3.5976
0.1199
0
0.3833
3.3711
0
-0.2169
-0.0072
0
-0.0231
-0.2033
0
0
0
0.0000
0
0
0.000
Table 4.3. Comparison of the final fundamental matrix of the evolution of the deviation vector of a sample orbit with y0 = −100 km and x˙ 0 = −20 m/s with MatLab and C implementations of RK8.
h = π/125 is used. Orbits that escape from Phobos’ vicinity (more than 1000 km) are assigned a value between 1 and 2, depending on how long they are able to stay near the moon. The reasoning behind this assignment is discussed later.
(a) Full View
(b) Zoom-In
Figure 4.1. FLI Behavior computed by the method implemented in MatLab and in C.
Although there is a small difference in the absolute value of the FLI, the results presented in figure 4.1 suggest that the FLI keeps its behavior in both implementations. There is a small bump in the FLI behavior in the C implementation but this can be neglected since it does not affect its capability to identify the chaotic nature of orbits. It is important to refer that the data for this graph took more than 118 minutes to compute in the MatLab implementation against less than 2 minutes in the C implementations illustrating the reason why we use the latter.
52
4.3. Computational Parameters & Methodology We now choose the method to use, RK4 or RK8, the time-step, and discuss how the FLI value for orbits that escape Phobos (that do not meet the stability requirements defined in chapter 1) is defined. The goal is to choose a method with a proper time-step that allows the smallest CPU times possible meeting defined performance thresholds.
4.3.1. Integration Method and Time-Step In order to analyze the performance of both methods for different time-steps we set reference values (table 4.4) computed with the RK8 method with a very small step-size to serve as a basis of comparison. There is not a reference value for the FLI since its value depends on the time-step size and, hence, it is not possible to compare it for different time-step sizes. Method
RK8
Time-Step
π/2000
y0
-100 km
x˙ 0
-20 m/s
xf
24.849701 km
yf
152.157434 km
zf
0.0 km
x˙ f
-15.745171 m/s
y˙ f
-5.729053 m/s
z˙f
0.0 m/s
dmin
74.639542 km
dmax
249.380685 km
Table 4.4. Parameters computed for an orbit with initial conditions y0 = −100 km and x˙ 0 = −20 m/s with a time step of π/2000.
The analysis of the performance of both methods, and its variation with the time-step size, is carried out by analyzing the norms of the error vectors and the differences in the minimum and maximum distances. Th error vector is defined as the difference vector between the reference vectors and the obtained vectors for each step-size. The FLI values and the CPU times for different time-steps are also computed for both methods. The results are presented in figures 4.2. The results suggest that the RK8 outperforms the RK4 in every parameter with the exception of the CPU time. However, the RK8 method converges faster so we define performance thresholds to assess which of the methods has smaller CPU times meeting these thresholds. Note that the analysis of the error in the computation of the minimum and maximum distances is a little more complex — with lesser steps, less points are analyzed and there is a higher probability of missing the points where the minimum and maximum distances are achieved by a fair distance. The result
53
(a) Position Error
(b) Velocity Error
(c) Minimum Distance Error
(d) Maximum Distance Error
(e) FLI
(f) CPU Time
Figure 4.2. Performance comparison between RK4 and RK8 for different time-steps. The errors are computed in relation to the reference values previously computed. Note that the performance values are evaluated as functions of the number of integration steps per period 2π and not directly as functions of the time-step.
is the observed periodicities of these errors but these can not be minimized since different orbits will present different periodicities depending on the location of the points where the minimum and maximum distances to the second primary are achieved. These parameters are not critical as we only use them to assess if an orbit collides against the moon (d < 15 km) or if it escapes from the QSO (d > 1000) km and not to perform any analysis on its stability. Nevertheless, these errors are fairly small and the minimum and maximum distance achieved can be used for orbit design purposes. We define the thresholds as a maximum position error of 1 m and a maximum velocity error of 0.001 m/s. Figures 4.2(a) and 4.2(b) present a plot of the error in these parameters with the number of steps per period for both methods. The RK4 needs almost 200 steps per period to meet the same performance than the RK8 with 20 steps per period. Figure 4.2(f) suggests that for these numbers of steps per period
54
the RK8 method presents smaller CPU times and, consequently, it is chosen as the integration method to be used for the computation of the FLI maps.
(a) Position Error (Zoom-In)
(b) Velocity Error (Zoom-In)
Figure 4.3. The zoom-in of both the position and velocity errors allows us to choose a proper step-size to meet the defined thresholds.
Figure 4.3 presents a zoom-in of the errors in the position and velocity for a more careful choice of the number of steps per period. The results suggest that 16 steps per period is enough for the RK8 to meet the defined thresholds. This way, we use a time-step h = π/8. In table 4.5 the used reference values and the computed values with the selected time-step are presented. Remember that the FLI values are not comparable when computed with different time-steps.
Method Time-Step
RK8 π/2000
π/8
y0 [km]
-100
x˙ 0 [m/s]
-20
xf [km]
24.849701
24.849769
yf [km]
152.157434
152.157531
zf [km]
0.0
0.0
x˙ f [m/s]
-15.745171
-15.745201
y˙ f [m/s]
-5.729053
-5.729068
z˙f [m/s]
0.0
0.0
dmin [km]
74.639542
74.658947
dmax [km]
249.380685
248.344936
0.007763
1.447331
73.57
0.29
FLI CPU Time [s]
Table 4.5. Parameter computation with chosen time-step and respective comparison with the reference values.
55
4.3.2. Orbit Escape The two objectives of the numerical exploration are to obtain the solution of orbits and to assess which of these are sufficiently stable by the study of the deviation vector. A probe that is ejected from orbit is not necessarily unstable (from the point of view of the concept of exponential divergence, measured by the FLI). In fact, when a probe escapes from Phobos’ vicinity it can enter in a Keplerian orbit around Mars which is more stable than a QSO. In this sense, the FLI is not suited to distinguish orbits that escape from Phobos’ vicinity and are not stable following our criteria defined in chapter 1. With the purpose of separate the ejected orbits from the remainder, we will define the FLI of an orbit that escaped Phobos in our algorithm as F LIesc = 20 − 10
fesc ff
(4.4)
where fesc is the normalized time when the third body escaped from Phobos and ff is the total time of integration, i.e., 100 revolutions of Phobos around Mars. This way, the range of FLI values 10–20 is reserved for the assessment of ejected orbits and in this case the FLI measures the time that a probe remained in Phobos’ vicinity without colliding against it. This provides the information if an orbit can be used for a smaller timescale or if it escapes or crashed Phobos too soon to be used for any mission purposes.
4.4. FLI Maps The representation of the values of the FLI over a set of initial conditions generates a FLI map (or stability map) where regions filled with sufficiently stable orbits can be identified to provide insight into the characteristics of these regions. We study the problem of the ER3BP which is defined by 6 parameters plus the origin of the normalized time f0 . The FLI maps represent the FLI value as a function of two varying parameters, hence, the set of initial conditions needs to have 5 fixed initial parameters. A complete study of the stability of the system is not possible but it is possible to study the stability in a plane of two varying parameters. This approach provides enough insight into the stability properties of the system to help in the choice of an entry point, i.e., the point at which a probe is maneuvered to enter in a sufficiently stable QSO. Following (Villac and Aiello, 2005), we introduce the modified Jacobi integral1 defined in 2.5.2 C=2
Ω0 −V2 1 + e cos f
(4.5)
where V = (x˙ 2 + y˙ 2 + z˙ 2 )1/2 is the spacecraft velocity and Ω0 is the amended potential Ω0 = 1 Strictly
1 1−µ µ 1 ((x + µ − 1)2 + y 2 − e z 2 cos f ) + + + µ(1 − µ) 2 r1 r2 2
(4.6)
speaking, this modification of the Jacobi integral is not an integral of motion in the elliptic case but, as it possesses an invariant relation (which is constant for a constant value of f ), we call it modified Jacobi integral for the sake of simplicity.
56
with r1 and r2 being the distance of the spacecraft to Mars and Phobos, respectively. The greatest influence in the amended potential Ω0 is Mars gravitational potential and this is almost constant in the vicinity of Phobos which causes the modified Jacobi integral variation to be predominated by the variation of the spacecraft’s velocity and by the variation of the true anomaly f . The characterization of FLI maps as a function of the modified Jacobi integral provides insight into the sensibility of the stability regions to the initial conditions. Orbits are integrated over 100 periods of Phobos around Mars and the FLI value is computed in respect to this time-span. The dark regions of the FLI maps correspond to regions of mostly regular motion whereas the near white regions correspond to orbits that escaped or collided in the moon. The orange regions correspond to orbits that were able to remain stable for some time before escaping or crashing the moon. We study the particular case of quasi-synchronous orbits. Resonant orbits are not considered since in our model of the dynamics we did not considered the perturbations that lock QSOs at resonances represented by two small whole integers, namely, Phobos’ irregular shape and rotation. Nevertheless, the relation between the mean motions of the QSO and Phobos is addressed.
4.4.1. Planar QSOs We begin with the simpler case of planar QSOs before extending the study to three-dimensional orbits.
x0 Vs. y˙ 0 Let us begin with the stability analysis of planar QSOs with the variation of x0 and y˙ 0 . This analysis is performed for the case of initial true anomaly f0 = 0, i.e., Phobos’ perigee passage. The other initial parameters y0 , z0 , x˙ 0 , z˙0 are set to zero. The orbital amplitude in the x-y plane, α, is best measured in the passage of the orbit by y0 = 0 because experimentation suggests that the displacement of the QSO in the x direction is much smaller than in the y direction, i.e., δx FLI = log ( alpha [ i ] ) ; } } printf ( ” \ n O r b i t Frequency : %L f \n ” , freq / 1 0 0 ) ; }
Num Exp.c
87
Bibliography Arnold, V. I., Weinstein, A. and Vogtmann, K.: 1978, Mathematical Methods of Classical Mechanics, Springer. Benest, D.: 1974, Effects of the Mass Ratio on the Existence of Retrograde Satellites in the Circular Plane Restricted Problem, Astronomy and Astrophysics 32, 39–46. Benest, D.: 1975, Effects of the mass ratio on the existence of retrograde satellites in the circular plane restricted problem. II, Astronomy and Astrophysics 45, 353–363. Benest, D.: 1976, Effects of the mass ratio on the existence of retrograde satellites in the circular plane restricted problem. III, Astronomy and Astrophysics 53, 231–236. Benest, D.: 1977a, Effects of the mass ratio on the existence of retrograde satellites in the circular restricted problem. IV - Three-dimensional stability of plane periodic orbits, Astronomy and Astrophysics 54, 563–568. Benest, D.: 1977b, Stable Orbits in the Circular Plane Restricted Three-Body Problem, Rev. Mexicana Astron. Astrofis. 3, 151–158. Benettin, G., Galgani, L., Giorgilli, A. and Strelcyn, J.-M.: 1980a, Lyapunov characteristic exponents for smooth dynamical systems and for hamiltonian systems; a method for computing all of them. part 1: Theory, Meccanica 15, 9–20. DOI: 10.1007/BF02128236. Benettin, G., Galgani, L., Giorgilli, A. and Strelcyn, J.-M.: 1980b, Lyapunov characteristic exponents for smooth dynamical systems and for hamiltonian systems; a method for computing all of them. part 2: Numerical application, Meccanica 15, 21–30. DOI: 10.1007/BF02128237. ´ C.: 2000, Simple tools to study global dynamics in non-axisymmetric galactic Cincotta, P. M. and Simo, potentials - I, Astronomy and Astrophysics, Supplement 147, 205–228. Danby, J. M. A.: 1962, Fundamentals of Celestial Mechanics, Willmann-Bell, Inc. de Broeck, P.: 1989, Stable relative motion as a solution of the restricted three-body problem, Working paper 399, European Space Operations Centre, ESA Orbit and Attitude Division, Darmstadt, Germany.
89
Efroimsky, M.: 2002, Preprint # 1844 of the ima equations for the keplerian elements: Hidden symmetry. ´ ´ ´ ´ Floquet, G.: 1883, Sur les equations differentielles lineaires a` coefficients periodiques, Ann. Sci. Ecole Norm. Sup. 12, 47+. ´ C. and Froeschle, ´ C.: 2002, On the Relationship Between Fast Fouchard, M., Lega, E., Froeschle, Lyapunov Indicator and Periodic Orbits for Continuous Flows, Celestial Mechanics and Dynamical Astronomy 83, 205–222. Froeschle, C.: 1970, Numerical Study of Dynamical Systems with Three Degrees of Freedom. I. Graphical Displays of Four-Dimensional Sections, Astronomy and Atrophysics 4, 115–128. Froeschle, C.: 1972, Numerical Study of a Four-Dimensional Mapping, Astronomy and Astrophysics 16, 172–189. ´ C. and Lega, E.: 2000, On the Structure of Symplectic Mappings. The Fast Lyapunov IndicaFroeschle, tor: a Very Sensitive Tool, Celestial Mechanics and Dynamical Astronomy 78, 167–195. ´ C., Lega, E. and Gonczi, R.: 1997, Fast Lyapunov Indicators. Application to Asteroidal Motion, Froeschle, Celestial Mechanics and Dynamical Astronomy 67, 41–62. Galimov, E. M.: 2010, Phobos sample return mission: Scientific substantiation, Solar System Research 44(1), 5–14. DOI: 10.1134/S0038094610010028. Gil, P. J. S. and Schwartz, J.: 2010, Simulations of quasi-satellite orbits around phobos, Journal of Guidance, Control, and Dynamics 33(3), 901–914. DOI: 10.2514/1.44434. Gradshteyn, I. S. and Ryzhik, I. M.: 1994, Table of Integrals, Series, and Products, Fifth Edition, 7th edn, Academic Press. Hamilton, D. P. and Burns, J. A.: 1992, Orbital stability zones about asteroids. II - The destabilizing effects of eccentric orbits and of solar radiation, Icarus 96, 43–64. ´ ´ ` ´ ´ Henon, M.: 1965a, Exploration numerique du probleme restreint. I. Masses egales ; orbites periodiques, Annales d’Astrophysique 28, 499–511. ´ ´ ` ´ Henon, M.: 1965b, Exploration numerique du probleme restreint. II. Masses egales, stabilite´ des orbites ´ periodiques, Annales d’Astrophysique 28, 992–1007. ´ Henon, M.: 1969, Numerical exploration of the restricted problem, V, Astronomy and Astrophysics 1, 223–238. ´ Henon, M.: 1970, Numerical exploration of the restricted problem. VI. Hill’s case: Non-periodic orbits., Astronomy and Astrophysics 9, 24–36. ´ Henon, M.: 1973, Vertical Stability of Periodic Orbits in the Restricted Problem, Celestial Mechanics 8, 269–272.
90
´ Henon, M.: 1974, Vertical Stability of Periodic Orbits in the Restricted Problem. II. Hill’s Case, Astronomy and Astrophysics 30, 317–321. ´ Henon, M. and Heiles, C.: 1964, The applicability of the third integral of motion: Some numerical experiments, The Astronomical Journal 69, 73–79. Hill, G. W.: 1878, Researches in the lunar theory, American Journal of Mathematics 1(1–3), 5–26, 129– 147, 245–260. Kogan, A. L.: 1989, Distant satellite orbits in the restricted circular three-body problem, Cosmic Research (Translation of Kosmicheskie Issledovaniya) 26, 705–710. ´ C.: 2001, on the relationship between fast lyapunov indicator and periodic orbits Lega, E. and Froeschle, for symplectic mappings, Celestial Mechanics and Dynamical Astronomy 81, 129–147. Lidov, M. L. and Vashkov’yak, M. A.: 1993, Theory of perturbations and analysis of the evolution of quasi-satellite orbits in the restricted three-body problem, Kosmicheskie Issledovaniia 31, 75–99. Lyapunov, A.: 1992, The general problem of the stability of motion, Control Theory and Applications Series, internacional journal of control edn, Tayor & Francis. Maffione, N. P., Darriba, L. A., Cincotta, P. M. and Giordano, C. M.: 2011, A comparison of different indicators of chaos based on the deviation vectors: application to symplectic mappings, Celestial Mechanics and Dynamical Astronomy 111, 285–307. DOI: 10.1007/s10569-011-9373-z. Marov, M. et al.: 2004, Phobos-grunt: Russian sample return mission, Advances in Space Research 33, 2276–2280. DOI: 10.1016/S0273-1177(03)00515-5. Meyer, R. M., Hall, G. R. and Offin, D.: 2009, Introduction to Hamiltonian Dynamical Systems and the N-Body Problem, Springer, New York, NY. DOI: 10.1007/978-0-387-09724-4. Mukhin, L., Sagdeev, R., Karavasili, K. and Zakharov, A.: 2000, Phobos, deimos mission, in N. G. Barlow (ed.), Concepts and Approaches for Mars Exploration, Lunar and Planetary Institute, pp. 230–232. NASA: 2010, Mars fact sheet, http://nssdc.gsfc.nasa.gov/planetary/factsheet/marsfact.html. Oseledec, V. I.: 1968, A multiplicative ergodic theorem. lyapunov characteristic numbers for dynamical systems., Trans. Moscow Math. Soc. 19, 197–231. Peale, S. J.: 1976, Orbital Resonances In The Solar System, Annual Review of Astronomy and Astrophysics 14, 215–246. ´ H.: 1993, New Methods of Celestial Mechanics, American Institute of Physics. Translation of: Poincare, Les methodes nouvelles de la mecanique celeste. Polyanin, A. D. and Zaitsev, V. F.: 1995, Handbook of exact solutions for ordinary differential equations, CRC Press, Boca Raton, FL.
91
Prince, P. J. and Dormand, J. R.: 1981, High order embedded Runge-Kutta formulae, Journal of Computational and Applied Mathematics 7(1), 67–75. Russell, D. L.: 2007, Inhomogeneous linear systems, Lecture Notes. Course Math 2214, Virginia Tech. ´ ´ Sandor, Z., Erdi, B. and Efthymiopoulos, C.: 2000, The Phase Space Structure Around L4 in the Restricted Three-Body Problem, Celestial Mechanics and Dynamical Astronomy 78, 113–123. Schaub, H. and Junkins, J. L.: 2003, Analytical mechanics of space systems, AIAA education series, American Institute of Aeronautics and Astronautics. Skokos, C.: 2010, The Lyapunov Characteristic Exponents and Their Computation, in J. Souchay & R. Dvorak (ed.), Lecture Notes in Physics, Berlin Springer Verlag, Vol. 790 of Lecture Notes in Physics, Berlin Springer Verlag, pp. 63–135. Skokos, C., Bountis, T. C. and Antonopoulos, C.: 2007, Geometrical properties of local dynamics in Hamiltonian systems: The Generalized Alignment Index (GALI) method, Physica D Nonlinear Phenomena 231, 30–54. Sussman, G. J. and Wisdom, J.: 2000, Structure and Interpretation of Classical Mechanics, The MIT Press, Cambridge, Massachussetts. Szebehely, V. G.: 1967, Theory of Orbits, The Restricted Problem of Three Bodies., Academic Press. Villac, B. and Aiello, J.: 2005, Mapping Long-Term Stability Regions Using The Fast Lyapunov Indicator, 15th AAS/AIAA Spaceflight Mechanics Meeting, Pasadena, CA : Jet Propulsion Laboratory, National Aeronautics and Space Administration, 2005, Copper Mountain, CO. URL: http://hdl.handle.net/ 2014/37515. Villac, B. F.: 2008, Using FLI maps for preliminary spacecraft trajectory design in multi-body environments, Celestial Mechanics and Dynamical Astronomy 102, 29–48. Villac, B. and Lara, M.: 2005, Stability Maps, Global Dynamics And Transfers, 2005 AAS/AIAA Astrodynamics Specialist Conference, Pasadena, CA : Jet Propulsion Laboratory, National Aeronautics and Space Administration, 2005, Lake Tahoe, CA. URL: http://hdl.handle.net/2014/37534. Wiegert, P., Innanen, K. and Mikkola, S.: 2000, The stability of quasi satellites in the outer solar system, Astronomical Journal 119, 1978–1984. DOI: 10.1086/301291. Wiesel, W. E.: 1993, Stable orbits about the martian moons, Journal of Guidance, Control, and Dynamics 16(3), 434–440. Wolfram Research, Inc.: 2011, Mathematica, version 8.0 edn, Wolfram Research, Inc., Champaign, Illinois. Zwillinger, D.: 1997, Handbook of Differential Equations, 3rd edn, Academic Press.
92
