Orbits

One area of physics where computers make an essential contribution is that of calculating planetary orbits. Whereas if two bodies orbit each other following Newton's law of gravity their trajectories can be calculated analytically, add a third body and the problem can only be solved computationally.

That with two bodies the problem is solvable analytically is a useful check for any program claiming the same capability.

The program shown here is quite long, but still under ninety lines so it is included in its entirety. It defines a new "class", or data type, to hold the essential elements of a celestial body, and uses python's turtle graphics to plot the orbits.

The underlying physics is that the gravitational force between two bodies is attractive, and proportional to the product of their masses divided by the square of their separation. This is then resolved into components parallel to the x and y axes, which requires further multiplication by the separation in x divided by the total separation, and similarly for y.

The integration algorithm chosen is the Verlet algorithm with velocities.

The code lacks any concept of units, for the qualitative results will not depend on whether masses are measured in kg, tons, or Earth masses, nor on the precise value of G. It is also 2D only, which is a reasonable first approximation to our solar system, slightly quicker to code, and a lot easier to visualise on a 2D display. The art of Physics is always to make approximations which make a problem easier, but do not destroy the behaviour one is trying to study.

#!/usr/bin/python3
import turtle
from math import sqrt

class Body:
    def __init__(self):
        self.mass=1
        self.x=self.y=0
        self.vx=self.vy=0
        self.ax=self.ay=0
        self.marker=turtle.Turtle()
        self.marker.speed(0)

def summary(bodies):
    E=0
    mx=my=0
    for a in bodies:
        E+=0.5*a.mass*(a.vx**2+a.vy**2)
        mx+=a.mass*a.vx
        my+=a.mass*a.vy
        for b in bodies:
            if (a==b): continue
            r=sqrt((a.x-b.x)**2+(a.y-b.y)**2)
            E+=-0.5*a.mass*b.mass/r
    print("Total energy:   %f"%E)
    print("Total momentum: (%f,%f)"%(mx,my))
    
bodies=[]
b=Body()
b.mass=10
b.x=-5
bodies.append(b)
b=Body()
b.x=10
b.vy=0.5
bodies.append(b)

# Adjust first body's velocity for zero total momentum
mx=my=0
for i in range(1,len(bodies)):
    mx+=bodies[i].mass*bodies[i].vx
    my+=bodies[i].mass*bodies[i].vy

bodies[0].vx=-mx/bodies[0].mass
bodies[0].vy=-my/bodies[0].mass
    
summary(bodies)

dt=0.02
t=0

turtle.setworldcoordinates(-20,-20,20,20)
turtle.tracer(10,0)
for a in bodies:
    a.marker.penup()
    a.marker.goto(a.x,a.y)
    a.marker.pendown();
    
while(t<300):
    for a in bodies:
        fx=fy=0
        for b in bodies:
            if (a==b): continue
            r=sqrt((a.x-b.x)**2+(a.y-b.y)**2)
            fx+=-a.mass*b.mass*(a.x-b.x)/r**3
            fy+=-a.mass*b.mass*(a.y-b.y)/r**3
        a.vx=a.vx+0.5*(a.ax+fx/a.mass)*dt
        a.vy=a.vy+0.5*(a.ay+fy/a.mass)*dt
        a.ax=fx/a.mass
        a.ay=fy/a.mass
    for a in bodies:
        a.x=a.x+a.vx*dt+0.5*a.ax*dt*dt
        a.y=a.y+a.vy*dt+0.5*a.ay*dt*dt
        a.marker.goto(a.x,a.y)
    t=t+dt

print("Finished")
summary(bodies)
input('Press ENTER to close window')

This program prints the total energy and momentum at the start and end. Both should be conserved, but, in practice, the energy will drift slightly.

The orbits for the initial conditions as above should be elliptical, and should repeat exactly. That test the program does seem to pass.

Try increasing the time-step used by the integrator, dt, from 0.02 to 0.2. Note that the energy drift is now much greater. Increase it further to 0.5. The result might now be many times faster than the first version, but it is not as correct.

This code can do more than just two-body orbits. It can certainly do three bodies ones (and more), and there is a further page devoted to three body orbits.

But before moving on to that, try a couple more tests to show that this code is producing credible results. After resetting dt to 0.02,

• Increase vy for the second mass from 0.5 to 1, and the maximum simulation time from 300 to 800. Note that the total energy is still less than zero, so the system is stil bound, but the orbit is now very large.
• Increase vy for the second mass further to 1.2. Note that the total energy is now positive. There should not be a bound orbit, but the bodies should move infinitely far apart.

Now consider the three body orbits.