# simulate Kepler's first law (*only*)

# using polar coordinates and formulas given here to calculate r and p:
# http://en.wikipedia.org/w/index.php?title=Kepler%27s_laws_of_planetary_motion&oldid=473666923#First_Law

import turtle
import math

turtle.tracer(5)

# size of ellipse
a = 100
# eccentricity - how squished the ellipse is
epsilon = 0.66
# simulation steps - stepping through theta values
NSTEPS = 1000000

# semi-latus rectum
# make your own joke
p = a * (1 - epsilon ** 2)

# sun
turtle.dot(42, 'yellow')

turtle.pu()
turtle.ht()

for theta in range(NSTEPS):
	theta_rad = math.radians(theta)
	denom = 1 + epsilon * math.cos(theta_rad)
	r = p / denom

	# go to polar coords (r, theta)
	# could use turtle.lt too, since turtle.home below resets heading
	turtle.seth(theta)
	turtle.fd(r)

	# only want the turtle to appear when it's on the orbital path
	turtle.st()
	turtle.dot(1, 'red')
	turtle.ht()
	turtle.home()