Branch with adimensionned variables

lib/LeapFrog.py

@@ -22,7 +22,7 @@ def Kick(dyn_syst, dt):
         for otherbody in dyn_syst.bodylist:
             if body != otherbody:
                 rij = np.linalg.norm(body.q - otherbody.q)
-                body.a = body.a - (body.q - otherbody.q) * G * otherbody.m / (rij ** 3)
+                body.a = body.a - (body.q - otherbody.q) * Ga * otherbody.m / (rij ** 3)
         body.v = body.v + dt * body.a
 
 

lib/objects.py

@@ -22,10 +22,10 @@ class Body:
         self.vp = np.zeros(3)
 
     def __repr__(self): # Called upon "print(body)"
-        return r"Body of mass: {0:.1e} $M_\odot$, position: {1}, velocity: {2}".format(self.m/Ms, self.q, self.v)
+        return r"Body of mass: {0:.1e} $M_\odot$, position: {1}, velocity: {2}".format(self.m, self.q, self.v)
 
     def __str__(self): # Called upon "str(body)"
-        return r"Body of mass: {0:.1e} $M_\odot$".format(self.m/Ms)
+        return r"Body of mass: {0:.1e} $M_\odot$".format(self.m)
 
     @property
     def p(self):
@@ -118,7 +118,7 @@ class System(Body):
             for otherbody in self.bodylist:
                 if body != otherbody:
                     rij = np.linalg.norm(body.q-otherbody.q)
-                    W = W - G*body.m*otherbody.m/rij
+                    W = W - Ga*body.m*otherbody.m/rij
         E = T + W
         return E
          
@@ -133,7 +133,7 @@ class System(Body):
     @property
     def ecc(self): #exentricity of two body sub system
         if len(self.bodylist) == 2 :
-            ecc = (2.*self.E*(np.linalg.norm(self.L)**2))/((G**2)*(self.M**2)*(self.mu**3)) + 1.
+            ecc = (2.*self.E*(np.linalg.norm(self.L)**2))/((Ga**2)*(self.M**2)*(self.mu**3)) + 1.
         else :
             ecc = np.nan
         return ecc
@@ -141,7 +141,7 @@ class System(Body):
     @property
     def sma(self): #semi major axis of two body sub system
         if len(self.bodylist) == 2 :
-            sma = -G*self.M*self.mu/(2.*self.E)
+            sma = -Ga*self.M*self.mu/(2.*self.E)
         else :
             sma = np.nan
         return sma

lib/plots.py

@@ -110,12 +110,12 @@ def display_parameters(E,L,sma,ecc,parameters,savename=""):
     bodies = ""
     for body in dyn_syst.bodylist:
         bodies += str(body)+" ; "
-    title1, title2 = "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)
+    title1, title2 = "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(title1.format("mechanical energy")+title2)
@@ -126,7 +126,7 @@ def display_parameters(E,L,sma,ecc,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(title1.format("kinetic moment")+title2)
@@ -134,8 +134,8 @@ def display_parameters(E,L,sma,ecc,parameters,savename=""):
 
     fig3 = plt.figure(figsize=(15,7))
     ax3 = fig3.add_subplot(111)
-    ax3.plot(np.arange(sma.shape[0])*step[i]/yr, sma/au, label="a (semi major axis)")
-    ax3.plot(np.arange(ecc.shape[0])*step[i]/yr, ecc, label="e (eccentricity)")
+    ax3.plot(np.arange(sma.shape[0])*step[i], sma, label="a (semi major axis)")
+    ax3.plot(np.arange(ecc.shape[0])*step[i], ecc, label="e (eccentricity)")
     ax3.set(xlabel=r"$t \, [yr]$", ylabel=r"$a \, [au] \, or \, e$")
     ax3.legend()
     fig3.suptitle("Semi major axis and eccentricity "+title2)

lib/units.py

@@ -8,3 +8,4 @@ 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()['Ga'] = G*Ms*yr**2/au**3 #Gravitational constant adimensionned

main.py

@@ -11,8 +11,8 @@ from lib.units import *
 
 def main():
     #initialisation
-    m = np.array([1., 1., 1e-5])*Ms  # Masses in Solar mass
-    a = np.array([1., 1., 5.])*au   # Semi-major axis in astronomical units
+    m = np.array([1., 1., 1e-5])*Ms/Ms  # Masses in Solar mass
+    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
 
@@ -21,13 +21,13 @@ def main():
     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.])
-    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.])
+    v1 = np.array([np.sqrt(Ga*m[1]**2/((m[0]+m[1])*np.sqrt(np.sum((q[0]-q[1])**2)))), 0., 0.])
+    v2 = np.array([-np.sqrt(Ga*m[0]**2/((m[0]+m[1])*np.sqrt(np.sum((q[0]-q[1])**2)))), 0., 0.])
+    v3 = np.array([0., np.sqrt(Ga*(m[0]+m[1])*(2./np.sqrt(np.sum(q[2]**2))-1./a[2])), 0.])
     v = np.array([v1, v2, v3])
 
     #integration parameters
-    duration, step = 100*yr, np.array([1./(365.25*2.), 1./(365.25*1.), 5./(365.25*1.)])*yr #integration time and step in years
+    duration, step = 100*yr/yr, np.array([1./(365.25*2.), 1./(365.25*1.), 5./(365.25*1.)])*yr/yr #integration time and step in years
     step = np.sort(step)[::-1]
     integrator = "leapfrog"
     n_bodies = 2
@@ -52,7 +52,7 @@ def main():
 
     parameters = [duration, step, dyn_syst, integrator]
     display_parameters(E, L, sma, ecc, parameters=parameters, savename=savename)
-    print(sma/au)
+    print(sma,ecc)
     return 0
 
 if __name__ == '__main__':