Correct computation of E and L in center of mass for binary system

lib/LeapFrog.py

@@ -55,7 +55,7 @@ def leapfrog(dyn_syst, bin_syst, duration, dt, recover_param=False, display=Fals
         sma[j] = bin_syst.sma
         ecc[j] = bin_syst.ecc
 
-        if display and j % 50 == 0:
+        if display and j % 10 == 0:
             # display progression
             if len(dyn_syst.bodylist) == 1:
                 d.on_running(dyn_syst, step=j, label="step {0:d}/{1:d}".format(j, N))

lib/hermite.py

@@ -102,7 +102,7 @@ def hermite(dyn_syst, duration, dt, recover_param=False, display=False, savename
         sma[j] = bin_syst.sma
         ecc[j] = bin_syst.ecc
 
-        if display and j % 100 == 0:
+        if display and j % 10 == 0:
             # display progression
             if len(dyn_syst.bodylist) == 1:
                 d.on_running(dyn_syst, step=j, label="step {0:d}/{1:d}".format(j, N))

lib/objects.py

@@ -81,6 +81,14 @@ class System(Body):
             mass = mass + body.m
         return mass
 
+    @property
+    def mu(self):
+        prod = 1
+        for body in self.bodylist:
+            prod = prod * body.m
+        mu = prod/self.M
+        return mu
+
     @property
     def COM(self): #return center of mass in cartesian np_array
         coord = np.zeros(3)
@@ -102,6 +110,31 @@ class System(Body):
             body.q = body.q - self.COM
             body.v = body.v - self.COMV
 
+    @property
+    def LCOM(self): #return angular momentum of the center of mass
+        LCOM = np.zeros(3)
+        dr = np.zeros(3)
+        dv = np.zeros(3)
+        for body in self.bodylist:
+            for otherbody in self.bodylist:
+                if body != otherbody:
+                    dr = body.q-otherbody.q
+                    dv = body.v-otherbody.v
+        LCOM = self.mu*np.cross(dr,dv)
+        return LCOM
+
+    @property
+    def ECOM(self): #return mechanical energy of the center of mass
+        dr = np.zeros(3)
+        dv = np.zeros(3)
+        for body in self.bodylist:
+            for otherbody in self.bodylist:
+                if body != otherbody:
+                    dr = body.q-otherbody.q
+                    dv = body.v-otherbody.v
+        ECOM = self.mu/2.*np.linalg.norm(dv)**2 - Ga*self.M*self.mu/np.linalg.norm(dr)
+        return ECOM
+
     @property
     def L(self): #return angular momentum of bodies in system
         L = np.zeros(3)
@@ -122,18 +155,10 @@ class System(Body):
         E = T + W
         return E
          
-    @property
-    def mu(self):
-        prod = 1
-        for body in self.bodylist:
-            prod = prod * body.m
-        mu = prod/self.M
-        return mu
-
     @property
     def ecc(self): #exentricity of two body sub system
         if len(self.bodylist) == 2 :
-            ecc = (2.*self.E*(np.linalg.norm(self.L)**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
@@ -141,7 +166,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.E)
+            sma = -Ga*self.M*self.mu/(2.*self.ECOM)
         else :
             sma = np.nan
         return sma

lib/plots.py

@@ -75,12 +75,12 @@ class DynamicUpdate():
 
     def on_running(self, dyn_syst, step=None, label=None):
         xdata, ydata, zdata = dyn_syst.get_positions
-        values = np.sqrt(np.sum((np.array((xdata,ydata,zdata))**2).T,axis=1))/au
+        values = np.sqrt(np.sum((np.array((xdata,ydata,zdata))**2).T,axis=1))
         self.min_x, self.max_x = -np.max([np.abs(values).max(),self.max_x]), np.max([np.abs(values).max(),self.max_x])
         self.set_lims()
         #Update data (with the new _and_ the old points)
         for i,body in enumerate(dyn_syst.bodylist):
-            x, y, z = body.q/au
+            x, y, z = body.q
             self.lines[i].set_data_3d([x], [y], [z])
         if not label is None:
             if self.blackstyle:
@@ -94,7 +94,7 @@ class DynamicUpdate():
         #We need to draw *and* flush
         self.fig.canvas.draw()
         self.fig.canvas.flush_events()
-        if not step is None and step%100==0:
+        if not step is None and step%50==0:
             self.fig.savefig("tmp/{0:06d}.png".format(step),bbox_inches="tight")
     
     def close(self):

main.py

@@ -31,7 +31,7 @@ def main():
     step = np.sort(step)[::-1]
     integrator = "leapfrog"
     n_bodies = 2
-    display = False
+    display = True
     savename = "{0:d}bodies_{1:s}".format(n_bodies, integrator)
 
     #simulation start
@@ -42,7 +42,9 @@ def main():
     dyn_syst = System(bodylist, main=True)
 
     E, L = [], []
-    for step0 in step:
+    for i,step0 in enumerate(step):
+        if i != 0:
+            display = False
         if integrator.lower() in ['leapfrog', 'frogleap', 'frog']:
             E0, L0, sma, ecc = leapfrog(dyn_syst, bin_syst, duration, step0, recover_param=True, display=display, savename=savename)
         elif integrator.lower() in ['hermite','herm']:
@@ -52,7 +54,6 @@ def main():
 
     parameters = [duration, step, dyn_syst, integrator]
     display_parameters(E, L, sma, ecc, parameters=parameters, savename=savename)
-    print(sma,ecc)
     return 0
 
 if __name__ == '__main__':