+
2-Body Gravitation

Developed by Walter Freeman and Kelly Roos - Published June 5, 2019

In this exercise set, the student builds a computational model of a gravitational 2-body system consisting of a mass *m* gravitationally bound to a larger, stationary mass *M*, and then graphs and analyzes the resulting orbital motion. The student is guided in exploring the accuracy and stability of the computational algorithms used. Two numerical approaches are employed, the simple Euler method and the Euler-Cromer method.
Subject Areas Mechanics and Astronomy First Year and Beyond the First Year MATLAB, Python, Mathematica, Glowscript, and Spreadsheet Students who complete this set of exercises will be able to - model a dynamical system as a set of coupled differential equations, and use numerical integration to solve them (**Exercises 1, 2**) - probe the stability and accuracy of numerical simulations, and recognize the advantages of sympletic integrators (**Exercise 2**) - connect the parameters of a computer simulation to empirical observations of real-world systems, like the Earth-Moon system (**Exercise 3**) - discuss basic orbital dynamics as manifest in their numerical simulations (**Exercises 1-5**) - use the conservation of energy to probe the stability of a numerical simulation (**Exercise 4**) - identify physical systems that are difficult to simulate because they involve large differences in timescales (**Exercise 5**). 90-120 min

These exercises are not tied to a specific programming language. Example implementations are provided under the Code tab, but the Exercises can be implemented in whatever platform you wish to use (e.g., Excel, Python, MATLAB, etc.).

### 0. Describing orbits In an elliptical orbit, the satellite's distance from the object that it orbits changes. The degree to which this happens is called the *eccentricity* of the orbit. A circular orbit has zero eccentricity, while a strongly elongated orbit like a comet has zero eccentricity. We will use a few words to describe points on the orbit. The points nearest and furthest from the body being orbited (the host) are called *apsides* (singular, *apsis*). The point of closest approach is called the *periapsis* or *pericenter*. If the object being orbited is a star, it may also be called the *periastron*; if it is the Sun, the *perihelion*; and if it is the Earth, the *perigee*. Likewise, the furthest point from the host is the *apoapsis* or *apocenter*; for specific hosts, this may also be called *apastron*, *aphelion*, or *apogee*. ### 1. **Creating the model** The physics of a small satellite orbiting a large body are the same, whether this is the Moon orbiting the Earth, or a planet orbiting the Sun. As a reminder, if the large body is at the origin and the satellite is at position $(x,y)$, the acceleration of the satellite is given by $$a_x = \frac{-GMx}{r^3}$$ $$a_y = \frac{-GMy}{r^3}$$ where $r = \sqrt{x^2 + y^2}$. **Write a program that allows you to simulate this motion numerically using the Euler-Cromer algorithm and visualize it.** Your program should be able to both generate a picture or animation of the orbit and a plot of the magnitude of the orbiting object's position vector vs. time. It is always a good idea to test any computer model by simulating a system whose behavior you are familiar with, to ensure that your simulation reproduces all the features that it should. To test your program, consider the most familiar orbiting system we have: Earth's motion around the Sun. This system has the following properties: * The Sun has a mass of $1.99 \times 10^{30}$ kg, and moves only slightly compared to its distance from the Earth. (Thus, you can consider it to be stationary.) * Earth orbits the Sun in a nearly circular orbit at a distance of 1 AU = $1.5 \times 10^{11}$ m First, **determine initial conditions for the Earth: what should its starting position and velocity vectors be?** It will be useful to draw a cartoon of Earth's orbit, place the Earth somewhere in it, and then draw an arrow representing the Earth's velocity vector at that point. To determine the magnitude of Earth's velocity, think about the following: * How much time does it take for the Earth to complete one orbit? * How far does Earth move during that time? **Simulate the Earth's orbit around the Sun and confirm that you see the behavior that you expect**: a closed orbit that is nearly circular and that takes the expected amount of time to travel around the Sun. ### 2. **The Euler and Euler-Cromer algorithms** Now, investigate the behavior of your simulation using by comparing your Euler-Cromer results from the first Exercise to the the those produced using the simple Euler method. **Experiment with different values of the time step $\Delta t$ for the two algorithms.** How do their behaviors differ for large $\Delta t$? What about at small $\Delta t$? (Which algorithm works better?) ### 3. **The Moon's orbit around the Earth** Another familiar system is the Moon's orbit around the Earth. However, its orbit is more eccentric. The Moon's orbit has the following parameters: * Mass of the Earth: $5.97 \times 10^{24}$ kg * Lunar apogee (point furthest from Earth): $4.05 \times 10^8$ km * Lunar perigee (point closest to Earth): $3.63 \times 10^8$ km As before, use initial conditions in which the Moon's displacement and velocity vectors are perpendicular. This corresponds to either apogee or perigee (why?) However, here you don't know the initial velocity that will reproduce the actual orbit of the Moon. So, you'll need to find it. **Place the Moon at either apogee or perigee, and then by trial and error, adjust the initial velocity so that you obtain the correct value for the other apsis.** Then, examine the orbital period of the Moon. Does this match with what you know about the Moon's motion? In particular, **investigate the length of the orbital period and compare it to what we actually observe from the Moon.** ### 4. **Monitoring the orbital energy** Modify your code to keep track of the specific kinetic energy, specific potential energy, and total specific orbital energy of the Moon. (The *specific energy* is the energy per unit mass of an orbiting body.) In particular: * $T_s = \frac{1}{2}(v_x^2 + v_y^2)$ * $U_s = -\frac{GM}{r}$ * $E_s = T_s + U_s$ **For an eccentric orbit like the Moon's, plot the kinetic energy, potential energy, and total energy vs. time. Can you identify the perigee and apogee points from your plots?** Then play with your initial conditions by slowly increasing the Moon's initial velocity. **What sorts of orbits do you see if the total energy is very negative? Only slightly negative? Positive?** ### 5. **Timesteps and timescales** Determine roughly how many *timesteps per orbit* you need to simulate the motion of the Moon around the Earth with acceptable results. (What this means is up to you!) Now, simulate the orbit of Halley's comet around the Sun. (You'll need to change the central mass back to $M_\odot = 1.99\times 10^{30}$ kg. It has an extremely eccentric orbit, with an aphelion of $5.25\times 10^{12}$ m and a perihelion of $8.77\times 10^{10}$ m. As before, you'll need to use trial and error to find the initial velocity at aphelion that gives you the correct perihelion distance. **How many timesteps per orbit are required to get an acceptable simulation of Halley's comet around the Sun? Compare this to the number needed for the Moon; why are they different?** One of the common difficulties in computer simulation is the presence of dissimilar timescales: when the simulation needs to simulate both processes that occur very quickly and those that occur very slowly. **How is that present in the Halley's comet simulation, and why does it require you to do so much more work to simulate one orbit well?**