diff --git a/lib/objects.py b/lib/objects.py index d5ee015..eb7500e 100755 --- a/lib/objects.py +++ b/lib/objects.py @@ -24,10 +24,10 @@ class Body: self.vp = np.zeros(3,dtype=np.longdouble) def __repr__(self): # Called upon "print(body)" - return r"Body of mass: {0:.1e} $M_\odot$, position: {1}, velocity: {2}".format(self.m, self.q, self.v) + return r"Body of mass: {0:.1e} $M_\odot$, position: {1}, velocity: {2}".format(self.m/Ms, self.q, self.v) def __str__(self): # Called upon "str(body)" - return r"Body of mass: {0:.1e} $M_\odot$".format(self.m) + return r"Body of mass: {0:.1e} $M_\odot$".format(self.m/Ms) @property def p(self): @@ -150,7 +150,7 @@ class System(Body): for otherbody in self.bodylist: if body != otherbody: rij = np.linalg.norm(body.q-otherbody.q) - W = W - Ga*body.m*otherbody.m/rij + W = W - Ga*otherbody.m/(body.m+otherbody.m)*body.m**2/rij E = T + W return E @@ -171,7 +171,7 @@ class System(Body): for otherbody in self.bodylist: if body != otherbody: rij = np.linalg.norm(body.q-otherbody.q) - W = W - Ga*body.m*otherbody.m/rij + W = W - Ga*otherbody.m/(body.m+otherbody.m)*body.m**2/rij E = T + W return E @@ -185,7 +185,7 @@ class System(Body): @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. + 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 @@ -193,7 +193,7 @@ class System(Body): @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) + sma = -Ga*self.mu*self.bodylist[0].m/(2.*self.ECOM) else : sma = np.nan return sma @@ -201,7 +201,7 @@ class System(Body): @property def ecc(self): #exentricity of two body sub system if len(self.bodylist) == 2 : - ecc = (2.*self.EBIN*(np.linalg.norm(self.LBIN)**2))/((Ga**2)*(self.M**2)*(self.mu**3)) + 1. + ecc = (2.*self.EBIN*(np.linalg.norm(self.LBIN)**2))/(Ga**2*self.M**2*self.mu**3) + 1. else : ecc = np.nan return ecc @@ -209,7 +209,7 @@ class System(Body): @property def sma(self): #semi major axis of two body sub system if len(self.bodylist) == 2 : - sma = -Ga*self.M*self.mu/(2.*self.EBIN) + sma = -Ga*self.mu*self.bodylist[0].m/(2.*self.EBIN) else : sma = np.nan return sma diff --git a/main.py b/main.py index 5eaa253..e21ee86 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([0.75, 0.75, 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,7 +27,7 @@ def main(): v = np.array([v1, v2, v3],dtype=np.longdouble) #integration parameters - duration, step = 500*yr, np.array([30./1.*86400.],dtype=np.longdouble) #integration time and step in seconds + duration, step = 5000*yr, np.array([30./1.*86400.],dtype=np.longdouble) #integration time and step in seconds step = np.sort(step)[::-1] integrator = "leapfrog" n_bodies = 3 diff --git a/plots/2bodies_leapfrog_E.png b/plots/2bodies_leapfrog_E.png index 51c1aa2..113f151 100644 Binary files a/plots/2bodies_leapfrog_E.png and b/plots/2bodies_leapfrog_E.png differ diff --git a/plots/2bodies_leapfrog_L.png b/plots/2bodies_leapfrog_L.png new file mode 100644 index 0000000..e1307c2 Binary files /dev/null and b/plots/2bodies_leapfrog_L.png differ diff --git a/plots/2bodies_leapfrog_a_e.png b/plots/2bodies_leapfrog_a_e.png index a29656a..b46b79b 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 9d4002f..81a89b9 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 3f268a9..a29f6d8 100644 Binary files a/plots/2bodies_leapfrog_dL2.png and b/plots/2bodies_leapfrog_dL2.png differ diff --git a/plots/2bodies_leapfrog_phi.png b/plots/2bodies_leapfrog_phi.png index 3bdd8f0..1761349 100644 Binary files a/plots/2bodies_leapfrog_phi.png and b/plots/2bodies_leapfrog_phi.png differ diff --git a/plots/3bodies_leapfrog_E.png b/plots/3bodies_leapfrog_E.png index bddea11..bff8b6a 100644 Binary files a/plots/3bodies_leapfrog_E.png and b/plots/3bodies_leapfrog_E.png differ diff --git a/plots/3bodies_leapfrog_L.png b/plots/3bodies_leapfrog_L.png index 9d2acd6..9551d51 100644 Binary files a/plots/3bodies_leapfrog_L.png and b/plots/3bodies_leapfrog_L.png differ diff --git a/plots/3bodies_leapfrog_a_e.png b/plots/3bodies_leapfrog_a_e.png index 9bf16a0..2bf1293 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 788a489..617dd37 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 0e0900a..bcbc142 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 62c2ba5..de8178a 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 04e248e..4887124 100644 Binary files a/plots/3bodies_leapfrog_phi.png and b/plots/3bodies_leapfrog_phi.png differ