Tibeuleu/KozaiLidov
1
0
Fork 0
Code Activity

tidy things up

Browse Source
This commit is contained in:
Thibault Barnouin
2022-01-18 22:33:33 +01:00
parent e10b4a380c
commit 5168a52ade
13 changed files with 111 additions and 139 deletions
Download Patch File Download Diff File Expand all files Collapse all files

3
.vscode/settings.json vendored
View File

@@ -1,3 +1,4 @@
{ {
"git.ignoreLimitWarning": true "git.ignoreLimitWarning": true,
"python.pythonPath": "/usr/bin/python"
} }

15
lib/LeapFrog.py
View File

@@ -21,18 +21,19 @@ def Kick(dyn_syst, dt):
body.a = np.zeros(3,dtype=np.longdouble) body.a = np.zeros(3,dtype=np.longdouble)
for otherbody in dyn_syst.bodylist: for otherbody in dyn_syst.bodylist:
if body != otherbody: if body != otherbody:
rij = np.linalg.norm(body.q - otherbody.q) 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)
body.v = body.v + dt * body.a body.v += dt*body.a
def LP(dyn_syst, dt): def LP(dyn_syst, dt):
Drift(dyn_syst, dt / 2) Drift(dyn_syst, dt/2)
Kick(dyn_syst, dt) Kick(dyn_syst, dt)
Drift(dyn_syst, dt / 2) Drift(dyn_syst, dt/2)
dyn_syst.time = dyn_syst.time + dt 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: if display:
try: try:
system("mkdir tmp") system("mkdir tmp")

39
lib/hermite.py
View File

@@ -17,7 +17,7 @@ def Update_a(dyn_syst): # update acceleration of bodies in system
for otherbody in dyn_syst.bodylist: for otherbody in dyn_syst.bodylist:
if body != otherbody: if body != otherbody:
rij = np.linalg.norm(body.q - otherbody.q) 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 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) rij = np.linalg.norm(body.q - otherbody.q)
deltav = (body.v - otherbody.v) deltav = (body.v - otherbody.v)
deltar = (body.q - otherbody.q) deltar = (body.q - otherbody.q)
vr = deltav + 3. * deltar * np.inner(deltav, deltar) / (rij ** 2) vr = deltav+3.*deltar*np.inner(deltav,deltar)/(rij**2)
body.j = body.j - G * otherbody.m / (rij ** 3) * vr 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: for body in dyn_syst.bodylist:
body.qp = body.q + dt * body.v + ((dt ** 2) * body.a / 2.) + ((dt ** 3) * body.j / 6.) 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.vp = body.v+dt*body.a+((dt**2)*body.j/2.)
def Update_ap(dyn_syst): # update acceleration of bodies in system 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) body.ap = np.zeros(3,dtype=np.longdouble)
for otherbody in dyn_syst.bodylist: for otherbody in dyn_syst.bodylist:
if body != otherbody: if body != otherbody:
rij = np.linalg.norm(body.qp - otherbody.qp) rij = np.linalg.norm(body.qp-otherbody.qp)
body.ap = body.ap - (body.qp - otherbody.qp) * G * otherbody.m / (rij ** 3) body.ap -= (body.qp-otherbody.qp)*G*otherbody.m/(rij**3)
def Update_jp(dyn_syst): # update jerk of bodies in system 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) body.jp = np.zeros(3,dtype=np.longdouble)
for otherbody in dyn_syst.bodylist: for otherbody in dyn_syst.bodylist:
if body != otherbody: if body != otherbody:
rij = np.linalg.norm(body.qp - otherbody.qp) rij = np.linalg.norm(body.qp-otherbody.qp)
deltav = (body.vp - otherbody.vp) deltav = (body.vp-otherbody.vp)
deltar = (body.qp - otherbody.qp) deltar = (body.qp-otherbody.qp)
vr = deltav + 3. * deltar * np.inner(deltav, deltar) / (rij ** 2) vr = deltav+3.*deltar*np.inner(deltav,deltar)/(rij**2)
body.jp = body.jp - G * otherbody.m / (rij ** 3) * vr body.jp -= G*otherbody.m/(rij**3)*vr
def Correct(dyn_syst, dt): # correct position and velocities of bodies in system def Correct(dyn_syst, dt): # correct position and velocities of bodies in system
for body in dyn_syst.bodylist: for body in dyn_syst.bodylist:
a2 = (6. * (body.a - body.ap) + dt * (4 * body.j + 2 * body.jp)) / (dt ** 2) 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) 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.q = body.qp+((dt**4)*a2/24.)+((dt**5)*a3/120.)
body.v = body.vp + ((dt ** 3) * a2 / 6.) + ((dt ** 4) * a3 / 24.) 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 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 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: if display:
try: try:
system("mkdir tmp") 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 = DynamicUpdate(dyn_syst)
d.launch(dyn_syst.blackstyle) 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) E = np.zeros(N+1,dtype=np.longdouble)
L = np.zeros((N+1, 3),dtype=np.longdouble) L = np.zeros((N+1, 3),dtype=np.longdouble)
sma = np.zeros(N+1,dtype=np.longdouble) sma = np.zeros(N+1,dtype=np.longdouble)

111
lib/objects.py
View File

@@ -14,8 +14,6 @@ class Body:
self.m = np.longdouble(mass) self.m = np.longdouble(mass)
self.q = np.longdouble(position) self.q = np.longdouble(position)
self.v = np.longdouble(velocity) self.v = np.longdouble(velocity)
self.qb = np.longdouble(position)
self.vb = np.longdouble(velocity)
self.a = np.zeros(3,dtype=np.longdouble) self.a = np.zeros(3,dtype=np.longdouble)
self.ap = np.zeros(3,dtype=np.longdouble) self.ap = np.zeros(3,dtype=np.longdouble)
self.j = 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) self.vp = np.zeros(3,dtype=np.longdouble)
def __repr__(self): # Called upon "print(body)" 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)" 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/Ms)
@@ -48,7 +47,6 @@ class System(Body):
self.m = self.M self.m = self.M
self.q = self.COM self.q = self.COM
self.v = self.COMV self.v = self.COMV
self.coordarray = []
def __repr__(self): # Called upon "print(system)" def __repr__(self): # Called upon "print(system)"
return str([print(body) for body in self.bodylist]) 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 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 def get_positions(self): #return the positions of the bodies
xdata = np.array([body.q[0] for body in self.bodylist],dtype=np.longdouble) xdata = np.array([body.q[0] for body in self.bodylist],
ydata = np.array([body.q[1] for body in self.bodylist],dtype=np.longdouble) dtype=np.longdouble)
zdata = np.array([body.q[2] 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 return xdata, ydata, zdata
def get_positionsCOM(self): #return the positions of the bodies in the center of mass frame def get_positionsCOM(self): #return the positions of the bodies
COM = self.COM 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) xdata = np.array([body.q[0]-COM[0] for body in self.bodylist],
ydata = np.array([body.q[1]-COM[1] for body in self.bodylist],dtype=np.longdouble) dtype=np.longdouble)
zdata = np.array([body.q[2]-COM[2] 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 return xdata, ydata, zdata
@property @property
@@ -105,45 +110,17 @@ class System(Body):
coord = coord/self.M coord = coord/self.M
return coord return coord
def COMShift(self): #Shift coordinates of bodies in system to COM frame and set COM at rest def COMShift(self): #Shift coordinates of bodies in system to
COM = self.COM COM = self.COM # COM frame and set COM at rest
COMV = self.COMV COMV = self.COMV
for body in self.bodylist: for body in self.bodylist:
body.q = body.q - COM body.q = body.q - COM
body.v = body.v - COMV 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 @property
def LBIN(self): #return angular momentum of inner binary def ECOM(self): #return total energy in COM frame
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()
T = np.longdouble(0.) T = np.longdouble(0.)
W = 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 COM, COMV = self.COM, self.COMV
for body in self.bodylist: for body in self.bodylist:
T = T + 1./2.*body.m*np.linalg.norm(body.v-COMV)**2 T = T + 1./2.*body.m*np.linalg.norm(body.v-COMV)**2
@@ -155,37 +132,18 @@ class System(Body):
return E return E
@property @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) L = np.zeros(3,dtype=np.longdouble)
COM, COMV = self.COM, self.COMV COM, COMV = self.COM, self.COMV
for body in self.bodylist: for body in self.bodylist:
L = L + np.cross(body.q-COM,body.p-body.m*COMV) L = L + np.cross(body.q-COM,body.p-body.m*COMV)
return L 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 @property
def eccCOM(self): #exentricity of two body sub system def eccCOM(self): #exentricity of two body sub system
if len(self.bodylist) == 2 : 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 : else :
ecc = np.nan ecc = np.nan
@@ -200,30 +158,15 @@ class System(Body):
return sma return sma
@property @property
def ecc(self): #exentricity of two body sub system def phi(self): #return angle formed by perturbator plan and reference plan
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
if len(self.bodylist) == 3 : if len(self.bodylist) == 3 :
body1 = self.bodylist[0] body1 = self.bodylist[0]
body2 = self.bodylist[2] body2 = self.bodylist[2]
n1 = np.cross(body1.q-self.COM, body1.v-self.COMV) n1 = np.cross(body1.q-self.COM, body1.v-self.COMV)
n2 = np.cross(body2.q-self.COM, body2.v-self.COMV) n2 = np.cross(body2.q-self.COM, body2.v-self.COMV)
n1 = np.array([0., 0., 1.], dtype=np.longdouble) 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 : else :
phi = np.nan phi = np.nan
return phi return phi

36
lib/plots.py
View File

@@ -82,7 +82,8 @@ class DynamicUpdate():
self.lines = [] self.lines = []
for i,body in enumerate(self.dyn_syst.bodylist): for i,body in enumerate(self.dyn_syst.bodylist):
x, y, z = body.q/au-self.dyn_syst.COM/au 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.append(lines)
self.lines = np.array(self.lines) self.lines = np.array(self.lines)
#Autoscale on unknown axis and known lims on the other #Autoscale on unknown axis and known lims on the other
@@ -104,7 +105,8 @@ class DynamicUpdate():
def on_running(self, step=None, label=None): def on_running(self, step=None, label=None):
xdata, ydata, zdata = self.dyn_syst.get_positionsCOM() xdata, ydata, zdata = self.dyn_syst.get_positionsCOM()
values = np.sqrt(np.sum((np.array((xdata,ydata,zdata))**2).T,axis=1))/au 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) self.set_lims(factor=self.lim_factor)
#Update data (with the new _and_ the old points) #Update data (with the new _and_ the old points)
for i,body in enumerate(self.dyn_syst.bodylist): 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): for i, body in enumerate(dyn_syst.bodylist):
bodies += str(body)+" ; " 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) 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)) fig1 = plt.figure(figsize=(15,7))
ax1 = fig1.add_subplot(111) ax1 = fig1.add_subplot(111)
for i in range(len(E)): 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.set(xlabel=r"$t \, [yr]$", ylabel=r"$\left|\frac{\delta E_m}{E_m(t=0)}\right|$", yscale='log')
ax1.legend() ax1.legend()
fig1.suptitle(title1.format("mechanical energy")+title2) 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)): for i in range(len(L)):
dL = ((L[i]-L[i][0])/L[i][0]) dL = ((L[i]-L[i][0])/L[i][0])
dL[np.isnan(dL)] = 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.plot(np.arange(L[i].shape[0]-1)*step[i]/yr,
ax2.set(xlabel=r"$t \, [yr]$", ylabel=r"$\left|\frac{\delta \vec{L}}{\vec{L}(t=0)}\right|$",yscale='log') 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() ax2.legend()
fig2.suptitle(title1.format("kinetic moment")+title2) fig2.suptitle(title1.format("kinetic moment")+title2)
fig2.savefig("plots/{0:s}dL2.png".format(savename),bbox_inches="tight") 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)) fig3 = plt.figure(figsize=(15,7))
ax3 = fig3.add_subplot(111) ax3 = fig3.add_subplot(111)
for i in range(len(E)): 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.set(xlabel=r"$t \, [yr]$", ylabel=r"$E \, [J]$")
ax3.legend() ax3.legend()
fig3.suptitle("Mechanical energy of the whole system "+title2) 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) ax4 = fig4.add_subplot(111)
for i in range(len(L)): for i in range(len(L)):
L2 = np.array([np.linalg.norm(Li)**2 for Li in L[i]]) 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.plot(np.arange(L[i].shape[0])*step[i]/yr, L2,
ax4.set(xlabel=r"$t \, [yr]$", ylabel=r"$\left|\vec{L}\right|^2 \, [kg^2 \cdot m^4 \cdot s^{-2}]$",yscale='log') 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() ax4.legend()
fig4.suptitle("Squared norm of the kinetic moment of the whole system "+title2) 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") fig4.savefig("plots/{0:s}L.png".format(savename),bbox_inches="tight")
fig5 = plt.figure(figsize=(15,7)) fig5 = plt.figure(figsize=(15,7))
ax5 = fig5.add_subplot(111) 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(sma[-1].shape[0])*step[-1]/yr, sma[-1]/au,
ax5.plot(np.arange(ecc[-1].shape[0])*step[-1]/yr, ecc[-1], label="e (step of {0:.2e}s)".format(step[-1])) 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.set(xlabel=r"$t \, [yr]$", ylabel=r"$a \, [au] \, or \, e$")
ax5.legend() ax5.legend()
fig5.suptitle("Semi major axis and eccentricity "+title2) 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)) fig6 = plt.figure(figsize=(15,7))
ax6 = fig6.add_subplot(111) ax6 = fig6.add_subplot(111)
for i in range(len(phi)): 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.set(xlabel=r"$t \, [yr]$", ylabel=r"$\phi \, [^{\circ}]$")
ax6.legend() ax6.legend()
fig6.suptitle("Inclination angle of the perturbator's orbital plane "+title2) fig6.suptitle("Inclination angle of the perturbator's orbital plane "+title2)

2
lib/units.py
View File

@@ -7,4 +7,4 @@ Units used in the project.
globals()['G'] = 6.67e-11 #Gravitational constant in SI units globals()['G'] = 6.67e-11 #Gravitational constant in SI units
globals()['Ms'] = 2e30 #Solar mass in kg globals()['Ms'] = 2e30 #Solar mass in kg
globals()['au'] = 1.5e11 #Astronomical unit in m globals()['au'] = 1.5e11 #Astronomical unit in m
globals()['yr'] = 3.15576e7 #year in seconds globals()['yr'] = 3.15576e7 #year in seconds

44
main.py
View File

@@ -12,29 +12,40 @@ from lib.units import *
def main(): def main():
#initialisation #initialisation
m = np.array([1., 1., 0.1],dtype=np.longdouble)*Ms#/Ms # Masses in Solar mass m = np.array([1., 1., 0.1],
a = np.array([1.00, 1.00, 10.0],dtype=np.longdouble)*au#/au # Semi-major axis in astronomical units dtype=np.longdouble)*Ms # Masses in Solar mass
e = np.array([0., 0., 0.25],dtype=np.longdouble) # Eccentricity a = np.array([1.00, 1.00, 7.0],
psi = np.array([0., 0., 80.],dtype=np.longdouble)*np.pi/180. # Inclination of the orbital plane in degrees 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]) x1 = a[0]*(1.+e[0])*np.array([0., -1., 0.],
x2 = np.array([0., 1., 0.],dtype=np.longdouble)*a[1]*(1.+e[1]) dtype=np.longdouble)
x3 = np.array([np.cos(psi[2]), 0., np.sin(psi[2])],dtype=np.longdouble)*a[2]*(1.+e[2]) 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) 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) 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.],
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) 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) 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) v = np.array([v1, v2, v3],dtype=np.longdouble)
#integration parameters #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" integrator = "leapfrog"
n_bodies = 3 n_bodies = 3
display = False display = False
gif = False gif = False
blackstyle = True 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 display_param = True
#simulation start #simulation start
@@ -47,15 +58,18 @@ def main():
t1 = time() t1 = time()
if integrator.lower() in ['leapfrog', 'frogleap', 'frog']: 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']: 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() t2=time()
print("...Integration end.\n Elapsed time {0:.3f} sec".format(t2-t1)) print("...Integration end.\n Elapsed time {0:.3f} sec".format(t2-t1))
parameters = [duration, [step], dyn_syst, integrator, [a, e, psi]] 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 return 0
if __name__ == '__main__': if __name__ == '__main__':

BIN
plots/2bodies_hermite_E.png
View File

Binary file not shown.
Side by Side Swipe Overlay

Before

Width:  |  Height:  |  Size: 52 KiB

After

Width:  |  Height:  |  Size: 54 KiB

BIN
plots/2bodies_hermite_L.png
View File

Binary file not shown.
Side by Side Swipe Overlay

Before

Width:  |  Height:  |  Size: 51 KiB

After

Width:  |  Height:  |  Size: 63 KiB

BIN
plots/2bodies_hermite_a_e.png
View File

Binary file not shown.
Side by Side Swipe Overlay

Before

Width:  |  Height:  |  Size: 50 KiB

After

Width:  |  Height:  |  Size: 39 KiB

BIN
plots/2bodies_hermite_dEm.png
View File

Binary file not shown.
Side by Side Swipe Overlay

Before

Width:  |  Height:  |  Size: 43 KiB

After

Width:  |  Height:  |  Size: 46 KiB

BIN
plots/2bodies_hermite_dL2.png
View File

Binary file not shown.
Side by Side Swipe Overlay

Before

Width:  |  Height:  |  Size: 44 KiB

After

Width:  |  Height:  |  Size: 56 KiB

BIN
plots/2bodies_leapfrog_phi.png
View File

Binary file not shown.
Side by Side

Before

Width:  |  Height:  |  Size: 36 KiB