In [1]:
%matplotlib inline

In [2]:
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.dates as mdates

from astropy.time import Time

import subprocess

# Larger figure size
fig_size = [10, 6]
plt.rcParams['figure.figsize'] = fig_size

In [3]:
mjd_unixtimestamp_offset = 10587.5
seconds_in_day = 3600 * 24

def mjd2unixtimestamp(m):
    return (m - mjd_unixtimestamp_offset) * seconds_in_day

def unixtimestamp2mjd(u):
    return u / seconds_in_day + mjd_unixtimestamp_offset

In [4]:
def load_orbit_file(path):
    ncols = 8
    data = np.fromfile(path, sep=' ')
    return data.reshape((data.size // ncols, ncols))

Keys for each of the columns in the orbit (Keplerian state) report.


In [5]:
utc = 0
sma = 1
ecc = 2
inc = 3
raan = 4
aop = 5
ma = 6
ta = 7

In [6]:
fig_peri = plt.figure(figsize = [15,8], facecolor='w')
fig_apo = plt.figure(figsize = [15,8], facecolor='w')
fig3 = plt.figure(figsize = [15,8], facecolor='w')
fig4 = plt.figure(figsize = [15,8], facecolor='w')
fig4_rap = plt.figure(figsize = [15,8], facecolor='w')
fig5 = plt.figure(figsize = [15,8], facecolor='w')
fig6 = plt.figure(figsize = [15,8], facecolor='w')

sub_peri = fig_peri.add_subplot(111)
sub_apo = fig_apo.add_subplot(111)
sub3 = fig3.add_subplot(111)
sub4 = fig4.add_subplot(111)
sub4_rap = fig4_rap.add_subplot(111)
sub5 = fig5.add_subplot(111)
sub6 = fig6.add_subplot(111)

subs = [sub_peri, sub_apo, sub3, sub4, sub4_rap, sub5, sub6]

for file in ['dslwp_report.txt', 'dslwp_report_simple.txt']:
    orbit = load_orbit_file(file)

    t = Time(mjd2unixtimestamp(orbit[:,utc]), format='unix')

    sub_peri.plot(t.datetime, orbit[:,sma]*(1-orbit[:,ecc]))

    sub_apo.plot(t.datetime, orbit[:,sma]*(1+orbit[:,ecc]))
    sub3.plot(t.datetime, orbit[:,ecc])
    sub4.plot(t.datetime, orbit[:,aop])
    sub4_rap.plot(t.datetime, np.fmod(orbit[:,aop] + orbit[:,raan],360))
    sub5.plot(t.datetime, orbit[:,inc])
    sub6.plot(t.datetime, orbit[:,raan])

month_locator = mdates.MonthLocator()
day_locator = mdates.DayLocator()

for sub in subs:
    sub.set_xlabel('Time')
    sub.xaxis.set_major_locator(month_locator)
    sub.xaxis.set_major_formatter(mdates.DateFormatter('%Y-%m'))
    sub.xaxis.set_tick_params(rotation=45)

sub_peri.set_ylabel('Periapsis radius (km)')
sub_apo.set_ylabel('Apoapsis radius (km)')
sub3.set_ylabel('ECC')
sub4.set_ylabel('AOP (deg)')
sub4_rap.set_ylabel('RAOP (deg)')
sub5.set_ylabel('INC (deg)')
sub6.set_ylabel('RAAN (deg)')

sub_peri.set_title('Periapsis radius')
sub_apo.set_title('Apoapsis radius')
sub3.set_title('Eccentricity')
sub4.set_title('Argument of periapsis')
sub4_rap.set_title('Right ascension of periapsis')
sub5.set_title('Inclination')
sub6.set_title('Right ascension of ascending node')

for sub in subs:
    sub.legend(['Full force model', 'Only Moon and Earth point masses'])



In [7]:
mu_1 = 3.986004418e14 # Earth's gravitational constant
mu_0 = 4.9048695e12 # Moon's gravitational parameter

In [8]:
orbit = load_orbit_file('dslwp_report_earth.txt')
t = Time(mjd2unixtimestamp(orbit[:,utc]), format='unix')

In [9]:
plt.plot(t.datetime, orbit[:,raan])


Out[9]:
[<matplotlib.lines.Line2D at 0x7f707e52b6d8>]

In [45]:
plt.figure(figsize = [15,8], facecolor='w')
plt.plot(t.datetime, orbit[:,aop])
plt.xlabel('Time')
plt.ylabel('$\omega$ (degrees)')
plt.title('Argument of periapsis');



In [11]:
plt.figure(figsize = [15,8], facecolor='w')
plt.plot(t.datetime, orbit[:,inc])
plt.xlabel('Time')
plt.ylabel('$i$ (degrees)')
plt.title('Inclination');



In [12]:
Omega = orbit[:,raan]/180*np.pi
w = orbit[:,aop]/180*np.pi
i = orbit[:,inc]/180*np.pi
xi_1 = np.cos(Omega)*np.cos(w) - np.sin(Omega)*np.cos(i)*np.sin(w)
xi_2 = -np.cos(Omega)*np.sin(w) - np.sin(Omega)*np.cos(i)*np.cos(w)
xi_3 = np.sin(Omega)*np.sin(i)
beta_1 = xi_1**2
beta_2 = xi_2**2
beta_3 = xi_1*xi_2
beta_4 = xi_2*xi_3
beta_5 = xi_1*xi_3

In [13]:
a = 8960 # SMA in km
P = 2*np.pi*np.sqrt((1e3*a)**3/mu_0)/(24*3600) # period in days
e = 0.67 # mean eccentricity
i_avg = 0.64 # mean inclination
r_1 = 384400 # mean Earth-Moon distance

In [14]:
def integrate_delta(x):
    return np.concatenate((np.array([0]), np.cumsum(x[:-1]/P*np.diff(orbit[:,utc]))))

In [15]:
plt.figure(figsize = [15,8], facecolor='w')
plt.plot(t.datetime, orbit[:,ecc])
delta_e_secular = 15/4*np.pi*mu_1/mu_0*(a/r_1)**3*e*np.sqrt(1-e**2)*np.sin(i_avg)**2*np.sin(2*w)
delta_e_complete= -15*np.pi*mu_1/mu_0*(a/r_1)**3*e*np.sqrt(1-e**2)*beta_3
plt.plot(t.datetime, orbit[0,ecc]+integrate_delta(delta_e_complete))
plt.plot(t.datetime, orbit[0,ecc]+integrate_delta(delta_e_secular))
plt.legend(['GMAT', 'Long period + secular', 'Secular'])
plt.xlabel('Time')
plt.ylabel('$e$')
plt.title('Eccentricity perturbations');



In [46]:
plt.figure(figsize = [15,8], facecolor='w')
plt.plot(t.datetime, w)
delta_w_secular = 3/2*np.pi*mu_1/mu_0*(a/r_1)**3*1/np.sqrt(1-orbit[:,ecc]**2)*(5*np.cos(i)**2*np.sin(w)**2 + (1-orbit[:,ecc]**2)*(2-5*np.sin(w)**2))
delta_w_secular_simple = 3/2*np.pi*mu_1/mu_0*(a/r_1)**3*1/np.sqrt(1-e**2)*(5*np.cos(i_avg)**2*np.sin(w)**2 + (1-e**2)*(2-5*np.sin(w)**2))
delta_w_complete = 3*np.pi*mu_1/mu_0*(a/r_1)**3*(np.sqrt(1-orbit[:,ecc]**2)*(4*beta_1-beta_2-1) - 1/(np.sqrt(1-orbit[:,ecc]**2)*np.sin(i))*((1+4*orbit[:,ecc]**2)*beta_5*np.sin(w)+(1-orbit[:,ecc]**2)*beta_4*np.sin(w)))
w_secular_est = w[0]+integrate_delta(delta_w_secular)
w_secular_simple_est = w[0]+integrate_delta(delta_w_secular_simple)
w_complete_est = w[0]+integrate_delta(delta_w_complete)
plt.plot(t.datetime, np.fmod(w_complete_est, 2*np.pi))
plt.plot(t.datetime, np.fmod(w_secular_est, 2*np.pi))
plt.plot(t.datetime, np.fmod(w_secular_simple_est, 2*np.pi))
plt.legend(['GMAT', 'Long period + secular', 'Secular', 'Secular (simplified)'])
plt.xlabel('Time')
plt.ylabel('$\omega$ (rad)')
plt.title('Argument of periapsis perturbations');