change sma and ecc computation to ECOM and LCOM

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
+        sma[j] = bin_syst.smaCOM
-        ecc[j] = bin_syst.ecc
+        ecc[j] = bin_syst.eccCOM
         phi[j] = dyn_syst.phi
 
         if display and j % 10 == 0:

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
+        sma[j] = bin_syst.smaCOM
-        ecc[j] = bin_syst.ecc
+        ecc[j] = bin_syst.eccCOM
         phi[j] = dyn_syst.phi
 
         if display and j % 10 == 0:

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.qb = body.q - COM
-            body.vb = body.vb - COMV
+            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
+        dr = self.bodylist[0].m/self.mu*self.bodylist[0].q#b
-        dv = self.bodylist[0].m/self.mu*self.bodylist[0].v
+        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
+    def ECOM(self): #return mechanical energy in the center of mass of a binary system
-        dr = self.bodylist[0].m/self.mu*self.bodylist[0].q
+        #self.COMShiftBin()
-        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
         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 :

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()

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__':