we did it !

lib/plots.py

@@ -134,11 +134,12 @@ def display_parameters(E,L,sma,ecc,phi,parameters,savename="",display_param=True
     Save integrated parameters plots to multiple png files.
     -----
     Inputs:
-    E, L, sma, ecc, phi : list of np.ndarray
+    E, L, sma, ecc, phi: list of np.ndarray
         list of integrated parameters value computed for each step length in
         parameters[step] list.
     parameters : list
-        list of simulation parameters : duration, steps, system, integrator.
+        list of simulation parameters : duration, steps, system, integrator,
+        and initialisation parameters.        
     savename : str
         default savename that will be prepend to each saved file path.
         Default to empty string.
@@ -148,11 +149,14 @@ def display_parameters(E,L,sma,ecc,phi,parameters,savename="",display_param=True
     """
     if savename != "":
         savename += "_"
-    duration, step, dyn_syst, integrator = parameters
+    duration, step, dyn_syst, integrator, init = parameters
+    a, e, psi = init
     bodies = ""
-    for body in dyn_syst.bodylist:
+    init_str = ""
+    for i, body in enumerate(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/yr)
+        init_str += r"a{0:d} = {1:.2f} au, e{0:d} = {2:.2f}, $\psi${0:d} = {3:.2f}° ; ".format(i+1,a[i]/au,e[i],psi[i]*180./np.pi)
+    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 with initial parameters\n {3:s}".format(bodies, integrator, duration/yr, init_str)
 
     fig1 = plt.figure(figsize=(15,7))
     ax1 = fig1.add_subplot(111)

main.py

@@ -12,8 +12,8 @@ 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([.75, .75, 5.],dtype=np.longdouble)*au#/au   # Semi-major axis in astronomical units
-    e = np.array([0., 0., 0.25],dtype=np.longdouble)   # Eccentricity
+    a = np.array([1.0, 1.0, 7.0],dtype=np.longdouble)*au#/au   # Semi-major axis in astronomical units
+    e = np.array([0., 0., 0.10],dtype=np.longdouble)   # Eccentricity
     psi = np.array([0., 0., 80.],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])
@@ -27,7 +27,7 @@ def main():
     v = np.array([v1, v2, v3],dtype=np.longdouble)
 
     #integration parameters
-    duration, step = 1000*yr, np.longdouble(1./2.*86400.) #integration time and step in seconds
+    duration, step = 5000*yr, np.longdouble(5./1.*86400.) #integration time and step in seconds
     integrator = "leapfrog"
     n_bodies = 3
     display = False
@@ -47,7 +47,7 @@ def main():
     elif integrator.lower() in ['hermite','herm']:
         E, L, sma, ecc, phi = hermite(dyn_syst, bin_syst, duration, step, recover_param=True, display=display, savename=savename, gif=gif)
 
-    parameters = [duration, [step], dyn_syst, integrator]
+    parameters = [duration, [step], dyn_syst, integrator, [a, e, psi]]
     display_parameters([E], [L], [sma], [ecc], [phi], parameters=parameters, savename=savename, display_param=display_param)
     return 0