diff --git a/lib/LeapFrog.py b/lib/LeapFrog.py index 1fc58d4..46a1046 100644 --- a/lib/LeapFrog.py +++ b/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 diff --git a/lib/objects.py b/lib/objects.py index 3d20332..aa516b7 100755 --- a/lib/objects.py +++ b/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 diff --git a/lib/plots.py b/lib/plots.py index 62b8ee6..b573d75 100755 --- a/lib/plots.py +++ b/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) diff --git a/lib/units.py b/lib/units.py index 642249f..667632f 100644 --- a/lib/units.py +++ b/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 \ No newline at end of file diff --git a/main.py b/main.py index 548003f..87ca220 100755 --- a/main.py +++ b/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__': diff --git a/plots/2bodies_leapfrog_a_e.png b/plots/2bodies_leapfrog_a_e.png index d477a79..0306913 100644 Binary files a/plots/2bodies_leapfrog_a_e.png and b/plots/2bodies_leapfrog_a_e.png differ diff --git a/plots/2bodies_leapfrog_dEm.png b/plots/2bodies_leapfrog_dEm.png index 309586d..d13cda3 100644 Binary files a/plots/2bodies_leapfrog_dEm.png and b/plots/2bodies_leapfrog_dEm.png differ diff --git a/plots/2bodies_leapfrog_dL2.png b/plots/2bodies_leapfrog_dL2.png index e9e06b7..096e0d1 100644 Binary files a/plots/2bodies_leapfrog_dL2.png and b/plots/2bodies_leapfrog_dL2.png differ