# Simulating Planetary Orbits¶

According to Isaac Newton, the force of gravitational attraction between two objects is given by:

F = G\frac{m_1m_2}{d^2}

where G is the gravitational constant, m_1 and m_2 are the masses of the two objects, and d is the distance between them. In SI units, G has the value 6.67428 \times 10^{-11} N(m/kg)^2, so d is measured in meters, the masses are measured in kilograms, and the resulting F is in newtons.

Using this equation, Newton determined a formula for calculating how long it took an object to complete an orbit around a central mass. However, when dealing with three or more objects, it’s generally not possible to find a tidy formula to calculate what the three bodies will do.

Instead, such problems are tackled by numeric integration, a brute-force approach where you take all the object positions and velocities at time T, calculate the forces they exert on each other, update the velocities, and calculate the new positions at time T+\epsilon. Then you repeat this in a loop, stepping forward through time, and output or plot the results.

## Approach¶

To implement this in Python, we’ll use the turtle module to provide a graphical display, subclassing the Turtle class to create a Body class that will have additional attributes: mass for the object’s mass, vx and vy for its velocity, and px and py for its position.

An added method on Body, attraction(), will take another Body instance and return the X and Y components of the force exerted by the other body.

## Solution¶

  1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 #!/usr/bin/env python3 import math from turtle import * # The gravitational constant G G = 6.67428e-11 # Assumed scale: 100 pixels = 1AU. AU = (149.6e6 * 1000) # 149.6 million km, in meters. SCALE = 250 / AU class Body(Turtle): """Subclass of Turtle representing a gravitationally-acting body. Extra attributes: mass : mass in kg vx, vy: x, y velocities in m/s px, py: x, y positions in m """ name = 'Body' mass = None vx = vy = 0.0 px = py = 0.0 def attraction(self, other): """(Body): (fx, fy) Returns the force exerted upon this body by the other body. """ # Report an error if the other object is the same as this one. if self is other: raise ValueError("Attraction of object %r to itself requested" % self.name) # Compute the distance of the other body. sx, sy = self.px, self.py ox, oy = other.px, other.py dx = (ox-sx) dy = (oy-sy) d = math.sqrt(dx**2 + dy**2) # Report an error if the distance is zero; otherwise we'll # get a ZeroDivisionError exception further down. if d == 0: raise ValueError("Collision between objects %r and %r" % (self.name, other.name)) # Compute the force of attraction f = G * self.mass * other.mass / (d**2) # Compute the direction of the force. theta = math.atan2(dy, dx) fx = math.cos(theta) * f fy = math.sin(theta) * f return fx, fy def update_info(step, bodies): """(int, [Body]) Displays information about the status of the simulation. """ print('Step #{}'.format(step)) for body in bodies: s = '{:<8} Pos.={:>6.2f} {:>6.2f} Vel.={:>10.3f} {:>10.3f}'.format( body.name, body.px/AU, body.py/AU, body.vx, body.vy) print(s) print() def loop(bodies): """([Body]) Never returns; loops through the simulation, updating the positions of all the provided bodies. """ timestep = 24*3600 # One day for body in bodies: body.penup() body.hideturtle() step = 1 while True: update_info(step, bodies) step += 1 force = {} for body in bodies: # Add up all of the forces exerted on 'body'. total_fx = total_fy = 0.0 for other in bodies: # Don't calculate the body's attraction to itself if body is other: continue fx, fy = body.attraction(other) total_fx += fx total_fy += fy # Record the total force exerted. force[body] = (total_fx, total_fy) # Update velocities based upon on the force. for body in bodies: fx, fy = force[body] body.vx += fx / body.mass * timestep body.vy += fy / body.mass * timestep # Update positions body.px += body.vx * timestep body.py += body.vy * timestep body.goto(body.px*SCALE, body.py*SCALE) body.dot(3) def main(): sun = Body() sun.name = 'Sun' sun.mass = 1.98892 * 10**30 sun.pencolor('yellow') earth = Body() earth.name = 'Earth' earth.mass = 5.9742 * 10**24 earth.px = -1*AU earth.vy = 29.783 * 1000 # 29.783 km/sec earth.pencolor('blue') # Venus parameters taken from # http://nssdc.gsfc.nasa.gov/planetary/factsheet/venusfact.html venus = Body() venus.name = 'Venus' venus.mass = 4.8685 * 10**24 venus.px = 0.723 * AU venus.vy = -35.02 * 1000 venus.pencolor('red') loop([sun, earth, venus]) if __name__ == '__main__': main() 

## Code Discussion¶

The system described in the code consists of the Sun, Earth, and Venus, so the main() function creates three Body instances for each body and passed to the loop() function.

The loop() function is the heart of the simulation, taking a list of Body instances and then performing simulation steps forever. The time step chosen is one day, which works well for our Sun/Earth/Venus example. When you run the program, you can see how long it takes for the plot to complete an entire orbit; for Earth it’s the expected 365 days and for Venus it’s 224 days.

## Lessons Learned¶

Each simulation step requires calculating N * (N-1) distances and attractions, so the time complexity is O(N^2). On a laptop or desktop, the display will be visible changing up to around 20 objects. More efficient coding would let us handle more objects; we could rewrite the calculations in C or parallelize the code to divide the work in each step among multiple threads or CPUs. You could also adjust the timestep dynamically: if objects are far apart, a larger timestep would introduce less error, and the timestep could be shortened when objects are interacting more closely.

These techniques would increase our practical limit to hundreds (10^3) or thousands (10^4) of objects, but this means we can’t simulate even a small galaxy, which might contain tens of millions of stars (10^7). (Our galaxy is estimated to have around 200 billion stars, 2 \times 10^{11}.) Entirely different approaches need to be taken for that problem size; for example, the attraction of distant particles is approximated and only nearby particles are calculated exactly. The references include a survey by Drs. Trenti and Hut that describes the techniques used for larger simulations.

## References¶

http://www.scholarpedia.org/article/N-body_simulations
This survey, by Dr. Michele Trenti and Dr. Piet Hut, describes how the serious scientific N-body simulators work, using trees to approximate the attraction at great distances. Such programs are able to run in O(N log(N)) time.
http://ssd.jpl.nasa.gov/horizons.cgi
NASA’s Jet Propulsion Laboratory provides a system called HORIZONS that returns accurate positions and velocities for objects within the solar system. In the example code, the values used are only rough approximations; the orbital distances and planet velocities are set to the mean distances and their relative positions don’t correspond to any actual point in time – but they produce reasonable output.