Modify main.py for single step length, concurrent.py to multiple steps

concurrent_steps.py

@@ -55,15 +55,16 @@ def main():
         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, gif=gif))
     
-    E, L, sma, ecc = [], [], [], []
+    E, L, sma, ecc, phi = [], [], [], [], []
     for future in future_ELae:
-        E0, L0, sma0, ecc0 = future.result()
+        E0, L0, sma0, ecc0, phi0 = future.result()
         E.append(E0)
         L.append(L0)
         sma.append(sma0)
         ecc.append(ecc0)
+        phi.append(phi0)
     parameters = [duration, step, dyn_syst, integrator]
-    display_parameters(E, L, sma, ecc, parameters=parameters, savename=savename)
+    display_parameters(E, L, sma, ecc, phi, parameters=parameters, savename=savename)
     return 0
 
 if __name__ == '__main__':

main.py

@@ -12,7 +12,7 @@ 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([0.75, 0.75, 5.],dtype=np.longdouble)*au#/au   # Semi-major axis in astronomical units
+    a = np.array([1., 1., 5.5],dtype=np.longdouble)*au#/au   # Semi-major axis in astronomical units
     e = np.array([0., 0., 0.25],dtype=np.longdouble)   # Eccentricity
     psi = np.array([0., 0., 80.],dtype=np.longdouble)*np.pi/180.    # Inclination of the orbital plane in degrees
 
@@ -27,8 +27,7 @@ def main():
     v = np.array([v1, v2, v3],dtype=np.longdouble)
 
     #integration parameters
-    duration, step = 5000*yr, np.array([30./1.*86400.],dtype=np.longdouble) #integration time and step in seconds
+    duration, step = 5000*yr, np.longdouble(1./2.*86400.) #integration time and step in seconds
-    step = np.sort(step)[::-1]
     integrator = "leapfrog"
     n_bodies = 3
     display = False
@@ -36,28 +35,19 @@ def main():
     savename = "{0:d}bodies_{1:s}".format(n_bodies, integrator)
 
     #simulation start
-    E, L, ecc, sma, phi = [], [], [], [], []
+    bodylist = []
-    for i,step0 in enumerate(step):
+    for j in range(n_bodies):
-        bodylist = []
+        bodylist.append(Body(m[j], q[j], v[j]))
-        for j in range(n_bodies):
+    bin_syst = System(bodylist[0:2])
-            bodylist.append(Body(m[j], q[j], v[j]))
+    dyn_syst = System(bodylist, main=True)
-        bin_syst = System(bodylist[0:2])
-        dyn_syst = System(bodylist, main=True)
 
-        if i != 0:
+    if integrator.lower() in ['leapfrog', 'frogleap', 'frog']:
-            display = False
+        E, L, sma, ecc, phi = leapfrog(dyn_syst, bin_syst, duration, step, recover_param=True, display=display, savename=savename, gif=gif)
-        if integrator.lower() in ['leapfrog', 'frogleap', 'frog']:
+    elif integrator.lower() in ['hermite','herm']:
-            E0, L0, sma0, ecc0, phi0 = leapfrog(dyn_syst, bin_syst, duration, step0, recover_param=True, display=display, savename=savename, gif=gif)
+        E, L, sma, ecc, phi = hermite(dyn_syst, bin_syst, duration, step, recover_param=True, display=display, savename=savename, gif=gif)
-        elif integrator.lower() in ['hermite','herm']:
-            E0, L0, sma0, ecc0, phi0 = hermite(dyn_syst, bin_syst, duration, step0, recover_param=True, display=display, savename=savename, gif=gif)
-        E.append(E0)
-        L.append(L0)
-        ecc.append(ecc0)
-        sma.append(sma0)
-        phi.append(phi0)
 
-    parameters = [duration, step, dyn_syst, integrator]
+    parameters = [duration, [step], dyn_syst, integrator]
-    display_parameters(E, L, sma, ecc, phi, parameters=parameters, savename=savename) #take the mean value of sma/ecc on given time interval (up to one period)
+    display_parameters([E], [L], [sma], [ecc], [phi], parameters=parameters, savename=savename) #take the mean value of sma/ecc on given time interval (up to one period)
     # np.convolve(sma, np.ones(int(period/step)))/int(period/step) -> moving average on the period duration
     return 0