diff --git a/lib/plots.py b/lib/plots.py index 09bd66c..ceb2968 100755 --- a/lib/plots.py +++ b/lib/plots.py @@ -114,22 +114,30 @@ def display_parameters(E,L,parameters,savename=""): if savename != "": savename += "_" duration, step, dyn_syst, integrator = parameters + if type(step) != list: + step = [step] + if (len(E) == duration//step[0]) and (len(L) == duration//step[0]): + E, L = [E], [L] bodies = "" for body in dyn_syst.bodylist: bodies += str(body)+" ; " - title = "Relative difference of the {} "+"for a system composed of {0:s}\n integrated with {1:s} for a duration of {2:.2f} years with a step of {3:.2e} years.".format(bodies, integrator, duration/yr, step/yr) + title = "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) - ax1.plot(np.arange(E.shape[0])*step/yr, np.abs((E-E[0])/E[0]), label=r"$\left|\frac{\delta E_m}{E_m(t=0)}\right|$") + for i in range(len(E)): + ax1.plot(np.arange(E[i].shape[0])*step[i]/yr, np.abs((E[i]-E[i][0])/E[i][0]), label="step of {0:.2e}yr".format(step[i]/yr)) ax1.set(xlabel=r"$t (yr)$", ylabel=r"$\left|\frac{\delta E_m}{E_m(t=0)}\right|$", yscale='log') ax1.legend() fig1.suptitle(title.format("mechanical energy")) fig1.savefig("plots/{0:s}dEm.png".format(savename),bbox_inches="tight") + fig2 = plt.figure(figsize=(15,7)) ax2 = fig2.add_subplot(111) - dL = ((L-L[0])/L[0]) - dL[np.isnan(dL)] = 0. - ax2.plot(np.arange(L.shape[0])*step/yr, np.abs(np.sum(dL,axis=1)), label=r"$\left|\frac{\delta \vec{L}}{\vec{L}(t=0)}\right|$") + 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])*step[i]/yr, np.abs(np.sum(dL,axis=1)), label="step of {0:.2e}yr".format(step[i]/yr)) ax2.set(xlabel=r"$t (yr)$", ylabel=r"$\left|\frac{\delta \vec{L}}{\vec{L}(t=0)}\right|$",yscale='log') ax2.legend() fig2.suptitle(title.format("kinetic moment")) diff --git a/main.py b/main.py index f471659..f63e80b 100755 --- a/main.py +++ b/main.py @@ -9,7 +9,7 @@ from lib.units import * def main(): #initialisation - m = np.array([1., 1., 0.])*Ms # Masses in Solar mass + m = np.array([1., 1., 0])*Ms # Masses in Solar mass a = np.array([1., 1., 5.])*au # Semi-major axis in astronomical units e = np.array([0., 0., 1./4.]) # Eccentricity psi = np.array([0., 0., 0.])*np.pi/180. # Inclination of the orbital plane in degrees @@ -30,11 +30,13 @@ def main(): dyn_syst = System(bodylist) dyn_syst.COMShift() - duration, step = 100*yr, 1e4 - #E, L = dyn_syst.leapfrog(duration, step, recover_param=True, display=True) - E, L = dyn_syst.hermite(duration,step, recover_param=True, display=True) - parameters = [duration, step, dyn_syst, "hermite"] - display_parameters(E, L, parameters=parameters, savename="3bodies_hermite") + duration, step1, step2 = 100*yr, 1e4, 1e5 + E1, L1 = dyn_syst.leapfrog(duration, step1, recover_param=True)#, display=True) + E2, L2 = dyn_syst.leapfrog(duration, step2, recover_param=True)#, display=True) + #E1, L1 = dyn_syst.hermite(duration, step1, recover_param=True)#, display=True) + #E2, L2 = dyn_syst.hermite(duration, step2, recover_param=True)#, display=True) + parameters = [duration, [step1, step2], dyn_syst, "leapfrog"] + display_parameters([E1, E2], [L1, L2], parameters=parameters, savename="3bodies_leapfrog") return 0 diff --git a/plots/2bodies_hermite_dEm.png b/plots/2bodies_hermite_dEm.png index f32d3fe..7852b4a 100644 Binary files a/plots/2bodies_hermite_dEm.png and b/plots/2bodies_hermite_dEm.png differ diff --git a/plots/2bodies_hermite_dL2.png b/plots/2bodies_hermite_dL2.png index a917c0d..5b49995 100644 Binary files a/plots/2bodies_hermite_dL2.png and b/plots/2bodies_hermite_dL2.png differ diff --git a/plots/2bodies_leapfrog_dEm.png b/plots/2bodies_leapfrog_dEm.png index 7135334..cb873e5 100644 Binary files a/plots/2bodies_leapfrog_dEm.png and b/plots/2bodies_leapfrog_dEm.png differ diff --git a/plots/2bodies_leapfrog_dL2.png b/plots/2bodies_leapfrog_dL2.png index a2f2629..57e4fc5 100644 Binary files a/plots/2bodies_leapfrog_dL2.png and b/plots/2bodies_leapfrog_dL2.png differ diff --git a/plots/3bodies_hermite_dEm.png b/plots/3bodies_hermite_dEm.png index 2fd1269..414048b 100644 Binary files a/plots/3bodies_hermite_dEm.png and b/plots/3bodies_hermite_dEm.png differ diff --git a/plots/3bodies_hermite_dL2.png b/plots/3bodies_hermite_dL2.png index 2f8c857..9d517a5 100644 Binary files a/plots/3bodies_hermite_dL2.png and b/plots/3bodies_hermite_dL2.png differ diff --git a/plots/3bodies_leapfrog_dEm.png b/plots/3bodies_leapfrog_dEm.png index 8a46925..e5ca97a 100644 Binary files a/plots/3bodies_leapfrog_dEm.png and b/plots/3bodies_leapfrog_dEm.png differ diff --git a/plots/3bodies_leapfrog_dL2.png b/plots/3bodies_leapfrog_dL2.png index b9f1b77..19c1ee1 100644 Binary files a/plots/3bodies_leapfrog_dL2.png and b/plots/3bodies_leapfrog_dL2.png differ