diff --git a/lib/LeapFrog.py b/lib/LeapFrog.py index cbe4ad0..67127da 100644 --- a/lib/LeapFrog.py +++ b/lib/LeapFrog.py @@ -53,8 +53,8 @@ def leapfrog(dyn_syst, bin_syst, duration, dt, recover_param=False, display=Fals E[j] = dyn_syst.E L[j] = dyn_syst.L - sma[j] = bin_syst.sma - ecc[j] = bin_syst.ecc + sma[j] = bin_syst.smaCOM + ecc[j] = bin_syst.eccCOM phi[j] = dyn_syst.phi if display and j % 10 == 0: diff --git a/lib/hermite.py b/lib/hermite.py index 3c6fbd4..49255d0 100644 --- a/lib/hermite.py +++ b/lib/hermite.py @@ -100,8 +100,8 @@ def hermite(dyn_syst, bin_syst, duration, dt, recover_param=False, display=False E[j] = dyn_syst.E L[j] = dyn_syst.L - sma[j] = bin_syst.sma - ecc[j] = bin_syst.ecc + sma[j] = bin_syst.smaCOM + ecc[j] = bin_syst.eccCOM phi[j] = dyn_syst.phi if display and j % 10 == 0: diff --git a/lib/objects.py b/lib/objects.py index 8e6083c..d7b6d43 100755 --- a/lib/objects.py +++ b/lib/objects.py @@ -123,8 +123,8 @@ class System(Body): COM = self.COM COMV = self.COMV for body in self.bodylist: - body.qb = body.qb - COM - body.vb = body.vb - COMV + body.qb = body.q - COM + body.vb = body.v - COMV @property def LBIN(self): #return angular momentum of inner binary @@ -149,24 +149,22 @@ class System(Body): return E @property - def LCOM(self): #return angular momentum of the center of mass + def LCOM(self): #return angular momentum in the center of mass of a binary system + #self.COMShiftBin() LCOM = np.zeros(3,dtype=np.longdouble) - dr = self.bodylist[0].m/self.mu*self.bodylist[0].q - dv = self.bodylist[0].m/self.mu*self.bodylist[0].v + dr = self.bodylist[0].m/self.mu*self.bodylist[0].q#b + dv = self.bodylist[0].m/self.mu*self.bodylist[0].v#b LCOM = self.mu*np.cross(dr,dv) - LCOM = self.L - return LCOM @property - def ECOM(self): #return mechanical energy of the center of mass - dr = self.bodylist[0].m/self.mu*self.bodylist[0].q - dv = self.bodylist[0].m/self.mu*self.bodylist[0].v + def ECOM(self): #return mechanical energy in the center of mass of a binary system + #self.COMShiftBin() + dr = self.bodylist[0].m/self.mu*self.bodylist[0].q#b + dv = self.bodylist[0].m/self.mu*self.bodylist[0].v#b ECOM = self.mu/2.*np.linalg.norm(dv)**2 - Ga*self.M*self.mu/np.linalg.norm(dr) - ECOM = self.E - return ECOM @property @@ -189,6 +187,22 @@ class System(Body): E = T + W return E + @property + def eccCOM(self): #exentricity of two body sub system + if len(self.bodylist) == 2 : + ecc = (2.*self.ECOM*(np.linalg.norm(self.LCOM)**2))/((Ga**2)*(self.M**2)*(self.mu**3)) + 1. + else : + ecc = np.nan + return ecc + + @property + def smaCOM(self): #semi major axis of two body sub system + if len(self.bodylist) == 2 : + sma = -Ga*self.M*self.mu/(2.*self.ECOM) + else : + sma = np.nan + return sma + @property def ecc(self): #exentricity of two body sub system if len(self.bodylist) == 2 : diff --git a/lib/plots.py b/lib/plots.py index 2ce0ed2..98dc150 100755 --- a/lib/plots.py +++ b/lib/plots.py @@ -134,7 +134,7 @@ def display_parameters(E,L,sma,ecc,phi,parameters,savename=""): fig3 = plt.figure(figsize=(15,7)) ax3 = fig3.add_subplot(111) - ax3.plot(np.arange(sma[-1].shape[0])*step[-1]/yr, sma[-1], label="a (step of {0:.2e}s)".format(step[-1])) + ax3.plot(np.arange(sma[-1].shape[0])*step[-1]/yr, sma[-1]/au, label="a (step of {0:.2e}s)".format(step[-1])) ax3.plot(np.arange(ecc[-1].shape[0])*step[-1]/yr, ecc[-1], label="e (step of {0:.2e}s)".format(step[-1])) ax3.set(xlabel=r"$t \, [yr]$", ylabel=r"$a \, [au] \, or \, e$") ax3.legend() diff --git a/main.py b/main.py index 022e231..7cd878d 100755 --- a/main.py +++ b/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([1., 1., 10.],dtype=np.longdouble)*au#/au # Semi-major axis in astronomical units + a = np.array([1., 1., 10.],dtype=np.longdouble)/2.*au#/au # Semi-major axis in astronomical units e = np.array([0., 0., 0.25],dtype=np.longdouble) # Eccentricity psi = np.array([0., 0., 60.],dtype=np.longdouble)*np.pi/180. # Inclination of the orbital plane in degrees @@ -27,7 +27,7 @@ def main(): v = np.array([v1, v2, v3],dtype=np.longdouble) #integration parameters - duration, step = 1000*yr, np.array([10.*86400.],dtype=np.longdouble) #integration time and step in seconds + duration, step = 5000*yr, np.array([30.*86400.],dtype=np.longdouble) #integration time and step in seconds step = np.sort(step)[::-1] integrator = "leapfrog" n_bodies = 3 @@ -57,7 +57,8 @@ def main(): phi.append(phi0) parameters = [duration, step, dyn_syst, integrator] - display_parameters(E, L, sma, ecc, phi, parameters=parameters, savename=savename) + 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 if __name__ == '__main__': diff --git a/plots/2bodies_hermite_E.png b/plots/2bodies_hermite_E.png new file mode 100644 index 0000000..2339f0a Binary files /dev/null and b/plots/2bodies_hermite_E.png differ diff --git a/plots/2bodies_hermite_a_e.png b/plots/2bodies_hermite_a_e.png index aa61dac..f54c803 100644 Binary files a/plots/2bodies_hermite_a_e.png and b/plots/2bodies_hermite_a_e.png differ diff --git a/plots/2bodies_hermite_dEm.png b/plots/2bodies_hermite_dEm.png index f2dcc01..5d024a2 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 1f1aaca..49e9105 100644 Binary files a/plots/2bodies_hermite_dL2.png and b/plots/2bodies_hermite_dL2.png differ diff --git a/plots/2bodies_hermite_dynsyst.gif b/plots/2bodies_hermite_dynsyst.gif index 16f3c60..48d7f39 100644 Binary files a/plots/2bodies_hermite_dynsyst.gif and b/plots/2bodies_hermite_dynsyst.gif differ diff --git a/plots/2bodies_leapfrog_E.png b/plots/2bodies_leapfrog_E.png new file mode 100644 index 0000000..b6738d8 Binary files /dev/null and b/plots/2bodies_leapfrog_E.png differ diff --git a/plots/2bodies_leapfrog_a_e.png b/plots/2bodies_leapfrog_a_e.png index bd08fa1..1798efc 100644 Binary files a/plots/2bodies_leapfrog_a_e.png and b/plots/2bodies_leapfrog_a_e.png differ diff --git a/plots/2bodies_leapfrog_dEm.png b/plots/2bodies_leapfrog_dEm.png index 904f7eb..5dee9fa 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 0e7359a..1c8ceb0 100644 Binary files a/plots/2bodies_leapfrog_dL2.png and b/plots/2bodies_leapfrog_dL2.png differ diff --git a/plots/2bodies_leapfrog_dynsyst.gif b/plots/2bodies_leapfrog_dynsyst.gif index 0807821..7a2dd63 100644 Binary files a/plots/2bodies_leapfrog_dynsyst.gif and b/plots/2bodies_leapfrog_dynsyst.gif differ diff --git a/plots/2bodies_leapfrog_phi.png b/plots/2bodies_leapfrog_phi.png new file mode 100644 index 0000000..e08a342 Binary files /dev/null and b/plots/2bodies_leapfrog_phi.png differ diff --git a/plots/3bodies_hermite_E.png b/plots/3bodies_hermite_E.png new file mode 100644 index 0000000..5753659 Binary files /dev/null and b/plots/3bodies_hermite_E.png differ diff --git a/plots/3bodies_hermite_a_e.png b/plots/3bodies_hermite_a_e.png index 1c79c2d..390d738 100644 Binary files a/plots/3bodies_hermite_a_e.png and b/plots/3bodies_hermite_a_e.png differ diff --git a/plots/3bodies_hermite_dEm.png b/plots/3bodies_hermite_dEm.png index d8a4fa6..c736699 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 631d303..0c60474 100644 Binary files a/plots/3bodies_hermite_dL2.png and b/plots/3bodies_hermite_dL2.png differ diff --git a/plots/3bodies_hermite_dynsyst.gif b/plots/3bodies_hermite_dynsyst.gif index d9d499d..1fb0178 100644 Binary files a/plots/3bodies_hermite_dynsyst.gif and b/plots/3bodies_hermite_dynsyst.gif differ diff --git a/plots/3bodies_hermite_phi.png b/plots/3bodies_hermite_phi.png new file mode 100644 index 0000000..187c25d Binary files /dev/null and b/plots/3bodies_hermite_phi.png differ diff --git a/plots/3bodies_leapfrog_E.png b/plots/3bodies_leapfrog_E.png index c601e6e..2c599a1 100644 Binary files a/plots/3bodies_leapfrog_E.png and b/plots/3bodies_leapfrog_E.png differ diff --git a/plots/3bodies_leapfrog_a_e.png b/plots/3bodies_leapfrog_a_e.png index 6314785..8fb684e 100644 Binary files a/plots/3bodies_leapfrog_a_e.png and b/plots/3bodies_leapfrog_a_e.png differ diff --git a/plots/3bodies_leapfrog_dEm.png b/plots/3bodies_leapfrog_dEm.png index d67398c..44a546c 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 3c8fe6d..509719c 100644 Binary files a/plots/3bodies_leapfrog_dL2.png and b/plots/3bodies_leapfrog_dL2.png differ diff --git a/plots/3bodies_leapfrog_dynsyst.gif b/plots/3bodies_leapfrog_dynsyst.gif index 7dba3c9..62c2ba5 100644 Binary files a/plots/3bodies_leapfrog_dynsyst.gif and b/plots/3bodies_leapfrog_dynsyst.gif differ diff --git a/plots/3bodies_leapfrog_phi.png b/plots/3bodies_leapfrog_phi.png index d7bc6f3..ad60f1e 100644 Binary files a/plots/3bodies_leapfrog_phi.png and b/plots/3bodies_leapfrog_phi.png differ