import ephem 
import numpy as np
import matplotlib.pyplot as plt
from matplotlib import cm
import datetime as dt
au=ephem.meters_per_au


sun=ephem.Sun()
mercury=ephem.Mercury()
venus = ephem.Venus()
#earth = ephem.Earth()
mars = ephem.Mars()
jupiter = ephem.Jupiter()
saturn = ephem.Saturn()
uranus = ephem.Uranus()
neptune = ephem.Neptune()
pluto = ephem.Pluto()

bodies=[sun,mercury,venus,mars,jupiter,saturn,uranus,neptune,pluto]

colors = [ cm.jet(x) for x in np.linspace(0,1, len(bodies) ) ]

plt.figure()
ax = plt.subplot(111, projection='polar')

dates=[dt.datetime(1950,1,1)+dt.timedelta(days=30*i) for i in range(12*100)]

tlons=[]
trs=[]
print dates[0],dates[-1]
for q in bodies:
  tlons.append([])
  trs.append([])


for date in dates:
  mydate=ephem.Date(date) #"%i/1/1"%year

  for idx,body in enumerate(bodies):
    body.compute(mydate)
    tlon,tr=float(body.hlon),body.sun_distance
    tlons[idx].append(tlon)
    trs[idx].append(tr)


for idx,body in enumerate(bodies):
  plt.plot(tlons[idx],trs[idx],color=colors[idx],label=body.name,lw=2)

ax.legend(bbox_to_anchor=(1.27, 1.12),fontsize=10)
#ax.set_rscale('log')
plt.savefig("planets.pdf")