diff --git a/lib/LeapFrog.py b/lib/LeapFrog.py index 3d882f3..cbe4ad0 100644 --- a/lib/LeapFrog.py +++ b/lib/LeapFrog.py @@ -47,6 +47,7 @@ def leapfrog(dyn_syst, bin_syst, duration, dt, recover_param=False, display=Fals L = np.zeros((N, 3),dtype=np.longdouble) sma = np.zeros(N,dtype=np.longdouble) ecc = np.zeros(N,dtype=np.longdouble) + phi = np.zeros(N,dtype=np.longdouble) for j in range(N): LP(dyn_syst,dt) @@ -54,6 +55,7 @@ def leapfrog(dyn_syst, bin_syst, duration, dt, recover_param=False, display=Fals L[j] = dyn_syst.L sma[j] = bin_syst.sma ecc[j] = bin_syst.ecc + phi[j] = dyn_syst.phi if display and j % 10 == 0: if gif: @@ -72,4 +74,4 @@ def leapfrog(dyn_syst, bin_syst, duration, dt, recover_param=False, display=Fals system("convert tmp/temp.gif -fuzz 10% -layers Optimize plots/{0:s}_dynsyst.gif".format(savename)) if recover_param: - return E, L, sma, ecc \ No newline at end of file + return E, L, sma, ecc, phi \ No newline at end of file diff --git a/lib/hermite.py b/lib/hermite.py index 4ae1aa5..3c6fbd4 100644 --- a/lib/hermite.py +++ b/lib/hermite.py @@ -93,6 +93,7 @@ def hermite(dyn_syst, bin_syst, duration, dt, recover_param=False, display=False L = np.zeros((N, 3),dtype=np.longdouble) sma = np.zeros(N,dtype=np.longdouble) ecc = np.zeros(N,dtype=np.longdouble) + phi = np.zeros(N,dtype=np.longdouble) for j in range(N): HPC(dyn_syst, dt) @@ -101,6 +102,7 @@ def hermite(dyn_syst, bin_syst, duration, dt, recover_param=False, display=False L[j] = dyn_syst.L sma[j] = bin_syst.sma ecc[j] = bin_syst.ecc + phi[j] = dyn_syst.phi if display and j % 10 == 0: if gif: @@ -119,4 +121,4 @@ def hermite(dyn_syst, bin_syst, duration, dt, recover_param=False, display=False system("convert tmp/temp.gif -fuzz 10% -layers Optimize plots/{0:s}_dynsyst.gif".format(savename)) if recover_param: - return E, L, sma, ecc \ No newline at end of file + return E, L, sma, ecc, phi \ No newline at end of file diff --git a/lib/objects.py b/lib/objects.py index 56e02cd..8e6083c 100755 --- a/lib/objects.py +++ b/lib/objects.py @@ -5,8 +5,6 @@ Class definition for physical attribute """ from os import system import numpy as np -from astropy.coordinates import Angle -from astropy import units as u from lib.plots import DynamicUpdate from lib.units import * @@ -191,10 +189,6 @@ class System(Body): E = T + W return E - - - - @property def ecc(self): #exentricity of two body sub system if len(self.bodylist) == 2 : @@ -212,15 +206,13 @@ class System(Body): return sma @property - def phi(self,body1,body2): #return angle in degree between plans formed by body1 and body2 (perurbator) trajectories + def phi(self): #return angle in degree between plans formed by body1 and body2 (perurbator) trajectories if len(self.bodylist) == 3 : body1 = self.bodylist[0] body2 = self.bodylist[2] n1 = np.cross(body1.q, body1.v) n2 = np.cross(body2.q, body2.v) - phi = np.arccos(np.dot(n1, n2) / (np.linalg.norm(n1) * np.linalg.norm(n2))) - phi = Angle(phi, u.radian) - phi = phi.dec + phi = np.arccos(np.dot(n1, n2) / (np.linalg.norm(n1) * np.linalg.norm(n2)))*180./np.pi else : phi = np.nan return phi diff --git a/lib/plots.py b/lib/plots.py index bbc0dd3..2ce0ed2 100755 --- a/lib/plots.py +++ b/lib/plots.py @@ -101,7 +101,7 @@ class DynamicUpdate(): plt.close() -def display_parameters(E,L,sma,ecc,parameters,savename=""): +def display_parameters(E,L,sma,ecc,phi,parameters,savename=""): """ """ if savename != "": @@ -150,5 +150,14 @@ def display_parameters(E,L,sma,ecc,parameters,savename=""): fig4.suptitle("Mechanical energy of the whole system "+title2) fig4.savefig("plots/{0:s}E.png".format(savename),bbox_inches="tight") + fig5 = plt.figure(figsize=(15,7)) + ax5 = fig5.add_subplot(111) + for i in range(len(phi)): + ax5.plot(np.arange(phi[i].shape[0])*step[-1]/yr, phi[i], label="step of {0:.2e}s".format(step[i])) + ax5.set(xlabel=r"$t \, [yr]$", ylabel=r"$\phi \, [^{\circ}]$") + ax5.legend() + fig5.suptitle("Inclination angle of the perturbator's orbital plane "+title2) + fig5.savefig("plots/{0:s}phi.png".format(savename),bbox_inches="tight") + plt.show(block=True) \ No newline at end of file diff --git a/main.py b/main.py index 4ce22e6..022e231 100755 --- a/main.py +++ b/main.py @@ -31,12 +31,12 @@ def main(): step = np.sort(step)[::-1] integrator = "leapfrog" n_bodies = 3 - display = True + display = False gif = False savename = "{0:d}bodies_{1:s}".format(n_bodies, integrator) #simulation start - E, L, ecc, sma = [], [], [], [] + E, L, ecc, sma, phi = [], [], [], [], [] for i,step0 in enumerate(step): bodylist = [] for j in range(n_bodies): @@ -47,16 +47,17 @@ def main(): if i != 0: display = False if integrator.lower() in ['leapfrog', 'frogleap', 'frog']: - E0, L0, sma0, ecc0 = leapfrog(dyn_syst, bin_syst, duration, step0, recover_param=True, display=display, savename=savename, gif=gif) + E0, L0, sma0, ecc0, phi0 = leapfrog(dyn_syst, bin_syst, duration, step0, recover_param=True, display=display, savename=savename, gif=gif) elif integrator.lower() in ['hermite','herm']: - E0, L0, sma0, ecc0 = hermite(dyn_syst, bin_syst, duration, step0, recover_param=True, display=display, savename=savename, gif=gif) + 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] - 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__': diff --git a/plots/3bodies_leapfrog_phi.png b/plots/3bodies_leapfrog_phi.png new file mode 100644 index 0000000..d7bc6f3 Binary files /dev/null and b/plots/3bodies_leapfrog_phi.png differ