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.