test run

lib/plots.py

@@ -134,9 +134,8 @@ def display_parameters(E,L,sma,ecc,parameters,savename=""):
 
     fig3 = plt.figure(figsize=(15,7))
     ax3 = fig3.add_subplot(111)
-    for i in range(len(L)):
-        ax3.plot(np.arange(sma[i].shape[0])*step[i]/yr, sma[i], label="a (step of {0:.2e}s)".format(step[i]))
-        ax3.plot(np.arange(ecc[i].shape[0])*step[i]/yr, ecc[i], label="e (step of {0:.2e}s)".format(step[i]))
+    ax3.plot(np.arange(sma[-1].shape[0])*step[-1]/yr, sma[i], label="a (step of {0:.2e}s)".format(step[-1]))
+    ax3.plot(np.arange(ecc[-1].shape[0])*step[-1]/yr, ecc[i], label="e (step of {0:.2e}s)".format(step[-1]))
     ax3.set(xlabel=r"$t \, [yr]$", ylabel=r"$a \, [au] \, or \, e$")
     ax3.legend()
     fig3.suptitle("Semi major axis and eccentricity "+title2)

main_concurrent.py

@@ -14,9 +14,9 @@ from lib.units import *
 def main():
     #initialisation
     m = np.array([1., 1., 1e-1],dtype=np.longdouble)*Ms/Ms  # Masses in Solar mass
-    a = np.array([1., 1., 5.],dtype=np.longdouble)*au/au   # Semi-major axis in astronomical units
-    e = np.array([0., 0., 0.],dtype=np.longdouble)   # Eccentricity
-    psi = np.array([0., 0., 0.],dtype=np.longdouble)*np.pi/180.    # Inclination of the orbital plane in degrees
+    a = np.array([1., 1., 5.],dtype=np.longdouble)/2.*au/au   # Semi-major axis in astronomical units
+    e = np.array([0., 0., 1./4.],dtype=np.longdouble)   # Eccentricity
+    psi = np.array([0., 0., 45.],dtype=np.longdouble)*np.pi/180.    # Inclination of the orbital plane in degrees
 
     x1 = np.array([0., -1., 0.],dtype=np.longdouble)*a[0]*(1.+e[0])
     x2 = np.array([0., 1., 0.],dtype=np.longdouble)*a[1]*(1.+e[1])
@@ -29,19 +29,20 @@ def main():
     v = np.array([v1, v2, v3],dtype=np.longdouble)
 
     #integration parameters
-    duration, step = 100*yr, np.array([1./(365.25*24.), 12./(365.25*24.), 24./(365.25*24.)],dtype=np.longdouble)*yr #integration time and step in seconds
+    duration, step = 5000.*yr, np.array([10./(365.25)],dtype=np.longdouble)*yr #integration time and step in seconds
     step = np.sort(step)[::-1]
     integrator = "leapfrog"
     n_bodies = 3
     display = False
-    savename = "{0:d}bodies_conc_{1:s}".format(n_bodies, integrator)
+    savename = "{0:d}bodies_psi45_{1:s}".format(n_bodies, integrator)
 
     bodies, bodysyst = [],[]
     for j in range(n_bodies):
         bodies.append(Body(m[j], q[j], v[j]))
     bin_syst = System(bodies[0:2])
     dyn_syst = System(bodies, main=True)
-    bodysyst = [[deepcopy(bin_syst), deepcopy(dyn_syst)] for _ in range(n_bodies)]
+    bodysyst = [[deepcopy(bin_syst), deepcopy(dyn_syst)] for _ in range(len(step))]
+
     #simulation start
     exe = ProcessPoolExecutor()
     future_ELae = []