Astro2
October 18, 2016
In [1]: import pandas
import altair
from astropy import units as u
from numpy import sin, cos, arccos, sqrt
First, we load the data and format the types accordingly:
In [2]: data = pandas.read_table('moon.dat', sep='\s+')
data['date'] = pandas.to_datetime(data['date'], format='%d-%m-%Y')
Now we construct a new column containing the respective angle shifts:
In [3]: def angle_between_points(i,j):
alpha_i = data.iloc[i]['rect'] * u.hourangle
alpha_j = data.iloc[j]['rect'] * u.hourangle
delta_i = data.iloc[i]['decl'] * u.degree
delta_j = data.iloc[j]['decl'] * u.degree
angle = arccos(sin(delta_i)*sin(delta_j) +
cos(delta_i)*cos(delta_j)*cos(alpha_i-alpha_j))
return angle.to(u.rad).value
angle_shifts = [None]
rows = len(list(data.iterrows()))
for i in range(1, rows):
angle_shifts.append(angle_between_points(i, i-1))
data['shifts'] = pandas.Series(angle_shifts)
data.head()
Out[3]: date rect decl shifts
0 2011-07-01 6.96722 20.91167 NaN
1 2011-07-02 7.91667 17.78417 0.240692
2 2011-07-03 8.84167 13.60361 0.244187
3 2011-07-04 9.74028 8.63583 0.246491
4 2011-07-05 10.61861 3.18389 0.247641
In [4]: altair.Chart(data).mark_line().encode(x='date', y='shifts')
1
Since by Kepler’s second law we have:
r2
max · ϕmin = r2
min · ϕmax
Let us denote the ratio of rmax and rmin as k:
k =
rmax
rmin
=
ϕmax
ϕmin
Then using rmin = a(1 − e) and rmax = a(1 + e) we get:
k =
1 + e
1 − e
and finally:
e =
k − 1
k + 1
In [5]: phi_min = min(data.dropna()['shifts'])
phi_max = max(data.dropna()['shifts'])
k = sqrt(phi_max/phi_min)
e = (k-1)/(k+1)
print e
2
0.0529938798589
3