In [1]:
from __future__ import print_function
%matplotlib inline
import numpy as np
import matplotlib as mpl
from mpl_toolkits.mplot3d import Axes3D
import matplotlib.pyplot as plt
from matplotlib import rc
rc('text', usetex=True)
plt.rc('font', family='serif')
from IPython.html.widgets import interact
In [2]:
%%bash
make
./LidovKozai
In [3]:
body_1_filename = 'body_001.dat'
body_2_filename = 'body_002.dat'
body_3_filename = 'body_003.dat'
t_1, x_1, y_1, z_1, v_x_1, v_y_1, v_z_1 = np.loadtxt(body_1_filename, unpack=True, skiprows=1)
t_2, x_2, y_2, z_2, v_x_2, v_y_2, v_z_2 = np.loadtxt(body_2_filename, unpack=True, skiprows=1)
t_3, x_3, y_3, z_3, v_x_3, v_y_3, v_z_3 = np.loadtxt(body_3_filename, unpack=True, skiprows=1)
In [4]:
num_steps = 10
plot_step_size =int(t_1[-1] - t_1[0]) / num_steps
def plot_solution( xy_angle=0, z_angle=45, third_body=True, COM_frame=True, plot_start=0):
fig = plt.figure()
ax = fig.gca(projection='3d', aspect="equal")
print("plot range: t= ",plot_start, " - ", plot_start + plot_step_size)
indices = np.where( ((t_1 > plot_start) & (t_1 < plot_start + plot_step_size)))
if COM_frame is True:
x_rel = 0
y_rel = 0
z_rel = 0
else:
x_rel = x_1[indices]
y_rel = y_1[indices]
z_rel = z_1[indices]
plt.scatter(0,0, marker="+")
ax.plot(x_1[indices] - x_rel, y_1[indices] - y_rel, z_1[indices] - z_rel, linestyle='-')
ax.plot(x_2[indices] - x_rel, y_2[indices] - y_rel, z_2[indices] - z_rel, linestyle='dotted')
if third_body is True:
ax.plot(x_3[indices] - x_rel, y_3[indices] - y_rel, z_3[indices] - z_rel, linestyle='--')
plt.xlim(-30,30)
plt.ylim(-30,30)
plt.legend((r"Body 1", r"Body 2", r"Body 3"))
else:
plt.xlim(-2, 2)
plt.ylim(-2, 2)
plt.legend((r"Body 1", r"Body 2"))
plt.xlabel(r"x")
plt.ylabel(r"y")
ax.view_init(z_angle, -xy_angle)
return
interact(plot_solution,
plot_start=(0,plot_step_size*np.floor(t_1[-1] / plot_step_size),plot_step_size),
xy_angle=(0.,90.),
z_angle=(0,90),
third_body=True,
COM_frame=True);
Below, we use the instantaneous $\mathbf{x}$ and $\mathbf{v}$ vectors to get the orbital elements at each time step.
This follows the equations of Montenbruck & Gill's "Satellite Orbits: Models, Methods and Applications" (section 2.2.4), which is clearer than most sources.
In [5]:
G = 1
mass_1 = 1
mass_2 = .01
mu = (mass_2 * mass_1 ) / (mass_2 + mass_1)
h_x = np.multiply((y_2 - y_1), (v_z_2 - v_z_1)) - np.multiply((z_2 - z_1), (v_y_2 - v_y_1))
h_y = np.multiply((z_2 - z_1), (v_x_2 - v_x_1)) - np.multiply((x_2 - x_1), (v_z_2 - v_z_1))
h_z = np.multiply((x_2 - x_1), (v_y_2 - v_y_1)) - np.multiply((y_2 - y_1), (v_x_2 - v_x_1))
h = np.sqrt( np.power(h_x, 2) + np.power(h_y,2) + np.power(h_z, 2))
W_x = h_x / h
W_y = h_y / h
W_z = h_z / h
r = np.sqrt( np.power( x_2 - x_1 ,2) + np.power( y_2 - y_1 ,2) + np.power( z_2 - z_1 ,2))
v = np.sqrt( np.power(v_x_2 - v_x_1 ,2) + np.power(v_y_2 - v_y_1 ,2) + np.power(v_z_2 - v_z_1 ,2))
rdotv = np.multiply(x_2 - x_1, v_x_2 - v_x_1) + np.multiply(y_2 - y_1, v_y_2 - v_y_1) + np.multiply(z_2 - z_1, v_z_2 - v_z_1)
i = np.arctan2(np.sqrt(np.power(W_x, 2) + np.power(W_y, 2)) , W_z)
Omega = np.unwrap(np.arctan2(W_x, -W_y))
p = np.divide(np.power(h, 2), G * mass_1)
a = np.power(np.divide(2,r) - np.divide(np.power(v,2), G * mass_1), -1)
n = np.sqrt( np.divide(G * mass_1, np.power(a,3) ) )
e = np.sqrt( 1 - np.divide(p, a))
E = np.arctan2( np.divide(rdotv , np.multiply(np.power(a,2), n)), 1 - np.divide(r, a) )
# M needs to be solved numerically
u = np.arctan2(z_2 - z_1, np.multiply((y_2-y_1), W_x) - np.multiply((x_2-x_1), W_y))
nu = np.arctan2( np.multiply(np.sqrt(1 - np.power(e,2)), np.sin(E)) , np.cos(E) - e)
omega = np.unwrap(u - nu)
Theta = np.multiply((1 - np.power(e,2)), np.power(np.cos(i),2))
In [6]:
plt.figure()
plt.plot(t_1, np.rad2deg(i))
plt.xlabel(r"$t$")
plt.ylabel(r"$i$")
plt.ylim(0, 90)
plt.figure()
plt.plot(t_1, np.rad2deg(Omega))
plt.ylabel(r"$\Omega$")
plt.figure()
plt.plot(t_1, v)
plt.ylim(ymin=0)
plt.xlabel(r"$t$")
plt.ylabel(r"$v$")
plt.figure()
plt.plot(t_1, a)
plt.ylim(ymin=0)
plt.xlabel(r"$t$")
plt.ylabel(r"$a$")
plt.figure()
plt.plot(t_1, e)
plt.ylim(0, 1)
plt.xlabel(r"$t$")
plt.ylabel(r"$e$")
plt.figure()
plt.plot(t_1, np.rad2deg(omega))
plt.ylim(-180, 180)
plt.xlabel(r"$t$")
plt.ylabel(r"$\omega$")
plt.figure()
plt.plot(np.rad2deg(omega), 1 - np.power(e,2))
plt.xlim(-180,180)
plt.ylim(0,1)
plt.xlabel(r"$\omega$")
plt.ylabel(r"$1 - e^2$")
plt.figure()
plt.plot(t_1, Theta)
plt.ylim(0,1)
plt.xlabel(r"$t$")
plt.ylabel(r"$\Theta$")
Out[6]: