import matplotlib.pyplot as plt
import numpy as np

coords = np.loadtxt('moon.dat', usecols=range(1, 3))
alpha = coords[:, 0]*15*np.pi/180   # rektascenzia v radianoch
delta = coords[:, 1]*np.pi/180   # deklinacia v radianoch

# Denny pohyb n je uhol o ktory sa teleso posunie na oblohe za den. Je 
# ekvivalentny uhlovej rychlosti omega, len je vyjadreny v inych jednotkach.
n = []
for i in range(1, len(alpha)):
    n.append(180/np.pi*np.arccos(np.sin(delta[i])*np.sin(delta[i - 1]) +
        np.cos(delta[i])*np.cos(delta[i - 1])*np.cos(alpha[i] - alpha[i-1])))
    
plt.plot(n)
plt.xlabel('time [days]')
plt.ylabel('n [°/day]')
plt.title("Moon's daily motion")

# Linearnu excentricitu odvodime z pomeru periheliovej a afeliovej uhlovej
# rychlosti a pouzitim druheho Keplerovho zakona, vid ukol1_8.jpg.
Q = np.sqrt(np.max(n)/np.min(n))
e = (Q - 1)/(Q + 1)
print('Vypocitana excentricita Mesiaca je {}.'.format(np.round(e, 3)))
print('Skutocna excentricita Mesiaca je 0.055.')