start adimension of problem

lib/objects.py

@@ -23,10 +23,10 @@ class Body:
         self.vp = np.zeros(3)
 
     def __repr__(self): # Called upon "print(body)"
-        return "Body of mass: {0:.2e}kg, position: {1}, velocity: {2}".format(self.m, self.q, self.v)
+        return r"Body of mass: {0:.2f} $M_\odot$, position: {1}, velocity: {2}".format(self.m, self.q, self.v)
 
     def __str__(self): # Called upon "str(body)"
-        return "Body of mass: {0:.2e}kg".format(self.m)
+        return r"Body of mass: {0:.2f} $M_\odot$".format(self.m)
 
 class System:
 
@@ -35,20 +35,30 @@ class System:
         self.bodylist = np.array(bodylist)
         self.time = 0
     
+    @property
     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
         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
-        return np.array([body.v for body in self.bodylist])
+        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
-        return np.array([body.p for body in self.bodylist])
+        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
 
     def Mass(self): #return total system mass
         mass = 0
@@ -133,7 +143,7 @@ class System:
             E[j] = self.Eval()
             L[j] = self.Lval()
 
-            if display and j%100==0:
+            if display and j%5==0:
                 # display progression
                 if len(self.bodylist) == 1:
                     d.on_running(self, step=j, label="step {0:d}/{1:d}".format(j,N))

lib/plots.py

@@ -10,8 +10,8 @@ from lib.units import *
 
 class DynamicUpdate():
     #Suppose we know the x range
-    min_x = -10
+    min_x = -1
-    max_x = 10
+    max_x = 1
 
     plt.ion()
 
@@ -64,11 +64,17 @@ class DynamicUpdate():
         self.ax.grid()
         if self.blackstyle:
             self.ax.legend(labelcolor='w', frameon=True, framealpha=0.2)
+            self.ax.set_xlabel('AU', color='w')
+            self.ax.set_ylabel('AU', color='w')
+            self.ax.set_zlabel('AU', color='w')
         else:
             self.ax.legend()
+            self.ax.set_xlabel('AU')
+            self.ax.set_ylabel('AU')
+            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_positions
         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()
@@ -88,25 +94,12 @@ class DynamicUpdate():
         #We need to draw *and* flush
         self.fig.canvas.draw()
         self.fig.canvas.flush_events()
-        if not step is None and step%1000==0:
+        if not step is None and step%10==0:
             self.fig.savefig("tmp/{0:06d}.png".format(step),bbox_inches="tight")
     
     def close(self):
         plt.close()
 
-    #Example
-    def __call__(self):
-        import numpy as np
-        import time
-        self.on_launch()
-        xdata = []
-        ydata = []
-        for x in np.arange(0,10,0.5):
-            xdata.append(x)
-            ydata.append(np.exp(-x**2)+10*np.exp(-(x-7)**2))
-            self.on_running(xdata, ydata)
-            time.sleep(1)
-        return xdata, ydata
 
 def display_parameters(E,L,parameters,savename=""):
     """
@@ -116,17 +109,18 @@ def display_parameters(E,L,parameters,savename=""):
     duration, step, dyn_syst, integrator = parameters
     if type(step) != list:
         step = [step]
-    if (len(E) == duration//step[0]) and (len(L) == duration//step[0]):
+    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)+" ; "
-    title = "Relative difference of the {0:s} "+"for a system composed of {0:s}\n integrated with {1:s} for a duration of {2:.2f} years ".format(bodies, integrator, duration/yr)
+    title = "Relative difference of the {0:s} "+"for a system composed of {0:s}\n integrated with {1:s} for a duration of {2:.2f} years ".format(bodies, integrator, duration)
 
     fig1 = plt.figure(figsize=(15,7))
     ax1 = fig1.add_subplot(111)
     for i in range(len(E)):
-        ax1.plot(np.arange(E[i].shape[0])*step[i]/yr, np.abs((E[i]-E[i][0])/E[i][0]), label="step of {0:.2e}yr".format(step[i]/yr))
+        ax1.plot(np.arange(E[i].shape[0])*step[i], np.abs((E[i]-E[i][0])/E[i][0]), label="step of {0:.2e}yr".format(step[i]))
     ax1.set(xlabel=r"$t (yr)$", ylabel=r"$\left|\frac{\delta E_m}{E_m(t=0)}\right|$", yscale='log')
     ax1.legend()
     fig1.suptitle(title.format("mechanical energy"))
@@ -137,7 +131,7 @@ def display_parameters(E,L,parameters,savename=""):
     for i in range(len(L)):
         dL = ((L[i]-L[i][0])/L[i][0])
         dL[np.isnan(dL)] = 0.
-        ax2.plot(np.arange(L[i].shape[0])*step[i]/yr, np.abs(np.sum(dL,axis=1)), label="step of {0:.2e}yr".format(step[i]/yr))
+        ax2.plot(np.arange(L[i].shape[0])*step[i], np.abs(np.sum(dL,axis=1)), label="step of {0:.2e}yr".format(step[i]))
     ax2.set(xlabel=r"$t (yr)$", ylabel=r"$\left|\frac{\delta \vec{L}}{\vec{L}(t=0)}\right|$",yscale='log')
     ax2.legend()
     fig2.suptitle(title.format("kinetic moment"))

lib/units.py

@@ -4,7 +4,7 @@
 Units used in the project.
 """
 
-globals()['G'] = 6.67e-11 #Gravitational constant in SI units
 globals()['Ms'] = 2e30 #Solar mass in kg
 globals()['au'] = 1.5e11 #Astronomical unit in m
 globals()['yr'] = 3.15576e7 #year in seconds
+globals()['G'] = 6.67e-11*yr**2 #Gravitational constant in SI units

main.py

@@ -9,8 +9,8 @@ from lib.units import *
 
 def main():
     #initialisation
-    m = np.array([1., 1., 1e-5])*Ms  # Masses in Solar mass
+    m = np.array([1., 1., 1e-5])*Ms/Ms  # Masses in Solar mass
-    a = np.array([1., 1., 5.])*au   # Semi-major axis in astronomical units
+    a = np.array([1., 1., 5.])*au/au   # Semi-major axis in astronomical units
     e = np.array([0., 0., 1./4.])   # Eccentricity
     psi = np.array([0., 0., 0.])*np.pi/180.    # Inclination of the orbital plane in degrees
 
@@ -25,7 +25,7 @@ def main():
     v = np.array([v1, v2, v3])
 
     #integration parameters
-    duration, step = 100*yr, [1e4, 1e5]
+    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