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