%pylab inline
import pylab as pl
import numpy as np
import math
tau = 2.0 * math.pi
halfpi = 0.5 * math.pi
degree = tau / 360.0
minute = degree / 60.
second = minute / 60.
import ephem
pole = ephem.FixedBody()
pole._ra = '0' # hours
pole._dec = '90' # degrees
def pole_positions(dates):
theta = np.zeros(len(dates))
r = np.zeros(len(dates))
for i, d in enumerate(dates):
pole.compute(d)
theta[i] = pole.ra
r[i] = (halfpi - pole.dec) / degree
return theta, r
start_date = ephem.Date((-162, 1, 1))
end_date = ephem.Date((2200, 1, 1))
print float(start_date), float(end_date)
dates = np.linspace(start_date, end_date, 365.0)
theta, r = pole_positions(dates)
plot = pl.figure().add_subplot(111, polar=True)
plot.plot(theta, r)
# When will it be closest?
# TODO: this makes no sense; re-do this as Polaris later?
print min(theta)
print np.where(theta==min(theta))
print dates[theta==min(theta)]
closest_dates = dates[theta==min(theta)]
print ephem.Date(closest_dates[0])
start_date = ephem.Date((1950, 1, 1))
end_date = ephem.Date((2050, 1, 1))
dates = range(int(start_date), int(end_date))
theta, r = pole_positions(dates)
plot = pl.figure().add_subplot(111, polar=True)
plot.plot(theta, r)
plot.axis([None, None, 0.0, 0.20])
plot = pl.figure().add_subplot(111, polar=True)
plot.plot(theta, r)
plot.axis([None, None, 0.0, 0.02])
plot = pl.figure().add_subplot(111)
plot.plot(theta, r)
plot.axis([89.0 * degree, 91.0 * degree,
0.0, 0.1])