Finally get correct formula for Energy, Semi Major Axis and Eccentricity

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

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