switch back to unit time in seconds, prepare for parallel computing

lib/LeapFrog.py

@@ -58,9 +58,9 @@ def leapfrog(dyn_syst, bin_syst, duration, dt, recover_param=False, display=Fals
         if display and j % 10 == 0:
             # display progression
             if len(dyn_syst.bodylist) == 1:
-                d.on_running(dyn_syst, step=j, label="{0:.2f} years".format(j*dt))
+                d.on_running(dyn_syst, step=j, label="{0:.2f} years".format(j*dt/yr))
             else:
-                d.on_running(dyn_syst, step=j, label="{0:.2f} years".format(j*dt))
+                d.on_running(dyn_syst, step=j, label="{0:.2f} years".format(j*dt/yr))
     if display:
         d.close()
         if not savename is None:

lib/plots.py

@@ -110,12 +110,12 @@ def display_parameters(E,L,sma,ecc,parameters,savename=""):
     bodies = ""
     for body in dyn_syst.bodylist:
         bodies += str(body)+" ; "
-    title1, title2 = "Relative difference of the {0:s} ","for a system composed of {0:s}\n integrated with {1:s} for a duration of {2:.2f} years ".format(bodies, integrator, duration)
+    title1, title2 = "Relative difference of the {0:s} ","for a system composed of {0:s}\n integrated with {1:s} for a duration of {2:.2f} years ".format(bodies, integrator, duration/yr)
 
     fig1 = plt.figure(figsize=(15,7))
     ax1 = fig1.add_subplot(111)
     for i in range(len(E)):
-        ax1.plot(np.arange(E[i].shape[0]-1)*step[i], np.abs((E[i][1:]-E[i][0])/E[i][0]), label="step of {0:.2e}yr".format(step[i]))
+        ax1.plot(np.arange(E[i].shape[0]-1)*step[i]/yr, np.abs((E[i][1:]-E[i][0])/E[i][0]), label="step of {0:.2e}s".format(step[i]))
     ax1.set(xlabel=r"$t \, [yr]$", ylabel=r"$\left|\frac{\delta E_m}{E_m(t=0)}\right|$", yscale='log')
     ax1.legend()
     fig1.suptitle(title1.format("mechanical energy")+title2)
@@ -126,7 +126,7 @@ def display_parameters(E,L,sma,ecc,parameters,savename=""):
     for i in range(len(L)):
         dL = ((L[i]-L[i][0])/L[i][0])
         dL[np.isnan(dL)] = 0.
-        ax2.plot(np.arange(L[i].shape[0]-1)*step[i], np.abs(np.sum(dL[1:],axis=1)), label="step of {0:.2e}yr".format(step[i]))
+        ax2.plot(np.arange(L[i].shape[0]-1)*step[i]/yr, np.abs(np.sum(dL[1:],axis=1)), label="step of {0:.2e}s".format(step[i]))
     ax2.set(xlabel=r"$t \, [yr]$", ylabel=r"$\left|\frac{\delta \vec{L}}{\vec{L}(t=0)}\right|$",yscale='log')
     ax2.legend()
     fig2.suptitle(title1.format("kinetic moment")+title2)
@@ -134,8 +134,9 @@ def display_parameters(E,L,sma,ecc,parameters,savename=""):
 
     fig3 = plt.figure(figsize=(15,7))
     ax3 = fig3.add_subplot(111)
-    ax3.plot(np.arange(sma.shape[0])*step[i], sma, label="a (semi major axis)")
-    ax3.plot(np.arange(ecc.shape[0])*step[i], ecc, label="e (eccentricity)")
+    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.set(xlabel=r"$t \, [yr]$", ylabel=r"$a \, [au] \, or \, e$")
     ax3.legend()
     fig3.suptitle("Semi major axis and eccentricity "+title2)

lib/units.py

@@ -8,4 +8,4 @@ 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
 globals()['yr'] = 3.15576e7 #year in seconds
-globals()['Ga'] = G*Ms*yr**2/au**3 #Gravitational constant dimensionless
+globals()['Ga'] = G*Ms/au**3 #Gravitational constant dimensionless

main.py

@@ -27,9 +27,9 @@ def main():
     v = np.array([v1, v2, v3],dtype=np.longdouble)
 
     #integration parameters
-    duration, step = 100*yr/yr, np.array([1./(365.25*2.), 1./(365.25*1.), 5./(365.25*1.)],dtype=np.longdouble)*yr/yr #integration time and step in years
+    duration, step = 100*yr, np.array([60.],dtype=np.longdouble) #integration time and step in seconds
     step = np.sort(step)[::-1]
-    integrator = "hermite"
+    integrator = "leapfrog"
     n_bodies = 3
     display = False
     savename = "{0:d}bodies_{1:s}".format(n_bodies, integrator)

main_concurrent.py

@@ -0,0 +1,68 @@
+#!/usr/bin/python
+# -*- coding:utf-8 -*-
+from sys import exit as sysexit
+from copy import deepcopy
+import numpy as np
+import matplotlib.pyplot as plt
+from concurrent.futures import ThreadPoolExecutor, ProcessPoolExecutor
+from lib.objects import Body, System
+from lib.LeapFrog import leapfrog
+from lib.hermite import hermite
+from lib.plots import display_parameters
+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
+
+    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])
+    x3 = np.array([np.cos(psi[2]), 0., np.sin(psi[2])],dtype=np.longdouble)*a[2]*(1.+e[2])
+    q = np.array([x1, x2, x3],dtype=np.longdouble)
+
+    v1 = np.array([np.sqrt(Ga*m[0]*m[1]/((m[0]+m[1])*np.sqrt(np.sum((q[0]-q[1])**2)))), 0., 0.],dtype=np.longdouble)
+    v2 = np.array([-np.sqrt(Ga*m[0]*m[1]/((m[0]+m[1])*np.sqrt(np.sum((q[0]-q[1])**2)))), 0., 0.],dtype=np.longdouble)
+    v3 = np.array([0., np.sqrt(Ga*(m[0]+m[1])*(2./np.sqrt(np.sum(q[2]**2))-1./a[2])), 0.],dtype=np.longdouble)
+    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
+    step = np.sort(step)[::-1]
+    integrator = "leapfrog"
+    n_bodies = 3
+    display = False
+    savename = "{0:d}bodies_conc_{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)]
+    #simulation start
+    exe = ProcessPoolExecutor()
+    future_ELae = []
+    for i,step0 in enumerate(step):
+        if i != 0:
+            display = False
+        if integrator.lower() in ['leapfrog', 'frogleap', 'frog']:
+            future_ELae.append(exe.submit(leapfrog, bodysyst[i][1], bodysyst[i][0], duration, step0, recover_param=True, display=display, savename=savename))
+        elif integrator.lower() in ['hermite','herm']:
+            future_ELae.append(exe.submit(hermite, bodysyst[i][1], bodysyst[i][0], duration, step0, recover_param=True, display=display, savename=savename))
+    
+    E, L, sma, ecc = [], [], [], []
+    for future in future_ELae:
+        E0, L0, sma0, ecc0 = future.result()
+        E.append(E0)
+        L.append(L0)
+        sma.append(sma0)
+        ecc.append(ecc0)
+    parameters = [duration, step, dyn_syst, integrator]
+    display_parameters(E, L, sma, ecc, parameters=parameters, savename=savename)
+    return 0
+
+if __name__ == '__main__':
+    sysexit(main())