Merge remote-tracking branch 'origin/main' into main

# Conflicts:
#	lib/objects.py
#	plots/2bodies_leapfrog_dEm.png
#	plots/2bodies_leapfrog_dL2.png

lib/objects.py

@@ -8,7 +8,6 @@ import numpy as np
 from lib.plots import DynamicUpdate
 from lib.units import *
 
-
 class Body:
 
     def __init__(self, mass, position, velocity):
@@ -23,19 +22,18 @@ class Body:
         self.qp = np.zeros(3)
         self.vp = np.zeros(3)
 
-    def __repr__(self):  # Called upon "print(body)"
-        return r"Body of mass: {0:.2f} $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:.2f} $M_\odot$".format(self.m / Ms)
+    def __repr__(self): # Called upon "print(body)"
+        return r"Body of mass: {0:.2f} $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:.2f} $M_\odot$".format(self.m/Ms)
 
 class System(Body):
 
     def __init__(self, bodylist, blackstyle=True):
-        self.blackstyle = blackstyle  # for dark mode in plot
+        self.blackstyle = blackstyle #for dark mode in plot
         self.bodylist = np.array(bodylist)
-        self.time = 0  # lifetime of system
+        self.time = 0 #lifetime of system
         self.m = self.M
         self.q = self.COM
         self.v = self.COMV
@@ -47,99 +45,99 @@ class System(Body):
         return str([str(body) for body in self.bodylist])
 
     @property
-    def get_masses(self):  # return the masses of each object
+    def get_masses(self): #return the masses of each object
         return np.array([body.m for body in self.bodylist])
-
+    
     @property
-    def get_positions(self):  # return the positions of the bodies
+    def get_positions(self): #return the positions of the bodies
         xdata = np.array([body.q[0] for body in self.bodylist])
         ydata = np.array([body.q[1] for body in self.bodylist])
         zdata = np.array([body.q[2] for body in self.bodylist])
         return xdata, ydata, zdata
-
+    
     @property
-    def get_velocities(self):  # return the positions of the bodies
+    def get_velocities(self): #return the positions of the bodies
         vxdata = np.array([body.v[0] for body in self.bodylist])
         vydata = np.array([body.v[1] for body in self.bodylist])
         vzdata = np.array([body.v[2] for body in self.bodylist])
         return vxdata, vydata, vzdata
-
+    
     @property
-    def get_momenta(self):  # return the momenta of the bodies
+    def get_momenta(self): #return the momenta of the bodies
         pxdata = np.array([body.p[0] for body in self.bodylist])
         pydata = np.array([body.p[1] for body in self.bodylist])
         pzdata = np.array([body.p[2] for body in self.bodylist])
         return pxdata, pydata, pzdata
 
     @property
-    def M(self):  # return total system mass
+    def M(self): #return total system mass
         mass = 0
         for body in self.bodylist:
             mass = mass + body.m
         return mass
 
     @property
-    def COM(self):  # return center of mass in cartesian np_array
+    def COM(self): #return center of mass in cartesian np_array
         coord = np.zeros(3)
         for body in self.bodylist:
-            coord = coord + body.m * body.q
-        coord = coord / self.M
+            coord = coord + body.m*body.q
+        coord = coord/self.M
         return coord
 
     @property
-    def COMV(self):  # return center of mass velocity in cartesian np_array
+    def COMV(self): #return center of mass velocity in cartesian np_array
         coord = np.zeros(3)
         for body in self.bodylist:
             coord = coord + body.p
-        coord = coord / self.M
+        coord = coord/self.M
         return coord
 
-    def COMShift(self):  # Shift coordinates of bodies in system to COM frame and set COM at rest
+    def COMShift(self): #Shift coordinates of bodies in system to COM frame and set COM at rest
         for body in self.bodylist:
             body.q = body.q - self.COM
             body.p = body.p - self.COMV
 
     @property
-    def L(self):  # return angular momentum of bodies in system
+    def L(self): #return angular momentum of bodies in system
         L = np.zeros(3)
         for body in self.bodylist:
-            L = L + np.cross(body.q, body.p)
+            L = L + np.cross(body.q,body.p)
         return L
 
     @property
-    def E(self):  # return total energy of bodies in system
+    def E(self): #return total energy of bodies in system
         T = 0
         W = 0
         for body in self.bodylist:
-            T = T + 1. / 2. * body.m * np.linalg.norm(body.v) ** 2
+            T = T + 1./2.*body.m*np.linalg.norm(body.v)**2
             for otherbody in self.bodylist:
                 if body != otherbody:
-                    rij = np.linalg.norm(body.q - otherbody.q)
-                    W = W - G * body.m * otherbody.m / rij
+                    rij = np.linalg.norm(body.q-otherbody.q)
+                    W = W - G*body.m*otherbody.m/rij
         E = T + W
         return E
-
+         
     @property
     def mu(self):
         sum = 0
         prod = 1
         for body in self.bodylist:
             prod = prod * body.m
-        mu = prod / self.M
+        mu = prod/self.M
         return mu
 
     @property
-    def ex(self):  # exentricity of system (if composed of 2 bodies)
-        if len(self.bodylist) != 2:
+    def ex(self): #exentricity of system (if composed of 2 bodies)
+        if len(self.bodylist) != 2 :
             return np.nan
         else:
-            k = (2. * self.E * (np.linalg.norm(self.L) ** 2)) / ((G ** 2) * (self.M ** 2) * (self.mu ** 3)) + 1.
+            k = (2.*self.E*(np.linalg.norm(self.L)**2))/((G**2)*(self.M**2)*(self.mu**3)) + 1.
             return k
 
     @property
-    def sma(self):  # semi major axis of system (if composed of 2 bodies)
-        if len(self.bodylist) != 2:
+    def sma(self): #semi major axis of system (if composed of 2 bodies)
+        if len(self.bodylist) != 2 :
             return np.nan
         else:
-            sma = -G * self.M * self.mu / (2. * self.E)
+            sma = -G*self.M*self.mu/(2.*self.E)
             return sma

main.py

@@ -27,7 +27,7 @@ def main():
     v = np.array([v1, v2, v3])
 
     #integration parameters
-    duration, step = 100*yr, np.array([1./(365.25*2.), 1./365.25])*yr #integration time and step in years
+    duration, step = 100*yr, np.array([1./(365.25*4.), 1./(365.25*2.), 1./365.25])*yr #integration time and step in years
     integrator = "leapfrog"
     n_bodies = 2
     display = False