Tibeuleu/KozaiLidov
1
0
Fork 0
Code Activity

Initial conditions for 3 bodies

Browse Source
This commit is contained in:
Thibault Barnouin
2021-11-05 14:31:30 +01:00
parent 0d4cb52421
commit 6987baaf3a
7 changed files with 33 additions and 22 deletions
Download Patch File Download Diff File Expand all files Collapse all files

22
main.py
View File

@@ -6,18 +6,24 @@ import matplotlib.pyplot as plt
from lib.integrator import frogleap
from lib.objects import Body, System
globals()['G'] = 6.67e-11 #Gravitational constant in SI units
globals()['Ms'] = 2e30 #Solar mass in kg
globals()['au'] = 1.5e11 #Astronomical unit in m
def main():
#initialisation
m = np.array([1, 1, 1e-5])
m = np.array([1., 1., 0.1])*Ms # Masses in Solar mass
a = np.array([1., 1., 5.])*au # Semi-major axis in astronomical units
psi = np.array([0., 0., 80.])*np.pi/180. # Inclination of the orbital plane in degrees
x1 = np.array([-1, 0, 0])
x2 = np.array([1, 0, 0])
x3 = np.array([100, 0, 20])
x1 = np.array([-1., 0., 0.])*a[0]
x2 = np.array([1., 0., 0.])*a[1]
x3 = np.array([np.cos(psi[2]), 0., np.sin(psi[2])])*a[2]
q = np.array([x1, x2, x3])
v1 = np.array([0, -0.35, 0])
v2 = np.array([0, 0.35, 0])
v3 = np.array([0, 20., 0])
v1 = np.array([0, -np.sqrt(G*Ms/np.sqrt(np.sum(x1**2))), 0])
v2 = np.array([0, np.sqrt(G*Ms/np.sqrt(np.sum(x2**2))), 0])
v3 = np.array([0, np.sqrt(G*Ms*(2./np.sqrt(np.sum(x3**2))-1./a[2])), 0])
v = np.array([v1, v2, v3])
bodylist = []
@@ -27,7 +33,7 @@ def main():
dyn_syst.COMShift()
duration, step = 100, 0.01
E, L = frogleap(duration, step, dyn_syst, recover_param=True)#, display=True)
E, L = frogleap(duration, step, dyn_syst, recover_param=True, display=True)
fig1 = plt.figure(figsize=(30,15))
ax1 = fig1.add_subplot(111)
ax1.plot(np.arange(E.shape[0])/duration, E, label=r"$E_m$")