diff --git a/lib/integrator.py b/lib/integrator.py
index c31076b..7a9c971 100755
--- a/lib/integrator.py
+++ b/lib/integrator.py
@@ -10,7 +10,9 @@ import time
 import numpy as np
 from lib.plots import DynamicUpdate
 
-globals()["G"] = 1. #Gravitationnal constant
+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
 
 def dp_dt(m_array, q_array):
     """
@@ -20,7 +22,7 @@ def dp_dt(m_array, q_array):
     dp_array = np.zeros(q_array.shape)
     for i in range(q_array.shape[0]):
         q_j = np.delete(q_array, i, 0)
-        m_j = np.delete(m_array, i, 0)#.reshape((q_j.shape[0],1))
+        m_j = np.delete(m_array, i, 0)
         dp_array[i] = -G*m_array[i]*np.sum(m_j/np.sum(np.sqrt(np.sum((q_j-q_array[i])**2, axis=0)))**3*(q_j-q_array[i]), axis=0)
     dp_array[np.isnan(dp_array)] = 0.
     return dp_array
diff --git a/lib/objects.py b/lib/objects.py
index ee535ec..42fdff7 100755
--- a/lib/objects.py
+++ b/lib/objects.py
@@ -5,6 +5,9 @@ Class definition for physical atribute
 """
 import numpy as np
 
+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:
 
@@ -67,7 +70,6 @@ class System:
         return L
 
     def Eval(self): #return total energy of bodies in system
-        G = 1. #Gravitational constant (here normalized)
         T = 0
         W = 0
         for body in self.bodylist:
diff --git a/lib/plots.py b/lib/plots.py
index f04f205..9a61177 100755
--- a/lib/plots.py
+++ b/lib/plots.py
@@ -30,17 +30,18 @@ class DynamicUpdate():
 
     def on_running(self, xdata, ydata, zdata, step=None, label=None):
         #Update data (with the new _and_ the old points)
-        self.lines.set_data_3d(xdata, ydata, zdata)
-        if not label is None:
-            self.ax.set_title(label)
-        #Need both of these in order to rescale
-        self.ax.relim()
-        self.ax.autoscale_view()
-        #We need to draw *and* flush
-        self.fig.canvas.draw()
-        self.fig.canvas.flush_events()
-        if not step is None and step%10==0:
-            self.fig.savefig("tmp/{0:05d}.png".format(step),bbox_inches="tight")
+        if not step is None and step%100==0:
+            self.lines.set_data_3d(xdata, ydata, zdata)
+            if not label is None:
+                self.ax.set_title(label)
+            #Need both of these in order to rescale
+            self.ax.relim()
+            self.ax.autoscale_view()
+            #We need to draw *and* flush
+            self.fig.canvas.draw()
+            self.fig.canvas.flush_events()
+            if not step is None and step%100==0:
+                self.fig.savefig("tmp/{0:05d}.png".format(step),bbox_inches="tight")
 
     #Example
     def __call__(self):
diff --git a/main.py b/main.py
index 036a05b..f3820aa 100755
--- a/main.py
+++ b/main.py
@@ -6,18 +6,24 @@ import matplotlib.pyplot as plt
 from lib.integrator import frogleap
 from lib.objects import Body, System
 
+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
+
 def main():
     #initialisation
-    m = np.array([1, 1, 1e-5])
+    m = np.array([1., 1., 0.1])*Ms  # Masses in Solar mass
+    a = np.array([1., 1., 5.])*au   # Semi-major axis in astronomical units
+    psi = np.array([0., 0., 80.])*np.pi/180.    # Inclination of the orbital plane in degrees
 
-    x1 = np.array([-1, 0, 0])
-    x2 = np.array([1, 0, 0])
-    x3 = np.array([100, 0, 20])
+    x1 = np.array([-1., 0., 0.])*a[0]
+    x2 = np.array([1., 0., 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([0, -0.35, 0])
-    v2 = np.array([0, 0.35, 0])
-    v3 = np.array([0, 20., 0])
+    v1 = np.array([0, -np.sqrt(G*Ms/np.sqrt(np.sum(x1**2))), 0])
+    v2 = np.array([0, np.sqrt(G*Ms/np.sqrt(np.sum(x2**2))), 0])
+    v3 = np.array([0, np.sqrt(G*Ms*(2./np.sqrt(np.sum(x3**2))-1./a[2])), 0])
     v = np.array([v1, v2, v3])
 
     bodylist = []
@@ -27,7 +33,7 @@ def main():
     dyn_syst.COMShift()
 
     duration, step = 100, 0.01
-    E, L = frogleap(duration, step, dyn_syst, recover_param=True)#, display=True)
+    E, L = frogleap(duration, step, dyn_syst, recover_param=True, display=True)
     fig1 = plt.figure(figsize=(30,15))
     ax1 = fig1.add_subplot(111)
     ax1.plot(np.arange(E.shape[0])/duration, E, label=r"$E_m$")
diff --git a/plots/Em.png b/plots/Em.png
index 495f36e..28a7d13 100644
Binary files a/plots/Em.png and b/plots/Em.png differ
diff --git a/plots/L2.png b/plots/L2.png
index 5417e7d..78e061f 100644
Binary files a/plots/L2.png and b/plots/L2.png differ
diff --git a/plots/dynsyst.gif b/plots/dynsyst.gif
index ed1772c..cd48ea4 100644
Binary files a/plots/dynsyst.gif and b/plots/dynsyst.gif differ