object = system

Properties

lib/objects.py

@@ -28,12 +28,15 @@ class Body:
     def __str__(self): # Called upon "str(body)"
         return r"Body of mass: {0:.2f} $M_\odot$".format(self.m)
 
-class System:
+class System(Body):
 
     def __init__(self, bodylist, blackstyle=True):
-        self.blackstyle = blackstyle
+        self.blackstyle = blackstyle #for dark mode in plot
         self.bodylist = np.array(bodylist)
-        self.time = 0
+        self.time = 0 #lifetime of system
+        self.m = self.M
+        self.q = self.COM
+        self.v = self.COMV
     
     @property
     def get_masses(self): #return the masses of each object
@@ -60,39 +63,43 @@ class System:
         pzdata = np.array([body.p[2] for body in self.bodylist])
         return pxdata, pydata, pzdata
 
-    def Mass(self): #return total system mass
+    @property
+    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
         coord = np.zeros(3)
         for body in self.bodylist:
             coord = coord + body.m*body.q
-        coord = coord/self.Mass()
+        coord = coord/self.M
         return coord
 
+    @property
     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.Mass()
+        coord = coord/self.M
         return coord
 
     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()
-        return 0
+            body.q = body.q - self.COM
+            body.p = body.p - self.COMV
 
-    def Lval(self): #return angular momentum of bodies in system
+    @property
+    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)
         return L
 
-    def Eval(self): #return total energy of bodies in system
+    @property
+    def E(self): #return total energy of bodies in system
         T = 0
         W = 0
         for body in self.bodylist:
@@ -101,7 +108,8 @@ class System:
                 if body != otherbody:
                     rij = np.linalg.norm(body.q-otherbody.q)
                     W = W - G*body.m*otherbody.m/rij
-        return T + W
+        E = T + W
+        return E
 
     def Drift(self, dt):
         for body in self.bodylist:
@@ -140,8 +148,8 @@ class System:
         for j in range(N):
             self.LP(dt)
 
-            E[j] = self.Eval()
-            L[j] = self.Lval()
+            E[j] = self.E
+            L[j] = self.L
 
             if display and j%5==0:
                 # display progression
@@ -237,8 +245,8 @@ class System:
         for j in range(N):
             self.HPC(dt)
 
-            E[j] = self.Eval()
-            L[j] = self.Lval()
+            E[j] = self.E
+            L[j] = self.L
 
             if display and j%100==0:
                 # display progression
@@ -261,3 +269,33 @@ class System:
     def __str__(self): # Called upon "str(system)"
         return str([str(body) for body in self.bodylist])
 
+    @property
+    def mu(self):
+        sum = 0
+        prod = 1
+        for body in self.bodylist:
+            prod = prod * body.m
+        mu = prod/self.M
+        return mu
+
+    @property
+    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.
+            return k
+
+    @property
+    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)
+            return sma
+
+
+
+
+
+

lib/plots.py

@@ -107,11 +107,6 @@ def display_parameters(E,L,parameters,savename=""):
     if savename != "":
         savename += "_"
     duration, step, dyn_syst, integrator = parameters
-    if type(step) != list:
-        step = [step]
-    print(E.shape, L.shape)
-    if (len(E.shape) == 1) and (len(L.shape) == 2):
-        E, L = [E], [L]
     bodies = ""
     for body in dyn_syst.bodylist:
         bodies += str(body)+" ; "

main.py

@@ -28,7 +28,7 @@ def main():
     duration, step = 100*yr/yr, np.array([1./(365.25*2.), 1./365.25])*yr/yr #integration time and step in years
     integrator = "leapfrog"
     n_bodies = 2
-    display = False
+    display = True
     savename = "{0:d}bodies_{1:s}".format(n_bodies, integrator)
 
     #simulation start