Rodney Dunning's Home Page | Research and Scholarship | VPython | Kepler's 1st & 2nd Laws
Introduction
Working from observational data obtained by Tycho Brahe, Kepler discovered three empirical rules obeyed by each planet:
- The law of orbits: Each planetary orbit is an ellipse with the Sun at one focus.
- The law of areas: A radius vector connecting the Sun and a given planet sweeps out equal areas in equal time periods.
- The law of periods: The ratio of the orbital periods of any two planets varies as the three-halves power of the ratio of their semi-major axes.
Kepler's laws provided a strong source of motivation to Isaac Newton in the later's search for a set of physical laws to explain and predict the motion of the planets.
Description
The animation shows the Sun, a yellow sphere, in the center of the window, and a planet, a red sphere, on the left side of the window. By default, the planet is always drawn 150 pixels from the Sun, and begins its orbit at the apocenter position. The animation ends when the planet returns to its initial position.
A table with basic orbital data appears in the upper-right-hand corner of the scene. The table records the elapsed time in days, the current distance of the planet from the Sun in AU, the current speed of the planet in AU/yr, the current acceleration of the planet in AU/yr2, and the number of orbits completed.
For each orbit, the program draws twelve radius vectors from the Sun to the planet, dividing the orbital period into twelve equal parts. According to Kepler's second law, each "wedge" has the same surface area. (Note: the radius lines are drawn only for the first complete orbit.)
The program uses a simple Verlet (leap-frog) algorithm to integrate the planet's orbit. The number of steps is equal to 1000n, where n is the number of orbits selected by the user. The orbit is viewed from a heliocentric reference frame.
When the program is executed, the user will be prompted for the following information:
- The semi-major axis of the planet, in AU.
- The eccentricity of the planet's orbit.
- The number of complete orbits to draw.
- Whether to draw the radius vectors that illustrate Kepler's 2nd law.
Screen Shot
The program draws the orbit, with radius vectors drawn for every twelfth part of the
orbital period.
Suggested Use
In my introductory astronomy course, I show the animation before describing Kepler's laws. On the initial viewing, I turn off the radius vectors. I ask students to qualitatively describe the shape of the orbit and how the speed of the planet varies throughout the orbit. On subsequent viewings, I turn on the radius vectors, and lead the students in a discussion and analysis of Kepler's laws.
Hints:
- The Verlet integration scheme will be accurate for any reasonable simulation. Integration errors will creep in if the planet gets too close to the Sun at pericenter, so beware high-eccentricity orbits. Errors will be expressed in a planet failing to maintain a closed orbit.
- Click anywhere in the animation scene to pause the animation. Click again to restart the animation.
Source Code
The source code appears between the bars below. Copy and paste the code into the Python IDLE environment, and hit F5 to run the code.
from visual import *
from __future__ import division
"""
An illustration of Kepler's first and second laws.
Copyright, 2006.
Rodney Dunning
Assistant Professor of Physics
Longwood University
"""
#-------------------------------------------
# Prompt user for semi-major axis and
# eccentricity.
#-------------------------------------------
print "This animation illustrates Kepler's first"
print "and second laws.\n\n"
a = -1
while(a <=0):
a = input("Enter the semi-major axis in Astronomical Units: ")
if(a < 0.1 or a > 1000):
print "\n*** INPUT ERROR ***\n"
print "The semi-major axis must be between 0.1 and 1000!\n"
e = -1
while(e < 0 or e >= 0.991):
e = input("\nEnter the eccentricity: ")
if(e < 0 or e >= 0.991):
print "\n*** INPUT ERROR ***\n"
print "The eccentricity must be between 0 and 0.99!\n"
n = -1
while(n <= 0):
n = input("\nEnter the number of full orbits to draw: ")
if n <= 0:
print "\n*** INPUT ERROR *** \n"
print "You must enter at least one (1) full orbit!\n"
lines = -1
while lines != 0 and lines != 1:
lines = input("\nDraw radial lines for Kepler's 2nd law?\n"
"1 = Yes, 0 = No: ")
if lines !=0 and lines != 1:
print "\n*** INPUT ERROR *** \n"
print "Please 1 for Yes or 0 for No!"
print "\n THANK YOU. Please find the animation window."
#--------------------------------------------------
# Scene attributes and functions
#--------------------------------------------------
scene.title = "Kepler's 1st and 2nd Laws"
scene.x = 0
scene.y = 0
scene.width = 600
scene.height = 600
scene.range = (1.5*a*(1+e),1.5*a*(1+e),1.5*a*(1+e))
scene.autoscale = 0 ##0 means autoscaling is OFF
scene.userzoom = 1 ##0 means user cannot zoom
scene.userspin = 1 ##0 means user cannot spin
scene.lights = [vector(0,0,1)]
scene.ambient = 0.5
scene.label = label(visible=1,
pos=(0,0,0),
xoffset = 0,
yoffset = 0,
text = "Kepler's 1st and 2nd Laws")
##scene.mouse_label = label(visible=1,
## pos=(0,100,0),
## text=("Mouse position"))
def return_mouse_pos(scene):
scene.mouse_label.text = ("mouse\n x: %3.2f y: %3.2f"
%(scene.mouse.pos.x, scene.mouse.pos.y))
def check_for_pause(scene): #checks for pause request
if scene.mouse.clicked:
scene.mouse.getclick()
pause(scene)
def pause(scene):
while 1:
## return_mouse_pos(scene)
if scene.mouse.clicked:
scene.mouse.getclick()
break
#---------------------------------------------------------
# Variables for keeping up with the elasped time
#---------------------------------------------------------
t = 0 # elapsed time
t_max = n * sqrt(a**3) # one full orbit
dt = t_max / (n*1000)
dt_2 = dt / 2 # for the Verlet loop below
t_line = 0 # keeps the time since the last line was drawn
num_lines = 0
#---------------------------------------------------------
# Dynamical variables
#---------------------------------------------------------
GM = 39.4784176 # gravitational constant times mass of the sun
g = GM / a
vx = 0 # initial x-component of the velocity
vy = sqrt(g * (1 - e) / (1 + e) ) # initial y-component of the velocity
vz = 0 # initial z-component of the velocity
x = a*(1+e) # aphelion distance; initial x position
y = 0 # initial y position
z = 0 # initial z position
#---------------------------------------------------------
# Create a sun and planet
#---------------------------------------------------------
"""
The sun is always centered in the animation window.
The planet always begins at x = 150. This distance
is scaled differently depending on the number of
AUs given by the user for the semi-major axis (SMA).
"""
sun = sphere(pos=vector(0,0,0),
radius = 0.1*a,
color = color.yellow)
planet = sphere(pos=vector(x,y,z),
vel=vector(vx,vy,vz),
radius = 0.05*a,
color = color.red,
orbit = curve(color=color.white))
#---------------------------------------------------------
# Create labels for the period, velocity, etc.
#---------------------------------------------------------
info_label = label(pos=vector(-a,0.75*a,0),
border=1,
box=0,
opacity = 0,
text="INFORMATION")
param_label = label(pos=vector(0,-1.25*a,0),
border = 1,
box = 0,
opacity = 0,
text="PARAMETERS")
#---------------------------------------------------------
# Animation loop
#---------------------------------------------------------
ready = 0
while not ready:
ready = scene.mouse.getclick()
scene.label.visible = 0
scene.label.text = "Finished.\n"
while (t <= t_max):
rate(100)
check_for_pause(scene)
t += dt
t_line += dt
# Using the Verlet algorithm
planet.pos += planet.vel * dt_2
r = sqrt(planet.pos.x**2 + planet.pos.y**2)
r3 = r**3
planet.acc = -GM * planet.pos / r3
planet.vel += planet.acc * dt
planet.pos += planet.vel * dt_2
planet.orbit.append(pos=planet.pos)
# Update labels
v = sqrt(planet.vel.x**2 + planet.vel.y**2)
acc = -GM / r**2
info_label.text = ("time: %3.3f days \n"
"distance: %3.3f AU \n"
"speed: %3.3f AU/yr\n"
"acc: %3.3f AU/yr^2\n"
"orbits: %3.3f \n"
%(t*365.24, r, v, acc,
n*t/t_max))
param_label.text = "a: %3.3f AU, e: %3.3f, %3.2 orbits" %(a,e,n)
# The second condition gaurentees only 12 lines
# drawn even if the user requests more than one
# full orbit--otherwise the animation becomes
# confusing.
if lines and num_lines < 12:
## The following code draws 12 = 1000/83 radial
## from the Sun to the planet
if(abs(t_line - 83*dt) < dt): #draw a line
line = curve(pos = [sun.pos,planet.pos],
color = color.white)
t_line = 0
num_lines += 1
scene.label.visible = 1