Update display with readable data

lib/objects.py

@@ -1,15 +1,12 @@
 #!/usr/bin/python
 # -*- coding:utf-8 -*-
 """
-Class definition for physical atribute
+Class definition for physical attribute
 """
 from os import system
 import numpy as np
 from lib.plots import DynamicUpdate
-
+from lib.units import *
-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
 
 class Body:
 
@@ -33,7 +30,8 @@ class Body:
 
 class System:
 
-    def __init__(self, bodylist):
+    def __init__(self, bodylist, blackstyle=True):
+        self.blackstyle = blackstyle
         self.bodylist = np.array(bodylist)
         self.time = 0
     
@@ -41,7 +39,10 @@ class System:
         return np.array([body.m for body in self.bodylist])
     
     def get_positions(self): #return the positions of the bodies
-        return np.array([body.q for body in self.bodylist])
+        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
     
     def get_velocities(self): #return the positions of the bodies
         return np.array([body.v for body in self.bodylist])
@@ -121,7 +122,7 @@ class System:
             except IOError:
                 system("rm tmp/*")
             d = DynamicUpdate(self)
-            d.on_launch()
+            d.launch(self.blackstyle)
 
         N = np.ceil(duration/dt).astype(int)
         E = np.zeros(N)
@@ -139,6 +140,7 @@ class System:
                 else:
                     d.on_running(self, step=j, label="step {0:d}/{1:d}".format(j,N))
         if display:
+            d.close()
             system("convert -delay 5 -loop 0 tmp/??????.png tmp/temp.gif && rm tmp/??????.png")
             system("convert tmp/temp.gif -fuzz 10% -layers Optimize plots/dynsyst.gif")# && rm tmp/temp.gif")
      
@@ -216,7 +218,7 @@ class System:
             except IOError:
                 system("rm tmp/*")
             d = DynamicUpdate(self)
-            d.on_launch()
+            d.launch(self.blackstyle)
 
         N = np.ceil(duration/dt).astype(int)
         E = np.zeros(N)
@@ -234,6 +236,7 @@ class System:
                 else:
                     d.on_running(self, step=j, label="step {0:d}/{1:d}".format(j,N))
         if display:
+            d.close()
             system("convert -delay 5 -loop 0 tmp/??????.png tmp/temp.gif && rm tmp/??????.png")
             system("convert tmp/temp.gif -fuzz 10% -layers Optimize plots/dynsyst.gif")# && rm tmp/temp.gif")
      

lib/plots.py

@@ -6,6 +6,7 @@ Implementation of the plotting and visualization functions.
 import numpy as np
 import time
 import matplotlib.pyplot as plt
+from lib.units import *
 
 class DynamicUpdate():
     #Suppose we know the x range
@@ -22,8 +23,7 @@ class DynamicUpdate():
         self.ax.set_ylim(factor*self.min_x, factor*self.max_x)
         self.ax.set_zlim(factor*self.min_x, factor*self.max_x)
     
-    def on_launch(self):
+    def set_blackstyle(self):
-        #Set up plot
         self.fig = plt.figure(figsize=(10,10), facecolor='k')
         self.ax = self.fig.add_subplot(projection='3d')
         self.ax.set_facecolor('k')
@@ -39,6 +39,18 @@ class DynamicUpdate():
         self.ax.w_xaxis.set_pane_color((0,0,0,0))
         self.ax.w_yaxis.set_pane_color((0,0,0,0))
         self.ax.w_zaxis.set_pane_color((0,0,0,0))
+        
+
+    def launch(self, blackstyle=True):
+        #Set up plot
+        if blackstyle:
+            self.blackstyle = True
+            self.set_blackstyle()
+        else:
+            self.blackstyle = False
+            self.fig = plt.figure(figsize=(10,10))
+            self.ax = self.fig.add_subplot(projection='3d')
+
         self.lines = []
         for i,body in enumerate(self.dyn_syst.bodylist):
             x, y, z = body.q
@@ -50,7 +62,10 @@ class DynamicUpdate():
         self.set_lims()
         #Other stuff
         self.ax.grid()
+        if self.blackstyle:
             self.ax.legend(labelcolor='w', frameon=True, framealpha=0.2)
+        else:
+            self.ax.legend()
 
     def on_running(self, dyn_syst, step=None, label=None):
         xdata, ydata, zdata = dyn_syst.get_positions()
@@ -62,7 +77,11 @@ class DynamicUpdate():
             x, y, z = body.q
             self.lines[i].set_data_3d([x], [y], [z])
         if not label is None:
+            if self.blackstyle:
                 self.ax.set_title(label,color='w')
+            else:
+                self.ax.set_title(label)
+
         #Need both of these in order to rescale
         self.ax.relim()
         self.ax.autoscale_view()
@@ -72,6 +91,9 @@ class DynamicUpdate():
         if not step is None and step%1000==0:
             self.fig.savefig("tmp/{0:06d}.png".format(step),bbox_inches="tight")
     
+    def close(self):
+        self.fig.close()
+
     #Example
     def __call__(self):
         import numpy as np
@@ -85,3 +107,32 @@ class DynamicUpdate():
             self.on_running(xdata, ydata)
             time.sleep(1)
         return xdata, ydata
+
+def display_parameters(E,L,parameters,savename=""):
+    """
+    """
+    if savename != "":
+        savename += "_"
+    duration, step, dyn_syst, integrator = parameters
+    bodies = ""
+    for body in dyn_syst.bodylist:
+        bodies += str(body)+" ; "
+    title = "Relative difference of the {} "+"for a system composed of {0:s} \nintegrated with {1:s} for a duration of {2:.2f} years with a step of {3:.2e} years.".format(bodies, integrator, duration/yr, step/yr)
+    fig1 = plt.figure(figsize=(15,7))
+    ax1 = fig1.add_subplot(111)
+    ax1.plot(np.arange(E.shape[0])*step/yr, np.abs((E-E[0])/E[0]), label=r"$\left|\frac{\delta E_m}{E_m(t=0)}\right|$")
+    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"))
+    fig1.savefig("plots/{0:s}dEm.png".format(savename),bbox_inches="tight")
+    fig2 = plt.figure(figsize=(15,7))
+    ax2 = fig2.add_subplot(111)
+    dL = ((L-L[0])/L[0])
+    dL[np.isnan(dL)] = 0.
+    ax2.plot(np.arange(L.shape[0])*step/yr, np.abs(np.sum(dL,axis=1)), label=r"$\left|\frac{\delta \vec{L}}{\vec{L}(t=0)}\right|$")
+    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"))
+    fig2.savefig("plots/{0:s}dL2.png".format(savename),bbox_inches="tight")
+    plt.show(block=True)
+  

lib/units.py

@@ -0,0 +1,10 @@
+#!/usr/bin/python
+# -*- coding:utf-8 -*-
+"""
+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

main.py

@@ -4,49 +4,37 @@ from sys import exit as sysexit
 import numpy as np
 import matplotlib.pyplot as plt
 from lib.objects import Body, System
-
+from lib.plots import display_parameters
-globals()['G'] = 6.67e-11 #Gravitational constant in SI units
+from lib.units import *
-globals()['Ms'] = 2e30 #Solar mass in kg
-globals()['au'] = 1.5e11 #Astronomical unit in m
 
 def main():
     #initialisation
-    m = np.array([1, 1, 0.1])*Ms  # Masses in Solar mass
+    m = np.array([1., 1., 0.])*Ms  # Masses in Solar mass
     a = np.array([1., 1., 5.])*au   # Semi-major axis in astronomical units
     e = np.array([0., 0., 1./4.])   # Eccentricity
-    psi = np.array([0., 0., 80.])*np.pi/180.    # Inclination of the orbital plane in degrees
+    psi = np.array([0., 0., 0.])*np.pi/180.    # Inclination of the orbital plane in degrees
 
     x1 = np.array([0., -1., 0.])*a[0]
     x2 = np.array([0., 1., 0.])*a[1]
     x3 = np.array([np.cos(psi[2]), 0., np.sin(psi[2])])*a[2]
     q = np.array([x1, x2, x3])
 
-    v1 = np.array([-np.sqrt(G*m[1]**2/((m[0]+m[1])*np.sqrt(np.sum((q[0]-q[1])**2)))), 0., 0.])
+    v1 = np.array([np.sqrt(G*m[1]**2/((m[0]+m[1])*np.sqrt(np.sum((q[0]-q[1])**2)))), 0., 0.])
-    v2 = np.array([np.sqrt(G*m[0]**2/((m[0]+m[1])*np.sqrt(np.sum((q[0]-q[1])**2)))), 0., 0.])
+    v2 = np.array([-np.sqrt(G*m[0]**2/((m[0]+m[1])*np.sqrt(np.sum((q[0]-q[1])**2)))), 0., 0.])
     v3 = np.array([0., np.sqrt(G*(m[0]+m[1])*(2./np.sqrt(np.sum(q[2]**2))-1./a[2])), 0.])
     v = np.array([v1, v2, v3])
 
     bodylist = []
-    for i in range(3):
+    for i in range(2):
         bodylist.append(Body(m[i], q[i], v[i]))
     dyn_syst = System(bodylist)
     dyn_syst.COMShift()
 
-    duration, step = 100*3e7, 1e4
+    duration, step = 100*yr, 1e4
-    #E, L = dyn_syst.leapfrog(duration, step, recover_param=True, display=True)
+    #E, L = dyn_syst.leapfrog(duration, step, recover_param=True)#, display=True)
-    E, L = dyn_syst.hermite(duration,step, recover_param=True, display=True)
+    E, L = dyn_syst.hermite(duration,step, recover_param=True)#, display=True)
-    plt.close()
+    parameters = [duration, step, dyn_syst, "hermite"]
-    fig1 = plt.figure(figsize=(30,15))
+    display_parameters(E, L, parameters=parameters, savename="2bodies_hermite")
-    ax1 = fig1.add_subplot(111)
-    ax1.plot(np.arange(E.shape[0])/duration, E, label=r"$E_m$")
-    ax1.legend()
-    fig1.savefig("plots/Em.png",bbox_inches="tight")
-    fig2 = plt.figure(figsize=(30,15))
-    ax2 = fig2.add_subplot(111)
-    ax2.plot(np.arange(L.shape[0])/duration, np.sum(L**2,axis=1), label=r"$L^2$")
-    ax2.legend()
-    fig2.savefig("plots/L2.png",bbox_inches="tight")
-    plt.show(block=True)
 
     return 0