get results with right Energy formula

lib/objects.py

@@ -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*otherbody.m/(body.m+otherbody.m)*body.m**2/rij
+                    W = W - Ga*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*otherbody.m/(body.m+otherbody.m)*body.m**2/rij
+                    W = W - Ga*otherbody.m*body.m/(2.*rij)
         E = T + W
         return E
 

main.py

@@ -12,7 +12,7 @@ from lib.units import *
 
 def main():
     #initialisation
-    m = np.array([1., 1., 0.],dtype=np.longdouble)*Ms#/Ms  # Masses in Solar mass
+    m = np.array([1., 1., 0.1],dtype=np.longdouble)*Ms#/Ms  # Masses in Solar mass
     a = np.array([1.00, 1.00, 10.0],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
@@ -28,12 +28,12 @@ def main():
     v = np.array([v1, v2, v3],dtype=np.longdouble)
 
     #integration parameters
-    duration, step = 5000*yr, np.longdouble(1./2.*86400.) #integration time and step in seconds
+    duration, step = 10000*yr, np.longdouble(1.0/2.*86400.) #integration time and step in seconds
     integrator = "leapfrog"
     n_bodies = 3
     display = False
     gif = False
-    savename = "{0:d}bodies_massless_{1:s}".format(n_bodies, integrator)
+    savename = "{0:d}bodies_1_{1:s}".format(n_bodies, integrator)
     display_param = True
 
     #simulation start