diff --git a/.vscode/settings.json b/.vscode/settings.json index 3b66410..77acce4 100644 --- a/.vscode/settings.json +++ b/.vscode/settings.json @@ -1,3 +1,4 @@ { - "git.ignoreLimitWarning": true + "git.ignoreLimitWarning": true, + "python.pythonPath": "/usr/bin/python" } \ No newline at end of file diff --git a/lib/LeapFrog.py b/lib/LeapFrog.py index 7a6b423..0cb3623 100644 --- a/lib/LeapFrog.py +++ b/lib/LeapFrog.py @@ -21,18 +21,19 @@ def Kick(dyn_syst, dt): body.a = np.zeros(3,dtype=np.longdouble) 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.v = body.v + dt * body.a + rij = np.linalg.norm(body.q-otherbody.q) + body.a -= (body.q-otherbody.q)*G*otherbody.m/(rij**3) + body.v += dt*body.a def LP(dyn_syst, dt): - Drift(dyn_syst, dt / 2) + Drift(dyn_syst, dt/2) Kick(dyn_syst, dt) - Drift(dyn_syst, dt / 2) - dyn_syst.time = dyn_syst.time + dt + Drift(dyn_syst, dt/2) + dyn_syst.time = dyn_syst.time+dt -def leapfrog(dyn_syst, bin_syst, duration, dt, recover_param=False, display=False, savename=None, gif=False): +def leapfrog(dyn_syst, bin_syst, duration, dt, recover_param=False, + display=False, savename=None, gif=False): if display: try: system("mkdir tmp") diff --git a/lib/hermite.py b/lib/hermite.py index 994898e..465bcbc 100644 --- a/lib/hermite.py +++ b/lib/hermite.py @@ -17,7 +17,7 @@ def Update_a(dyn_syst): # update acceleration of bodies in system 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.q-otherbody.q)*G*otherbody.m/(rij**3) def Update_j(dyn_syst): # update jerk of bodies in system @@ -28,14 +28,14 @@ def Update_j(dyn_syst): # update jerk of bodies in system rij = np.linalg.norm(body.q - otherbody.q) deltav = (body.v - otherbody.v) deltar = (body.q - otherbody.q) - vr = deltav + 3. * deltar * np.inner(deltav, deltar) / (rij ** 2) - body.j = body.j - G * otherbody.m / (rij ** 3) * vr + vr = deltav+3.*deltar*np.inner(deltav,deltar)/(rij**2) + body.j -= G*otherbody.m/(rij**3)*vr -def Predict(dyn_syst, dt): # update predicted position and velocities of bodies in system +def Predict(dyn_syst, dt): # update predicted q and v of bodies in system for body in dyn_syst.bodylist: - body.qp = body.q + dt * body.v + ((dt ** 2) * body.a / 2.) + ((dt ** 3) * body.j / 6.) - body.vp = body.v + dt * body.a + ((dt ** 2) * body.j / 2.) + body.qp = body.q+dt*body.v+((dt**2)*body.a/2.)+((dt**3)*body.j/ 6.) + body.vp = body.v+dt*body.a+((dt**2)*body.j/2.) def Update_ap(dyn_syst): # update acceleration of bodies in system @@ -43,8 +43,8 @@ def Update_ap(dyn_syst): # update acceleration of bodies in system body.ap = np.zeros(3,dtype=np.longdouble) for otherbody in dyn_syst.bodylist: if body != otherbody: - rij = np.linalg.norm(body.qp - otherbody.qp) - body.ap = body.ap - (body.qp - otherbody.qp) * G * otherbody.m / (rij ** 3) + rij = np.linalg.norm(body.qp-otherbody.qp) + body.ap -= (body.qp-otherbody.qp)*G*otherbody.m/(rij**3) def Update_jp(dyn_syst): # update jerk of bodies in system @@ -52,20 +52,20 @@ def Update_jp(dyn_syst): # update jerk of bodies in system body.jp = np.zeros(3,dtype=np.longdouble) for otherbody in dyn_syst.bodylist: if body != otherbody: - rij = np.linalg.norm(body.qp - otherbody.qp) - deltav = (body.vp - otherbody.vp) - deltar = (body.qp - otherbody.qp) - vr = deltav + 3. * deltar * np.inner(deltav, deltar) / (rij ** 2) - body.jp = body.jp - G * otherbody.m / (rij ** 3) * vr + rij = np.linalg.norm(body.qp-otherbody.qp) + deltav = (body.vp-otherbody.vp) + deltar = (body.qp-otherbody.qp) + vr = deltav+3.*deltar*np.inner(deltav,deltar)/(rij**2) + body.jp -= G*otherbody.m/(rij**3)*vr def Correct(dyn_syst, dt): # correct position and velocities of bodies in system for body in dyn_syst.bodylist: - a2 = (6. * (body.a - body.ap) + dt * (4 * body.j + 2 * body.jp)) / (dt ** 2) - a3 = (12. * (body.a - body.ap) + dt * 6. * (body.j + body.jp)) / (dt ** 3) + a2 = (6.*(body.a-body.ap)+dt*(4.*body.j+2.*body.jp))/(dt**2) + a3 = (12.*(body.a-body.ap)+dt*6.*(body.j+body.jp))/(dt**3) - body.q = body.qp + ((dt ** 4) * a2 / 24.) + ((dt ** 5) * a3 / 120.) - body.v = body.vp + ((dt ** 3) * a2 / 6.) + ((dt ** 4) * a3 / 24.) + body.q = body.qp+((dt**4)*a2/24.)+((dt**5)*a3/120.) + body.v = body.vp+((dt**3)*a2/6.)+((dt**4)*a3/24.) def HPC(dyn_syst, dt): # update position and velocities of bodies in system with hermite predictor corrector @@ -78,7 +78,8 @@ def HPC(dyn_syst, dt): # update position and velocities of bodies in system wit dyn_syst.time = dyn_syst.time + dt -def hermite(dyn_syst, bin_syst, duration, dt, recover_param=False, display=False, savename=None, gif=False): +def hermite(dyn_syst, bin_syst, duration, dt, recover_param=False, + display=False, savename=None, gif=False): if display: try: system("mkdir tmp") @@ -87,7 +88,7 @@ def hermite(dyn_syst, bin_syst, duration, dt, recover_param=False, display=False d = DynamicUpdate(dyn_syst) d.launch(dyn_syst.blackstyle) - N = np.ceil(duration / dt).astype(int) + N = np.ceil(duration/dt).astype(int) E = np.zeros(N+1,dtype=np.longdouble) L = np.zeros((N+1, 3),dtype=np.longdouble) sma = np.zeros(N+1,dtype=np.longdouble) diff --git a/lib/objects.py b/lib/objects.py index 9f8517e..b18834e 100755 --- a/lib/objects.py +++ b/lib/objects.py @@ -14,8 +14,6 @@ class Body: self.m = np.longdouble(mass) self.q = np.longdouble(position) self.v = np.longdouble(velocity) - self.qb = np.longdouble(position) - self.vb = np.longdouble(velocity) self.a = np.zeros(3,dtype=np.longdouble) self.ap = np.zeros(3,dtype=np.longdouble) self.j = np.zeros(3,dtype=np.longdouble) @@ -24,7 +22,8 @@ class Body: self.vp = np.zeros(3,dtype=np.longdouble) 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/Ms, self.q, self.v) def __str__(self): # Called upon "str(body)" return r"Body of mass: {0:.1e} $M_\odot$".format(self.m/Ms) @@ -48,7 +47,6 @@ class System(Body): self.m = self.M self.q = self.COM self.v = self.COMV - self.coordarray = [] def __repr__(self): # Called upon "print(system)" return str([print(body) for body in self.bodylist]) @@ -58,20 +56,27 @@ class System(Body): def get_masses(self): #return the masses of each object - return np.array([body.m for body in self.bodylist],dtype=np.longdouble) + return np.array([body.m for body in self.bodylist], + dtype=np.longdouble) def get_positions(self): #return the positions of the bodies - xdata = np.array([body.q[0] for body in self.bodylist],dtype=np.longdouble) - ydata = np.array([body.q[1] for body in self.bodylist],dtype=np.longdouble) - zdata = np.array([body.q[2] for body in self.bodylist],dtype=np.longdouble) + xdata = np.array([body.q[0] for body in self.bodylist], + dtype=np.longdouble) + ydata = np.array([body.q[1] for body in self.bodylist], + dtype=np.longdouble) + zdata = np.array([body.q[2] for body in self.bodylist], + dtype=np.longdouble) return xdata, ydata, zdata - def get_positionsCOM(self): #return the positions of the bodies in the center of mass frame - COM = self.COM - xdata = np.array([body.q[0]-COM[0] for body in self.bodylist],dtype=np.longdouble) - ydata = np.array([body.q[1]-COM[1] for body in self.bodylist],dtype=np.longdouble) - zdata = np.array([body.q[2]-COM[2] for body in self.bodylist],dtype=np.longdouble) + def get_positionsCOM(self): #return the positions of the bodies + COM = self.COM # in the center of mass frame + xdata = np.array([body.q[0]-COM[0] for body in self.bodylist], + dtype=np.longdouble) + ydata = np.array([body.q[1]-COM[1] for body in self.bodylist], + dtype=np.longdouble) + zdata = np.array([body.q[2]-COM[2] for body in self.bodylist], + dtype=np.longdouble) return xdata, ydata, zdata @property @@ -105,45 +110,17 @@ class System(Body): coord = coord/self.M return coord - def COMShift(self): #Shift coordinates of bodies in system to COM frame and set COM at rest - COM = self.COM + def COMShift(self): #Shift coordinates of bodies in system to + COM = self.COM # COM frame and set COM at rest COMV = self.COMV for body in self.bodylist: body.q = body.q - COM body.v = body.v - COMV - def COMShiftBin(self): #Shift coordinates of inner binary system to COM frame and set COM at rest - COM = self.COM - COMV = self.COMV - for body in self.bodylist: - body.qb = body.q - COM - body.vb = body.v - COMV - @property - def LBIN(self): #return angular momentum of inner binary - self.COMShiftBin() - L = np.zeros(3,dtype=np.longdouble) - for body in self.bodylist: - L = L + np.cross(body.qb,body.pb) - return L - - @property - def EBIN(self): #return total energy of inner binary - self.COMShiftBin() + def ECOM(self): #return total energy in COM frame T = np.longdouble(0.) W = np.longdouble(0.) - for body in self.bodylist: - T = T + 1./2.*body.m*np.linalg.norm(body.vb)**2 - for otherbody in self.bodylist: - if body != otherbody: - rij = np.linalg.norm(body.qb-otherbody.qb) - W = W - G*body.m*otherbody.m/rij - E = T + W - return E - - @property - def ECOM(self): #return total energy of bodies in system in the center of mass frame - T, W = np.longdouble(0.), np.longdouble(0.) COM, COMV = self.COM, self.COMV for body in self.bodylist: T = T + 1./2.*body.m*np.linalg.norm(body.v-COMV)**2 @@ -155,37 +132,18 @@ class System(Body): return E @property - def LCOM(self): #return angular momentum of bodies in system + def LCOM(self): #return angular momentum of bodies in system in COM frame L = np.zeros(3,dtype=np.longdouble) COM, COMV = self.COM, self.COMV for body in self.bodylist: L = L + np.cross(body.q-COM,body.p-body.m*COMV) return L - @property - def E(self): #return total energy of bodies in system - T = np.longdouble(0.) - W = np.longdouble(0.) - for body in self.bodylist: - T = T + 1./2.*body.m*np.linalg.norm(body.v)**2 - for otherbody in self.bodylist: - if body != otherbody: - rij = np.linalg.norm(body.q-otherbody.q) - W = W - G*otherbody.m*body.m/(2.*rij) - E = T + W - return E - - @property - def L(self): #return angular momentum of bodies in system in the center of mass frame - L = np.zeros(3,dtype=np.longdouble) - for body in self.bodylist: - L = L + np.cross(body.q,body.p) - return L - @property def eccCOM(self): #exentricity of two body sub system if len(self.bodylist) == 2 : - ecc = (2.*self.ECOM*(np.linalg.norm(self.LCOM)**2))/(G**2*self.M**2*self.mu**3) + 1. + ecc = np.sqrt((2.*self.ECOM*(np.linalg.norm(self.LCOM)**2)) + /(G**2*self.M**2*self.mu**3) + 1.) else : ecc = np.nan @@ -200,30 +158,15 @@ class System(Body): return sma @property - def ecc(self): #exentricity of two body sub system - if len(self.bodylist) == 2 : - ecc = (2.*self.EBIN*(np.linalg.norm(self.LBIN)**2))/(G**2*self.M**2*self.mu**3) + 1. - else : - ecc = np.nan - return ecc - - @property - def sma(self): #semi major axis of two body sub system - if len(self.bodylist) == 2 : - sma = -G*self.mu*self.bodylist[0].m/(2.*self.EBIN) - else : - sma = np.nan - return sma - - @property - def phi(self): #return angle in degree between plans formed by perturbator plan and reference plan + def phi(self): #return angle formed by perturbator plan and reference plan if len(self.bodylist) == 3 : body1 = self.bodylist[0] body2 = self.bodylist[2] n1 = np.cross(body1.q-self.COM, body1.v-self.COMV) n2 = np.cross(body2.q-self.COM, body2.v-self.COMV) n1 = np.array([0., 0., 1.], dtype=np.longdouble) - phi = np.arccos(np.dot(n1, n2) / (np.linalg.norm(n1) * np.linalg.norm(n2)))*180./np.pi + phi = np.arccos(np.dot(n1, n2) / (np.linalg.norm(n1) + * np.linalg.norm(n2)))*180./np.pi else : phi = np.nan return phi diff --git a/lib/plots.py b/lib/plots.py index aed7b87..a9dbed6 100755 --- a/lib/plots.py +++ b/lib/plots.py @@ -82,7 +82,8 @@ class DynamicUpdate(): self.lines = [] for i,body in enumerate(self.dyn_syst.bodylist): x, y, z = body.q/au-self.dyn_syst.COM/au - lines, = self.ax.plot([x],[y],[z],'o',color="C{0:d}".format(i),label="{0:s}".format(str(body))) + lines, = self.ax.plot([x],[y],[z],'o',color="C{0:d}".format(i), + label="{0:s}".format(str(body))) self.lines.append(lines) self.lines = np.array(self.lines) #Autoscale on unknown axis and known lims on the other @@ -104,7 +105,8 @@ class DynamicUpdate(): def on_running(self, step=None, label=None): xdata, ydata, zdata = self.dyn_syst.get_positionsCOM() values = np.sqrt(np.sum((np.array((xdata,ydata,zdata))**2).T,axis=1))/au - 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.min_x = -np.max([np.abs(values).max(),self.max_x]) + self.max_x = np.max([np.abs(values).max(),self.max_x]) self.set_lims(factor=self.lim_factor) #Update data (with the new _and_ the old points) for i,body in enumerate(self.dyn_syst.bodylist): @@ -156,12 +158,15 @@ def display_parameters(E,L,sma,ecc,phi,parameters,savename="",display_param=True for i, body in enumerate(dyn_syst.bodylist): bodies += str(body)+" ; " init_str += r"a{0:d} = {1:.2f} au, e{0:d} = {2:.2f}, $\psi${0:d} = {3:.2f}° ; ".format(i+1,a[i]/au,e[i],psi[i]*180./np.pi) - 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 with initial parameters\n {3:s}".format(bodies, integrator, duration/yr, init_str) + title1 = "Relative difference of the {0:s} " + title2 = "for a system composed of {0:s}\n integrated with {1:s} for a duration of {2:.2f} years with initial parameters\n {3:s}".format(bodies, integrator, duration/yr, init_str) 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]-1)*step[i]/yr, np.abs((E[i][1:]-E[i][0])/E[i][0]), label="step of {0:.2e}s".format(step[i])) + ax1.plot(np.arange(E[i].shape[0]-1)*step[i]/yr, + np.abs((E[i][1:]-E[i][0])/E[i][0]), + label="step of {0:.2e}s".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) @@ -172,8 +177,10 @@ def display_parameters(E,L,sma,ecc,phi,parameters,savename="",display_param=True 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]-1)*step[i]/yr, np.abs(np.sum(dL[1:],axis=1)), label="step of {0:.2e}s".format(step[i])) - ax2.set(xlabel=r"$t \, [yr]$", ylabel=r"$\left|\frac{\delta \vec{L}}{\vec{L}(t=0)}\right|$",yscale='log') + ax2.plot(np.arange(L[i].shape[0]-1)*step[i]/yr, + np.abs(np.sum(dL[1:],axis=1)), + label="step of {0:.2e}s".format(step[i])) + ax2.set(xlabel=r"$t \, [yr]$", ylabel=r"$\left|\frac{\delta \vec{L}}{\vec{L}(t=0)}\right|$") ax2.legend() fig2.suptitle(title1.format("kinetic moment")+title2) fig2.savefig("plots/{0:s}dL2.png".format(savename),bbox_inches="tight") @@ -181,7 +188,8 @@ def display_parameters(E,L,sma,ecc,phi,parameters,savename="",display_param=True fig3 = plt.figure(figsize=(15,7)) ax3 = fig3.add_subplot(111) for i in range(len(E)): - ax3.plot(np.arange(E[i].shape[0])*step[-1]/yr, E[i], label="step of {0:.2e}s".format(step[i])) + ax3.plot(np.arange(E[i].shape[0])*step[-1]/yr, E[i], + label="step of {0:.2e}s".format(step[i])) ax3.set(xlabel=r"$t \, [yr]$", ylabel=r"$E \, [J]$") ax3.legend() fig3.suptitle("Mechanical energy of the whole system "+title2) @@ -191,16 +199,19 @@ def display_parameters(E,L,sma,ecc,phi,parameters,savename="",display_param=True ax4 = fig4.add_subplot(111) for i in range(len(L)): L2 = np.array([np.linalg.norm(Li)**2 for Li in L[i]]) - ax4.plot(np.arange(L[i].shape[0])*step[i]/yr, L2, label=r"$L^2$ for step of {0:.2e}s".format(step[i])) - ax4.set(xlabel=r"$t \, [yr]$", ylabel=r"$\left|\vec{L}\right|^2 \, [kg^2 \cdot m^4 \cdot s^{-2}]$",yscale='log') + ax4.plot(np.arange(L[i].shape[0])*step[i]/yr, L2, + label=r"$L^2$ for step of {0:.2e}s".format(step[i])) + ax4.set(xlabel=r"$t \, [yr]$", ylabel=r"$\left|\vec{L}\right|^2 \, [kg^2 \cdot m^4 \cdot s^{-2}]$") ax4.legend() fig4.suptitle("Squared norm of the kinetic moment of the whole system "+title2) fig4.savefig("plots/{0:s}L.png".format(savename),bbox_inches="tight") fig5 = plt.figure(figsize=(15,7)) ax5 = fig5.add_subplot(111) - ax5.plot(np.arange(sma[-1].shape[0])*step[-1]/yr, sma[-1]/au, label="a (step of {0:.2e}s)".format(step[-1])) - ax5.plot(np.arange(ecc[-1].shape[0])*step[-1]/yr, ecc[-1], label="e (step of {0:.2e}s)".format(step[-1])) + ax5.plot(np.arange(sma[-1].shape[0])*step[-1]/yr, sma[-1]/au, + label="a (step of {0:.2e}s)".format(step[-1])) + ax5.plot(np.arange(ecc[-1].shape[0])*step[-1]/yr, ecc[-1], + label="e (step of {0:.2e}s)".format(step[-1])) ax5.set(xlabel=r"$t \, [yr]$", ylabel=r"$a \, [au] \, or \, e$") ax5.legend() fig5.suptitle("Semi major axis and eccentricity "+title2) @@ -209,7 +220,8 @@ def display_parameters(E,L,sma,ecc,phi,parameters,savename="",display_param=True fig6 = plt.figure(figsize=(15,7)) ax6 = fig6.add_subplot(111) for i in range(len(phi)): - ax6.plot(np.arange(phi[i].shape[0])*step[-1]/yr, phi[i], label="step of {0:.2e}s".format(step[i])) + ax6.plot(np.arange(phi[i].shape[0])*step[-1]/yr, phi[i], + label="step of {0:.2e}s".format(step[i])) ax6.set(xlabel=r"$t \, [yr]$", ylabel=r"$\phi \, [^{\circ}]$") ax6.legend() fig6.suptitle("Inclination angle of the perturbator's orbital plane "+title2) diff --git a/lib/units.py b/lib/units.py index e23b866..642249f 100644 --- a/lib/units.py +++ b/lib/units.py @@ -7,4 +7,4 @@ 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 \ No newline at end of file +globals()['yr'] = 3.15576e7 #year in seconds diff --git a/main.py b/main.py index 4f8c251..d2c02d1 100755 --- a/main.py +++ b/main.py @@ -12,29 +12,40 @@ from lib.units import * def main(): #initialisation - m = np.array([1., 1., 0.1],dtype=np.longdouble)*Ms#/Ms # Masses in Solar mass - a = np.array([1.00, 1.00, 10.0],dtype=np.longdouble)*au#/au # Semi-major axis in astronomical units - e = np.array([0., 0., 0.25],dtype=np.longdouble) # Eccentricity - psi = np.array([0., 0., 80.],dtype=np.longdouble)*np.pi/180. # Inclination of the orbital plane in degrees + m = np.array([1., 1., 0.1], + dtype=np.longdouble)*Ms # Masses in Solar mass + a = np.array([1.00, 1.00, 7.0], + dtype=np.longdouble)*au # Semi-major axis in astronomical units + e = np.array([0., 0., 0.10], + dtype=np.longdouble) # Eccentricity + psi = np.array([0., 0., 35.], + dtype=np.longdouble)*np.pi/180. # Inclination of the orbital plane in degrees - x1 = np.array([0., -1., 0.],dtype=np.longdouble)*a[0]*(1.+e[0]) - x2 = np.array([0., 1., 0.],dtype=np.longdouble)*a[1]*(1.+e[1]) - x3 = np.array([np.cos(psi[2]), 0., np.sin(psi[2])],dtype=np.longdouble)*a[2]*(1.+e[2]) + x1 = a[0]*(1.+e[0])*np.array([0., -1., 0.], + dtype=np.longdouble) + x2 = a[1]*(1.+e[1])*np.array([0., 1., 0.], + dtype=np.longdouble) + x3 = a[2]*(1.+e[2])*np.array([np.cos(psi[2]), 0., np.sin(psi[2])], + dtype=np.longdouble) q = np.array([x1, x2, x3],dtype=np.longdouble) - v1 = np.array([np.sqrt(G*m[0]*m[1]/((m[0]+m[1])*np.sqrt(np.sum((q[0]-q[1])**2)))), 0., 0.],dtype=np.longdouble) - v2 = np.array([-np.sqrt(G*m[0]*m[1]/((m[0]+m[1])*np.sqrt(np.sum((q[0]-q[1])**2)))), 0., 0.],dtype=np.longdouble) - v3 = np.array([0., np.sqrt(G*(m[0]+m[1])*(2./np.sqrt(np.sum(q[2]**2))-1./a[2])), 0.],dtype=np.longdouble) + v1 = np.array([np.sqrt(G*m[0]*m[1]/((m[0]+m[1])*np.sqrt(np.sum((q[0]-q[1])**2)))), 0., 0.], + dtype=np.longdouble) + v2 = np.array([-np.sqrt(G*m[0]*m[1]/((m[0]+m[1])*np.sqrt(np.sum((q[0]-q[1])**2)))), 0., 0.], + dtype=np.longdouble) + v3 = np.array([0., np.sqrt(G*(m[0]+m[1])*(2./np.sqrt(np.sum(q[2]**2))-1./a[2])), 0.], + dtype=np.longdouble) v = np.array([v1, v2, v3],dtype=np.longdouble) #integration parameters - duration, step = 5000*yr, np.longdouble(1.0/1.*86400.) #integration time and step in seconds + duration = 2000*yr #integration time in seconds + step = np.longdouble(1.0/1.*86400.) #integration step in seconds integrator = "leapfrog" n_bodies = 3 display = False gif = False blackstyle = True - savename = "{0:d}bodies_{1:s}".format(n_bodies, integrator) + savename = "{0:d}bodies_test_{1:s}".format(n_bodies, integrator) display_param = True #simulation start @@ -47,15 +58,18 @@ def main(): t1 = time() if integrator.lower() in ['leapfrog', 'frogleap', 'frog']: - E, L, sma, ecc, phi = leapfrog(dyn_syst, bin_syst, duration, step, recover_param=True, display=display, savename=savename, gif=gif) + E, L, sma, ecc, phi = leapfrog(dyn_syst, bin_syst, duration, + step, recover_param=True, display=display, savename=savename, gif=gif) elif integrator.lower() in ['hermite','herm']: - E, L, sma, ecc, phi = hermite(dyn_syst, bin_syst, duration, step, recover_param=True, display=display, savename=savename, gif=gif) + E, L, sma, ecc, phi = hermite(dyn_syst, bin_syst, duration, + step, recover_param=True, display=display, savename=savename, gif=gif) t2=time() print("...Integration end.\n Elapsed time {0:.3f} sec".format(t2-t1)) parameters = [duration, [step], dyn_syst, integrator, [a, e, psi]] - display_parameters([E], [L], [sma], [ecc], [phi], parameters=parameters, savename=savename, display_param=display_param) + display_parameters([E], [L], [sma], [ecc], [phi], parameters=parameters, + savename=savename, display_param=display_param) return 0 if __name__ == '__main__': diff --git a/plots/2bodies_hermite_E.png b/plots/2bodies_hermite_E.png index aa2080f..3f566f7 100644 Binary files a/plots/2bodies_hermite_E.png and b/plots/2bodies_hermite_E.png differ diff --git a/plots/2bodies_hermite_L.png b/plots/2bodies_hermite_L.png index fcf577a..3ee9a24 100644 Binary files a/plots/2bodies_hermite_L.png and b/plots/2bodies_hermite_L.png differ diff --git a/plots/2bodies_hermite_a_e.png b/plots/2bodies_hermite_a_e.png index 8df4a68..dfe0120 100644 Binary files a/plots/2bodies_hermite_a_e.png and b/plots/2bodies_hermite_a_e.png differ diff --git a/plots/2bodies_hermite_dEm.png b/plots/2bodies_hermite_dEm.png index dc1f889..33a63c1 100644 Binary files a/plots/2bodies_hermite_dEm.png and b/plots/2bodies_hermite_dEm.png differ diff --git a/plots/2bodies_hermite_dL2.png b/plots/2bodies_hermite_dL2.png index b077d96..05f2e84 100644 Binary files a/plots/2bodies_hermite_dL2.png and b/plots/2bodies_hermite_dL2.png differ diff --git a/plots/2bodies_leapfrog_phi.png b/plots/2bodies_leapfrog_phi.png deleted file mode 100644 index e004769..0000000 Binary files a/plots/2bodies_leapfrog_phi.png and /dev/null differ