5

I'm trying to write a code that plots the elliptical paths of an object using the equation for the ellipse r=a(1-e^2)/(1+e*cos(theta)). I'd also like this data to be put into an array for other use.

from numpy import *#Imports Python mathematical functions library
import matplotlib.pyplot as plt #Imports plot library
from pylab import *

a = 5
e = 0.3
theta = 0
while theta <= 2*pi:
    r = (a*(1-e**2))/(1+e*cos(theta))
    print("r = ",r,"theta = ",theta)
    plt.polar(theta, r)
    theta += pi/180

plt.show()

The code spits out correct values for r and theta, but the plot is blank. The polar plot window appears, but there is nothing plotted.

Please help. Thanks in advance.

maelstromscientist
  • 551
  • 9
  • 17
  • One quickly obvious flaw is theta being incremented by 180 degrees (in radians) - don't you want smaller steps, say 1 degree? – DarenW Oct 13 '12 at 02:36

2 Answers2

9

Do not call plt.polar once for every point. Instead, call it once, with all the data as input:

import numpy as np #Imports Python mathematical functions library
import matplotlib.pyplot as plt #Imports plot library
cos = np.cos
pi = np.pi

a = 5
e = 0.3
theta = np.linspace(0,2*pi, 360)
r = (a*(1-e**2))/(1+e*cos(theta))
plt.polar(theta, r)

print(np.c_[r,theta])

plt.show()

enter image description here

By the way, numpy can do the calculation as a two-liner, instead of using a while-loop:

theta = np.linspace(0,2*pi, 360)   # 360 equally spaced values between 0 and 2*pi
r = (a*(1-e**2))/(1+e*cos(theta))

This defines theta and r as numpy arrays (rather than single values).

unutbu
  • 842,883
  • 184
  • 1,785
  • 1,677
2

I think you need to do points.append([theta,r]) then at the end plt.polar(points) ... that makes a kinda neat design too

from numpy import *#Imports Python mathematical functions library
import matplotlib.pyplot as plt #Imports plot library
from pylab import *

a = 5
e = 0.3
theta = 0

points = []
while theta <= 2*pi:
    r = (a*(1-e**2))/(1+e*cos(theta))
    print("r = ",r,"theta = ",theta)
    points.append((theta, r))
    theta += pi/180
#plt.polar(points) #this is cool but probably not what you want
plt.polar(*zip(*points))
plt.show()
Joran Beasley
  • 110,522
  • 12
  • 160
  • 179
  • I was hoping this was Jupiter, but the eccentricity is *way* too big for that ... – mgilson Sep 11 '12 at 18:34
  • Our answers are basically the same, but our plots are completely different! :) I'm a little perplexed about how to fix (yours or mine?), but since the theta's go from 0 to 2pi, shouldn't the plot make a single orbit? – unutbu Sep 11 '12 at 18:50
  • Ah, our graphs would be the same if you use `plt.polar(*zip(*points))` so `plt.polar` receives two arguments instead of a single list. I'm not sure what `plt.polar` is doing when given a list of tuples... – unutbu Sep 11 '12 at 18:54
  • lol oh I thought that was just the plot ... I thought it looked neato – Joran Beasley Sep 11 '12 at 18:59