%pylab inline
import ephem
from mpl_toolkits.basemap import Basemap

names = ['Dubhe', 'Merak', 'Phecda', 'Megrez',
         'Alioth', 'Mizar', 'Alcor', 'Alcaid',]
stars = [ephem.star(name) for name in names]

def dipper_radec(*args):
    for star in stars:
        star.compute(*args)
    ra = [-degrees(star.a_ra) for star in stars]
    dec = [degrees(star.a_dec) for star in stars]
    return ra, dec
    
def plot_dipper(*args):
    ra, dec = dipper_radec(*args)
    width = 4000000
    m = Basemap(width=width, height=width, projection='aeqd',
                lat_0=dec[3], lon_0=ra[3])
    m.drawparallels(np.arange(-80,81,10))
    m.drawmeridians(np.arange(-180,180,10))
    x, y = m(ra, dec)
    m.plot(x, y, marker='o', linestyle='--', color='b')
    plt.title('date={} epoch={}'.format(*args))

plot_dipper('2000', '-6000')
show()
plot_dipper('2000', '2000')
show()
plot_dipper('2000', '10000')

plot_dipper('-50000', '2000')
show()
plot_dipper('2000', '2000')
show()
plot_dipper('50000', '2000')

plot_dipper('-50000', '-50000')
show()
plot_dipper('2000', '2000')
show()
plot_dipper('50000', '50000')

polaris = ephem.star('Polaris')

width = 9000000

m = Basemap(width=width, height=width, projection='aeqd',
            lat_0=70.0, lon_0=280.0)
m.drawparallels(np.arange(-80,81,10))
m.drawmeridians(np.arange(-180,180,10))

for year in range(-13001, 13000, 500):
    polaris.compute(str(year), epoch=str(year))
    x, y = m(degrees(polaris.a_ra), degrees(polaris.a_dec))
    m.plot(x, y, marker='o', linestyle='--', color='b')

width = 18000000

m = Basemap(width=width, height=width, projection='aeqd',
            lat_0=70.0, lon_0=280.0)
m.drawparallels(np.arange(-80,81,10))
m.drawmeridians(np.arange(-180,180,10))

for year in range(-50001, 50000, 500):
    polaris.compute(str(year), epoch=str(year))
    x, y = m(degrees(polaris.a_ra), degrees(polaris.a_dec))
    m.plot(x, y, marker='o', linestyle='--', color='b')

from ctypes import *
cdll.LoadLibrary('../ephem/_libastro.so')
lib = CDLL('../ephem/_libastro.so', use_errno=True)

lib.precess.argtypes = [
    c_double, c_double, POINTER(c_double), POINTER(c_double)
    ]

ra = []
dec = []

ra_double = c_double()
dec_double = c_double()

epoch1 = ephem.J2000
year = arange(-50001.0, 50000.0, 100)
for y in year:
    epoch2 = epoch1 + 365.25 * y
    ra_double.value = 0.0
    dec_double.value = 0.0
    lib.precess(epoch1, epoch2,
                byref(ra_double), byref(dec_double))
    ra.append(ra_double.value)
    dec.append(dec_double.value)

plot(year, ra)
show()
plot(year, dec)

lib.nutation.argtypes = [
    c_double, POINTER(c_double), POINTER(c_double)
    ]

deps = []
dpsi = []

deps_double = c_double()
dpsi_double = c_double()

epoch1 = ephem.J2000
year = arange(-50001.0, 50000.0, 100)
for y in year:
    epoch2 = epoch1 + 365.25 * y
    deps_double.value = 0.0
    dpsi_double.value = 0.0
    lib.nutation(epoch2, byref(deps_double), byref(dpsi_double))
    deps.append(deps_double.value)
    dpsi.append(dpsi_double.value)

plot(year, deps)
show()
plot(year, dpsi)

from skyfield.api import JulianDate
from skyfield import precessionlib

years = arange(-50000, 50000)
jd = JulianDate(tai=(years,))
plot(years, precessionlib.compute_precession(jd.tdb)[0][0])