remove COMShift everywhere

lib/LeapFrog.py

@@ -27,7 +27,6 @@ def Kick(dyn_syst, dt):
 
 
 def LP(dyn_syst, dt):
-    dyn_syst.COMShift()
     Drift(dyn_syst, dt / 2)
     Kick(dyn_syst, dt)
     Drift(dyn_syst, dt / 2)
@@ -43,16 +42,23 @@ def leapfrog(dyn_syst, bin_syst, duration, dt, recover_param=False, display=Fals
         d.launch(dyn_syst.blackstyle)
 
     N = np.ceil(duration / dt).astype(int)
-    E = np.zeros(N,dtype=np.longdouble)
-    L = np.zeros((N, 3),dtype=np.longdouble)
-    sma = np.zeros(N,dtype=np.longdouble)
-    ecc = np.zeros(N,dtype=np.longdouble)
-    phi = np.zeros(N,dtype=np.longdouble)
-    for j in range(N):
+    E = np.zeros(N+1,dtype=np.longdouble)
+    L = np.zeros((N+1, 3),dtype=np.longdouble)
+    sma = np.zeros(N+1,dtype=np.longdouble)
+    ecc = np.zeros(N+1,dtype=np.longdouble)
+    phi = np.zeros(N+1,dtype=np.longdouble)
+
+    E[0] = dyn_syst.ECOM
+    L[0] = dyn_syst.LCOM
+    sma[0] = bin_syst.smaCOM
+    ecc[0] = bin_syst.eccCOM
+    phi[0] = dyn_syst.phi
+
+    for j in range(1,N+1):
         LP(dyn_syst,dt)
 
-        E[j] = dyn_syst.E
-        L[j] = dyn_syst.L
+        E[j] = dyn_syst.ECOM
+        L[j] = dyn_syst.LCOM
         sma[j] = bin_syst.smaCOM
         ecc[j] = bin_syst.eccCOM
         phi[j] = dyn_syst.phi

lib/hermite.py

@@ -69,7 +69,6 @@ def Correct(dyn_syst, dt):  # correct position and velocities of bodies in syste
 
 
 def HPC(dyn_syst, dt):  # update position and velocities of bodies in system with hermite predictor corrector
-    dyn_syst.COMShift()
     Update_a(dyn_syst)
     Update_j(dyn_syst)
     Predict(dyn_syst, dt)
@@ -89,17 +88,23 @@ def hermite(dyn_syst, bin_syst, duration, dt, recover_param=False, display=False
         d.launch(dyn_syst.blackstyle)
 
     N = np.ceil(duration / dt).astype(int)
-    E = np.zeros(N,dtype=np.longdouble)
-    L = np.zeros((N, 3),dtype=np.longdouble)
-    sma = np.zeros(N,dtype=np.longdouble)
-    ecc = np.zeros(N,dtype=np.longdouble)
-    phi = np.zeros(N,dtype=np.longdouble)
+    E = np.zeros(N+1,dtype=np.longdouble)
+    L = np.zeros((N+1, 3),dtype=np.longdouble)
+    sma = np.zeros(N+1,dtype=np.longdouble)
+    ecc = np.zeros(N+1,dtype=np.longdouble)
+    phi = np.zeros(N+1,dtype=np.longdouble)
 
-    for j in range(N):
+    E[0] = dyn_syst.ECOM
+    L[0] = dyn_syst.LCOM
+    sma[0] = bin_syst.smaCOM
+    ecc[0] = bin_syst.eccCOM
+    phi[0] = dyn_syst.phi
+
+    for j in range(1,N+1):
         HPC(dyn_syst, dt)
 
-        E[j] = dyn_syst.E
-        L[j] = dyn_syst.L
+        E[j] = dyn_syst.ECOM
+        L[j] = dyn_syst.LCOM
         sma[j] = bin_syst.smaCOM
         ecc[j] = bin_syst.eccCOM
         phi[j] = dyn_syst.phi

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/Ms, self.q, self.v)
+        return r"Body of mass: {0:.1e} $M_\odot$, position: {1}, velocity: {2}".format(self.m, self.q, self.v)
 
     def __str__(self): # Called upon "str(body)"
-        return r"Body of mass: {0:.1e} $M_\odot$".format(self.m/Ms)
+        return r"Body of mass: {0:.1e} $M_\odot$".format(self.m)
 
     @property
     def p(self):
@@ -67,19 +67,12 @@ class System(Body):
         zdata = np.array([body.q[2] for body in self.bodylist],dtype=np.longdouble)
         return xdata, ydata, zdata
 
-
-    def get_velocities(self): #return the positions of the bodies
-        vxdata = np.array([body.v[0] for body in self.bodylist],dtype=np.longdouble)
-        vydata = np.array([body.v[1] for body in self.bodylist],dtype=np.longdouble)
-        vzdata = np.array([body.v[2] for body in self.bodylist],dtype=np.longdouble)
-        return vxdata, vydata, vzdata
-
-
-    def get_momenta(self): #return the momenta of the bodies
-        pxdata = np.array([body.p[0] for body in self.bodylist],dtype=np.longdouble)
-        pydata = np.array([body.p[1] for body in self.bodylist],dtype=np.longdouble)
-        pzdata = np.array([body.p[2] for body in self.bodylist],dtype=np.longdouble)
-        return pxdata, pydata, pzdata
+    def get_positionsCOM(self): #return the positions of the bodies in the center of mass frame
+        COM = self.COM
+        xdata = np.array([body.q[0]-COM[0] for body in self.bodylist],dtype=np.longdouble)
+        ydata = np.array([body.q[1]-COM[1] for body in self.bodylist],dtype=np.longdouble)
+        zdata = np.array([body.q[2]-COM[2] for body in self.bodylist],dtype=np.longdouble)
+        return xdata, ydata, zdata
 
     @property
     def M(self): #return total system mass
@@ -149,29 +142,24 @@ class System(Body):
         return E
 
     @property
-    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#b
-        dv = self.bodylist[0].m/self.mu*self.bodylist[0].v#b
-        LCOM = self.mu*np.cross(dr,dv)
-
-        return LCOM
-
-    @property
-    def ECOM(self): #return mechanical energy in the center of mass of a binary system
-        #self.COMShiftBin()
-        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)
-
-        return ECOM
-
-    @property
-    def L(self): #return angular momentum of bodies in system
-        L = np.zeros(3,dtype=np.longdouble)
+    def ECOM(self): #return total energy of bodies in system in the center of mass frame
+        T, W = 0, 0
+        COM, COMV = self.COM, self.COMV
         for body in self.bodylist:
-            L = L + np.cross(body.q,body.p)
+            T = T + 1./2.*body.m*np.linalg.norm(body.v-COMV)**2
+            for otherbody in self.bodylist:
+                if body != otherbody:
+                    rij = np.linalg.norm(body.q-otherbody.q)
+                    W = W - Ga*body.m*otherbody.m/rij
+        E = T + W
+        return E
+
+    @property
+    def LCOM(self): #return angular momentum of bodies in system
+        L = np.zeros(3,dtype=np.longdouble)
+        COM, COMV = self.COM, self.COMV
+        for body in self.bodylist:
+            L = L + np.cross(body.q-COM,body.p-body.m*COMV)
         return L
 
     @property
@@ -187,6 +175,13 @@ class System(Body):
         E = T + W
         return E
 
+    @property
+    def L(self): #return angular momentum of bodies in system in the center of mass frame
+        L = np.zeros(3,dtype=np.longdouble)
+        for body in self.bodylist:
+            L = L + np.cross(body.q,body.p)
+        return L
+
     @property
     def eccCOM(self): #exentricity of two body sub system
         if len(self.bodylist) == 2 :

lib/plots.py

@@ -53,7 +53,7 @@ class DynamicUpdate():
 
         self.lines = []
         for i,body in enumerate(self.dyn_syst.bodylist):
-            x, y, z = body.q
+            x, y, z = body.q/au-self.dyn_syst.COM/au
             lines, = self.ax.plot([x],[y],[z],'o',color="C{0:d}".format(i),label="{0:s}".format(str(body)))
             self.lines.append(lines)
         self.lines = np.array(self.lines)
@@ -74,7 +74,7 @@ class DynamicUpdate():
             self.ax.set_zlabel('AU')
 
     def on_running(self, dyn_syst, step=None, label=None):
-        xdata, ydata, zdata = dyn_syst.get_positions()
+        xdata, ydata, zdata = dyn_syst.get_positionsCOM()
         values = np.sqrt(np.sum((np.array((xdata,ydata,zdata))**2).T,axis=1))/au
         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()
@@ -134,30 +134,40 @@ 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]/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$")
+    for i in range(len(E)):
+        ax3.plot(np.arange(E[i].shape[0])*step[-1]/yr, E[i], label="step of {0:.2e}s".format(step[i]))
+    ax3.set(xlabel=r"$t \, [yr]$", ylabel=r"$E \, [J]$")
     ax3.legend()
-    fig3.suptitle("Semi major axis and eccentricity "+title2)
-    fig3.savefig("plots/{0:s}a_e.png".format(savename),bbox_inches="tight")
-
+    fig3.suptitle("Mechanical energy of the whole system "+title2)
+    fig3.savefig("plots/{0:s}E.png".format(savename),bbox_inches="tight")
+    
     fig4 = plt.figure(figsize=(15,7))
     ax4 = fig4.add_subplot(111)
-    for i in range(len(E)):
-        ax4.plot(np.arange(E[i].shape[0])*step[-1]/yr, E[i], label="step of {0:.2e}s".format(step[i]))
-    ax4.set(xlabel=r"$t \, [yr]$", ylabel=r"$E \, [J]$")
+    for i in range(len(L)):
+        L2 = np.array([np.linalg.norm(Li)**2 for Li in L[i]])
+        ax4.plot(np.arange(L[i].shape[0])*step[i]/yr, L2, label=r"$L^2$ for step of {0:.2e}s".format(step[i]))
+    ax4.set(xlabel=r"$t \, [yr]$", ylabel=r"$\left|\vec{L}\right|^2 \, [kg^2 \cdot m^4 \cdot s^{-2}]$",yscale='log')
     ax4.legend()
-    fig4.suptitle("Mechanical energy of the whole system "+title2)
-    fig4.savefig("plots/{0:s}E.png".format(savename),bbox_inches="tight")
-    
+    fig4.suptitle("Squared norm of the kinetic moment of the whole system "+title2)
+    fig4.savefig("plots/{0:s}L.png".format(savename),bbox_inches="tight")
+
     fig5 = plt.figure(figsize=(15,7))
     ax5 = fig5.add_subplot(111)
-    for i in range(len(phi)):
-        ax5.plot(np.arange(phi[i].shape[0])*step[-1]/yr, phi[i], label="step of {0:.2e}s".format(step[i]))
-    ax5.set(xlabel=r"$t \, [yr]$", ylabel=r"$\phi \, [^{\circ}]$")
+    ax5.plot(np.arange(sma[-1].shape[0])*step[-1]/yr, sma[-1]/au, label="a (step of {0:.2e}s)".format(step[-1]))
+    ax5.plot(np.arange(ecc[-1].shape[0])*step[-1]/yr, ecc[-1], label="e (step of {0:.2e}s)".format(step[-1]))
+    ax5.set(xlabel=r"$t \, [yr]$", ylabel=r"$a \, [au] \, or \, e$")
     ax5.legend()
-    fig5.suptitle("Inclination angle of the perturbator's orbital plane "+title2)
-    fig5.savefig("plots/{0:s}phi.png".format(savename),bbox_inches="tight")
+    fig5.suptitle("Semi major axis and eccentricity "+title2)
+    fig5.savefig("plots/{0:s}a_e.png".format(savename),bbox_inches="tight")
+
+    fig6 = plt.figure(figsize=(15,7))
+    ax6 = fig6.add_subplot(111)
+    for i in range(len(phi)):
+        ax6.plot(np.arange(phi[i].shape[0])*step[-1]/yr, phi[i], label="step of {0:.2e}s".format(step[i]))
+    ax6.set(xlabel=r"$t \, [yr]$", ylabel=r"$\phi \, [^{\circ}]$")
+    ax6.legend()
+    fig6.suptitle("Inclination angle of the perturbator's orbital plane "+title2)
+    fig6.savefig("plots/{0:s}phi.png".format(savename),bbox_inches="tight")
     
     plt.show(block=True)
   

main.py

@@ -12,9 +12,9 @@ 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)/2.*au#/au   # Semi-major axis in astronomical units
+    a = np.array([1., 1., 10.],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., 60.],dtype=np.longdouble)*np.pi/180.    # Inclination of the orbital plane in degrees
+    psi = np.array([0., 0., 80.],dtype=np.longdouble)*np.pi/180.    # Inclination of the orbital plane in degrees
 
     x1 = np.array([0., -1., 0.],dtype=np.longdouble)*a[0]*(1.+e[0])
     x2 = np.array([0., 1., 0.],dtype=np.longdouble)*a[1]*(1.+e[1])
@@ -27,7 +27,7 @@ def main():
     v = np.array([v1, v2, v3],dtype=np.longdouble)
 
     #integration parameters
-    duration, step = 5000*yr, np.array([30.*86400.],dtype=np.longdouble) #integration time and step in seconds
+    duration, step = 500*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