There's at least one mistake here. You have this:
r = (1+1.414*math.sin(theta)*math.cos(theta)*-0.5*math.cos(theta)*math.cos(theta))**(1/6)*math.exp(-0.4714*theta)
and, according to the link, it should be this:
r = (1+1.414*math.sin(theta)*math.cos(theta) -0.5*math.cos(theta)*math.cos(theta))**(1/6)*math.exp(-0.4714*theta)
Also, be careful with expressions like (1/6)
in Python. If you're using Python2.X, then this will perform integer division and the result is 0
. To get floating point division as standard, you need to put
from __future__ import division
at the top of your script.