file management
This commit is contained in:
Alex_Hubert
70
lib/LeapFrog.py
Normal file
70
lib/LeapFrog.py
Normal file
@@ -0,0 +1,70 @@
|
||||
#!/usr/bin/python
|
||||
# -*- coding:utf-8 -*-
|
||||
"""
|
||||
Hermite integrator
|
||||
"""
|
||||
|
||||
|
||||
from os import system
|
||||
import numpy as np
|
||||
from lib.plots import DynamicUpdate
|
||||
from lib.units import *
|
||||
|
||||
|
||||
def Drift(dyn_syst, dt):
|
||||
for body in dyn_syst.bodylist:
|
||||
body.q = body.q + dt * body.v
|
||||
|
||||
|
||||
def Kick(dyn_syst, dt):
|
||||
for body in dyn_syst.bodylist:
|
||||
body.a = np.zeros(3)
|
||||
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
|
||||
|
||||
|
||||
def LP(dyn_syst, dt):
|
||||
dyn_syst.COMShift()
|
||||
Drift(dyn_syst, dt / 2)
|
||||
Kick(dyn_syst, dt)
|
||||
Drift(dyn_syst, dt / 2)
|
||||
dyn_syst.time = dyn_syst.time + dt
|
||||
for body in dyn_syst.bodylist:
|
||||
body.p = body.v * body.m
|
||||
|
||||
|
||||
def leapfrog(dyn_syst, duration, dt, recover_param=False, display=False, savename=None):
|
||||
if display:
|
||||
try:
|
||||
system("mkdir tmp")
|
||||
except IOError:
|
||||
system("rm tmp/*")
|
||||
d = DynamicUpdate(dyn_syst)
|
||||
d.launch(dyn_syst.blackstyle)
|
||||
|
||||
N = np.ceil(duration / dt).astype(int)
|
||||
E = np.zeros(N)
|
||||
L = np.zeros((N, 3))
|
||||
for j in range(N):
|
||||
LP(dyn_syst,dt)
|
||||
|
||||
E[j] = dyn_syst.E
|
||||
L[j] = dyn_syst.L
|
||||
|
||||
if display and j % 5 == 0:
|
||||
# display progression
|
||||
if len(dyn_syst.bodylist) == 1:
|
||||
d.on_running(dyn_syst, step=j, label="step {0:d}/{1:d}".format(j, N))
|
||||
else:
|
||||
d.on_running(dyn_syst, step=j, label="step {0:d}/{1:d}".format(j, N))
|
||||
if display:
|
||||
d.close()
|
||||
if not savename is None:
|
||||
system("convert -delay 5 -loop 0 tmp/??????.png tmp/temp.gif && rm tmp/??????.png")
|
||||
system("convert tmp/temp.gif -fuzz 10% -layers Optimize plots/{0:s}_dynsyst.gif".format(savename))
|
||||
|
||||
if recover_param:
|
||||
return E, L
|
||||
115
lib/hermite.py
Normal file
115
lib/hermite.py
Normal file
@@ -0,0 +1,115 @@
|
||||
#!/usr/bin/python
|
||||
# -*- coding:utf-8 -*-
|
||||
"""
|
||||
Hermite integrator
|
||||
"""
|
||||
|
||||
|
||||
from os import system
|
||||
import numpy as np
|
||||
from lib.plots import DynamicUpdate
|
||||
from lib.units import *
|
||||
|
||||
|
||||
def Update_a(dyn_syst): # update acceleration of bodies in system
|
||||
for body in dyn_syst.bodylist:
|
||||
body.a = np.zeros(3)
|
||||
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)
|
||||
|
||||
|
||||
def Update_j(dyn_syst): # update jerk of bodies in system
|
||||
for body in dyn_syst.bodylist:
|
||||
body.j = np.zeros(3)
|
||||
for otherbody in dyn_syst.bodylist:
|
||||
if body != otherbody:
|
||||
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
|
||||
|
||||
|
||||
def Predict(dyn_syst, dt): # update predicted position and velocities 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.)
|
||||
|
||||
|
||||
def Update_ap(dyn_syst): # update acceleration of bodies in system
|
||||
for body in dyn_syst.bodylist:
|
||||
body.ap = np.zeros(3)
|
||||
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)
|
||||
|
||||
|
||||
def Update_jp(dyn_syst): # update jerk of bodies in system
|
||||
for body in dyn_syst.bodylist:
|
||||
body.jp = np.zeros(3)
|
||||
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
|
||||
|
||||
|
||||
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)
|
||||
|
||||
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
|
||||
COMShift(dyn_syst)
|
||||
Update_a(dyn_syst)
|
||||
Update_j(dyn_syst)
|
||||
Predict(dyn_syst, dt)
|
||||
Update_ap(dyn_syst)
|
||||
Update_jp(dyn_syst)
|
||||
Correct(dyn_syst, dt)
|
||||
dyn_syst.time = dyn_syst.time + dt
|
||||
for body in dyn_syst.bodylist:
|
||||
body.p = body.v * body.m
|
||||
|
||||
|
||||
def hermite(dyn_syst, duration, dt, recover_param=False, display=False, savename=None):
|
||||
if display:
|
||||
try:
|
||||
system("mkdir tmp")
|
||||
except IOError:
|
||||
system("rm tmp/*")
|
||||
d = DynamicUpdate(dyn_syst)
|
||||
d.launch(dyn_syst.blackstyle)
|
||||
|
||||
N = np.ceil(duration / dt).astype(int)
|
||||
E = np.zeros(N)
|
||||
L = np.zeros((N, 3))
|
||||
for j in range(N):
|
||||
HPC(dyn_syst, dt)
|
||||
|
||||
E[j] = dyn_syst.E
|
||||
L[j] = dyn_syst.L
|
||||
|
||||
if display and j % 100 == 0:
|
||||
# display progression
|
||||
if len(dyn_syst.bodylist) == 1:
|
||||
d.on_running(dyn_syst, step=j, label="step {0:d}/{1:d}".format(j, N))
|
||||
else:
|
||||
d.on_running(dyn_syst, step=j, label="step {0:d}/{1:d}".format(j, N))
|
||||
if display:
|
||||
d.close()
|
||||
if not savename is None:
|
||||
system("convert -delay 5 -loop 0 tmp/??????.png tmp/temp.gif && rm tmp/??????.png")
|
||||
system("convert tmp/temp.gif -fuzz 10% -layers Optimize plots/{0:s}_dynsyst.gif".format(savename))
|
||||
|
||||
if recover_param:
|
||||
return E, L
|
||||
147
lib/objects.py
147
lib/objects.py
@@ -8,185 +8,138 @@ import numpy as np
|
||||
from lib.plots import DynamicUpdate
|
||||
from lib.units import *
|
||||
|
||||
|
||||
class Body:
|
||||
|
||||
def __init__(self, mass, position, velocity):
|
||||
self.m = mass
|
||||
self.q = position
|
||||
self.v = velocity
|
||||
self.p = velocity*mass
|
||||
self.p = velocity * mass
|
||||
self.a = np.zeros(3)
|
||||
self.ap = np.zeros(3)
|
||||
self.j = np.zeros(3)
|
||||
self.jp = np.zeros(3)
|
||||
self.qp = np.zeros(3)
|
||||
self.vp = np.zeros(3)
|
||||
|
||||
def __repr__(self): # Called upon "print(body)"
|
||||
return r"Body of mass: {0:.2f} $M_\odot$, position: {1}, velocity: {2}".format(self.m/Ms, self.q, self.v)
|
||||
def __repr__(self): # Called upon "print(body)"
|
||||
return r"Body of mass: {0:.2f} $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:.2f} $M_\odot$".format(self.m / Ms)
|
||||
|
||||
def __str__(self): # Called upon "str(body)"
|
||||
return r"Body of mass: {0:.2f} $M_\odot$".format(self.m/Ms)
|
||||
|
||||
class System(Body):
|
||||
|
||||
def __init__(self, bodylist, blackstyle=True):
|
||||
self.blackstyle = blackstyle #for dark mode in plot
|
||||
self.blackstyle = blackstyle # for dark mode in plot
|
||||
self.bodylist = np.array(bodylist)
|
||||
self.time = 0 #lifetime of system
|
||||
self.time = 0 # lifetime of system
|
||||
self.m = self.M
|
||||
self.q = self.COM
|
||||
self.v = self.COMV
|
||||
|
||||
|
||||
def __repr__(self): # Called upon "print(system)"
|
||||
return str([print(body) for body in self.bodylist])
|
||||
|
||||
def __str__(self): # Called upon "str(system)"
|
||||
return str([str(body) for body in self.bodylist])
|
||||
|
||||
@property
|
||||
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])
|
||||
|
||||
|
||||
@property
|
||||
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])
|
||||
ydata = np.array([body.q[1] for body in self.bodylist])
|
||||
zdata = np.array([body.q[2] for body in self.bodylist])
|
||||
return xdata, ydata, zdata
|
||||
|
||||
|
||||
@property
|
||||
def get_velocities(self): #return the positions of the bodies
|
||||
def get_velocities(self): # return the positions of the bodies
|
||||
vxdata = np.array([body.v[0] for body in self.bodylist])
|
||||
vydata = np.array([body.v[1] for body in self.bodylist])
|
||||
vzdata = np.array([body.v[2] for body in self.bodylist])
|
||||
return vxdata, vydata, vzdata
|
||||
|
||||
|
||||
@property
|
||||
def get_momenta(self): #return the momenta of the bodies
|
||||
def get_momenta(self): # return the momenta of the bodies
|
||||
pxdata = np.array([body.p[0] for body in self.bodylist])
|
||||
pydata = np.array([body.p[1] for body in self.bodylist])
|
||||
pzdata = np.array([body.p[2] for body in self.bodylist])
|
||||
return pxdata, pydata, pzdata
|
||||
|
||||
@property
|
||||
def M(self): #return total system mass
|
||||
def M(self): # return total system mass
|
||||
mass = 0
|
||||
for body in self.bodylist:
|
||||
mass = mass + body.m
|
||||
return mass
|
||||
|
||||
@property
|
||||
def COM(self): #return center of mass in cartesian np_array
|
||||
def COM(self): # return center of mass in cartesian np_array
|
||||
coord = np.zeros(3)
|
||||
for body in self.bodylist:
|
||||
coord = coord + body.m*body.q
|
||||
coord = coord/self.M
|
||||
coord = coord + body.m * body.q
|
||||
coord = coord / self.M
|
||||
return coord
|
||||
|
||||
@property
|
||||
def COMV(self): #return center of mass velocity in cartesian np_array
|
||||
def COMV(self): # return center of mass velocity in cartesian np_array
|
||||
coord = np.zeros(3)
|
||||
for body in self.bodylist:
|
||||
coord = coord + body.p
|
||||
coord = coord/self.M
|
||||
coord = coord / self.M
|
||||
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 frame and set COM at rest
|
||||
for body in self.bodylist:
|
||||
body.q = body.q - self.COM
|
||||
body.p = body.p - self.COMV
|
||||
|
||||
@property
|
||||
def L(self): #return angular momentum of bodies in system
|
||||
def L(self): # return angular momentum of bodies in system
|
||||
L = np.zeros(3)
|
||||
for body in self.bodylist:
|
||||
L = L + np.cross(body.q,body.p)
|
||||
L = L + np.cross(body.q, body.p)
|
||||
return L
|
||||
|
||||
@property
|
||||
def E(self): #return total energy of bodies in system
|
||||
def E(self): # return total energy of bodies in system
|
||||
T = 0
|
||||
W = 0
|
||||
for body in self.bodylist:
|
||||
T = T + 1./2.*body.m*np.linalg.norm(body.v)**2
|
||||
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*body.m*otherbody.m/rij
|
||||
rij = np.linalg.norm(body.q - otherbody.q)
|
||||
W = W - G * body.m * otherbody.m / rij
|
||||
E = T + W
|
||||
return E
|
||||
|
||||
|
||||
@property
|
||||
def mu(self):
|
||||
sum = 0
|
||||
prod = 1
|
||||
for body in self.bodylist:
|
||||
prod = prod * body.m
|
||||
mu = prod/self.M
|
||||
mu = prod / self.M
|
||||
return mu
|
||||
|
||||
@property
|
||||
def ex(self): #exentricity of system (if composed of 2 bodies)
|
||||
if len(self.bodylist) != 2 :
|
||||
def ex(self): # exentricity of system (if composed of 2 bodies)
|
||||
if len(self.bodylist) != 2:
|
||||
return np.nan
|
||||
else:
|
||||
k = (2.*self.E*(np.linalg.norm(self.L)**2))/((G**2)*(self.M**2)*(self.mu**3)) + 1.
|
||||
k = (2. * self.E * (np.linalg.norm(self.L) ** 2)) / ((G ** 2) * (self.M ** 2) * (self.mu ** 3)) + 1.
|
||||
return k
|
||||
|
||||
@property
|
||||
def sma(self): #semi major axis of system (if composed of 2 bodies)
|
||||
if len(self.bodylist) != 2 :
|
||||
def sma(self): # semi major axis of system (if composed of 2 bodies)
|
||||
if len(self.bodylist) != 2:
|
||||
return np.nan
|
||||
else:
|
||||
sma = -G*self.M*self.mu/(2.*self.E)
|
||||
sma = -G * self.M * self.mu / (2. * self.E)
|
||||
return sma
|
||||
|
||||
def Drift(self, dt):
|
||||
for body in self.bodylist:
|
||||
body.q = body.q + dt*body.v
|
||||
|
||||
def Kick(self, dt):
|
||||
for body in self.bodylist:
|
||||
body.a = np.zeros(3)
|
||||
for otherbody in self.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
|
||||
|
||||
def LP(self, dt):
|
||||
self.COMShift()
|
||||
self.Drift(dt/2)
|
||||
self.Kick(dt)
|
||||
self.Drift(dt/2)
|
||||
self.time = self.time + dt
|
||||
for body in self.bodylist:
|
||||
body.p = body.v*body.m
|
||||
|
||||
def leapfrog(self, duration, dt, recover_param=False, display=False, savename=None):
|
||||
if display:
|
||||
try:
|
||||
system("mkdir tmp")
|
||||
except IOError:
|
||||
system("rm tmp/*")
|
||||
d = DynamicUpdate(self)
|
||||
d.launch(self.blackstyle)
|
||||
|
||||
N = np.ceil(duration/dt).astype(int)
|
||||
E = np.zeros(N)
|
||||
L = np.zeros((N,3))
|
||||
for j in range(N):
|
||||
self.LP(dt)
|
||||
|
||||
E[j] = self.E
|
||||
L[j] = self.L
|
||||
|
||||
if display and j%5==0:
|
||||
# display progression
|
||||
if len(self.bodylist) == 1:
|
||||
d.on_running(self, step=j, label="step {0:d}/{1:d}".format(j,N))
|
||||
else:
|
||||
d.on_running(self, step=j, label="step {0:d}/{1:d}".format(j,N))
|
||||
if display:
|
||||
d.close()
|
||||
if not savename is None:
|
||||
system("convert -delay 5 -loop 0 tmp/??????.png tmp/temp.gif && rm tmp/??????.png")
|
||||
system("convert tmp/temp.gif -fuzz 10% -layers Optimize plots/{0:s}_dynsyst.gif".format(savename))
|
||||
|
||||
if recover_param:
|
||||
return E, L
|
||||
|
||||
def __repr__(self): # Called upon "print(system)"
|
||||
return str([print(body) for body in self.bodylist])
|
||||
|
||||
def __str__(self): # Called upon "str(system)"
|
||||
return str([str(body) for body in self.bodylist])
|
||||
|
||||
@@ -1,238 +0,0 @@
|
||||
#!/usr/bin/python
|
||||
# -*- coding:utf-8 -*-
|
||||
"""
|
||||
Class definition for physical attribute
|
||||
"""
|
||||
from os import system
|
||||
import numpy as np
|
||||
from lib.plots import DynamicUpdate
|
||||
from lib.units import *
|
||||
|
||||
class Body:
|
||||
|
||||
def __init__(self, mass, position, velocity):
|
||||
self.m = mass
|
||||
self.q = position
|
||||
self.v = velocity
|
||||
self.p = velocity*mass
|
||||
self.a = np.zeros(3)
|
||||
self.ap = np.zeros(3)
|
||||
self.j = np.zeros(3)
|
||||
self.jp = np.zeros(3)
|
||||
self.qp = np.zeros(3)
|
||||
self.vp = np.zeros(3)
|
||||
|
||||
def __repr__(self): # Called upon "print(body)"
|
||||
return r"Body of mass: {0:.2f} $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:.2f} $M_\odot$".format(self.m/Ms)
|
||||
|
||||
class System(Body):
|
||||
|
||||
def __init__(self, bodylist, blackstyle=True):
|
||||
self.blackstyle = blackstyle #for dark mode in plot
|
||||
self.bodylist = np.array(bodylist)
|
||||
self.time = 0 #lifetime of system
|
||||
self.m = self.M
|
||||
self.q = self.COM
|
||||
self.v = self.COMV
|
||||
|
||||
@property
|
||||
def get_masses(self): #return the masses of each object
|
||||
return np.array([body.m for body in self.bodylist])
|
||||
|
||||
@property
|
||||
def get_positions(self): #return the positions of the bodies
|
||||
xdata = np.array([body.q[0] for body in self.bodylist])
|
||||
ydata = np.array([body.q[1] for body in self.bodylist])
|
||||
zdata = np.array([body.q[2] for body in self.bodylist])
|
||||
return xdata, ydata, zdata
|
||||
|
||||
@property
|
||||
def get_velocities(self): #return the positions of the bodies
|
||||
vxdata = np.array([body.v[0] for body in self.bodylist])
|
||||
vydata = np.array([body.v[1] for body in self.bodylist])
|
||||
vzdata = np.array([body.v[2] for body in self.bodylist])
|
||||
return vxdata, vydata, vzdata
|
||||
|
||||
@property
|
||||
def get_momenta(self): #return the momenta of the bodies
|
||||
pxdata = np.array([body.p[0] for body in self.bodylist])
|
||||
pydata = np.array([body.p[1] for body in self.bodylist])
|
||||
pzdata = np.array([body.p[2] for body in self.bodylist])
|
||||
return pxdata, pydata, pzdata
|
||||
|
||||
@property
|
||||
def M(self): #return total system mass
|
||||
mass = 0
|
||||
for body in self.bodylist:
|
||||
mass = mass + body.m
|
||||
return mass
|
||||
|
||||
@property
|
||||
def COM(self): #return center of mass in cartesian np_array
|
||||
coord = np.zeros(3)
|
||||
for body in self.bodylist:
|
||||
coord = coord + body.m*body.q
|
||||
coord = coord/self.M
|
||||
return coord
|
||||
|
||||
@property
|
||||
def COMV(self): #return center of mass velocity in cartesian np_array
|
||||
coord = np.zeros(3)
|
||||
for body in self.bodylist:
|
||||
coord = coord + body.p
|
||||
coord = coord/self.M
|
||||
return coord
|
||||
|
||||
def COMShift(self): #Shift coordinates of bodies in system to COM frame and set COM at rest
|
||||
for body in self.bodylist:
|
||||
body.q = body.q - self.COM
|
||||
body.p = body.p - self.COMV
|
||||
|
||||
@property
|
||||
def L(self): #return angular momentum of bodies in system
|
||||
L = np.zeros(3)
|
||||
for body in self.bodylist:
|
||||
L = L + np.cross(body.q,body.p)
|
||||
return L
|
||||
|
||||
@property
|
||||
def E(self): #return total energy of bodies in system
|
||||
T = 0
|
||||
W = 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*body.m*otherbody.m/rij
|
||||
E = T + W
|
||||
|
||||
def Update_a(self): #update acceleration of bodies in system
|
||||
for body in self.bodylist:
|
||||
body.a = np.zeros(3)
|
||||
for otherbody in self.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)
|
||||
|
||||
def Update_j(self): #update jerk of bodies in system
|
||||
for body in self.bodylist:
|
||||
body.j = np.zeros(3)
|
||||
for otherbody in self.bodylist:
|
||||
if body != otherbody:
|
||||
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
|
||||
|
||||
def Predict(self,dt): # update predicted position and velocities of bodies in system
|
||||
for body in self.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.)
|
||||
|
||||
def Update_ap(self): #update acceleration of bodies in system
|
||||
for body in self.bodylist:
|
||||
body.ap = np.zeros(3)
|
||||
for otherbody in self.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)
|
||||
|
||||
def Update_jp(self): #update jerk of bodies in system
|
||||
for body in self.bodylist:
|
||||
body.jp = np.zeros(3)
|
||||
for otherbody in self.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
|
||||
|
||||
def Correct(self,dt): # correct position and velocities of bodies in system
|
||||
for body in self.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)
|
||||
|
||||
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(self, dt): # update position and velocities of bodies in system with hermite predictor corrector
|
||||
self.COMShift()
|
||||
self.Update_a()
|
||||
self.Update_j()
|
||||
self.Predict(dt)
|
||||
self.Update_ap()
|
||||
self.Update_jp()
|
||||
self.Correct(dt)
|
||||
self.time = self.time + dt
|
||||
for body in self.bodylist:
|
||||
body.p = body.v*body.m
|
||||
|
||||
def hermite(self, duration, dt, recover_param=False, display=False, savename=None):
|
||||
if display:
|
||||
try:
|
||||
system("mkdir tmp")
|
||||
except IOError:
|
||||
system("rm tmp/*")
|
||||
d = DynamicUpdate(self)
|
||||
d.launch(self.blackstyle)
|
||||
|
||||
N = np.ceil(duration/dt).astype(int)
|
||||
E = np.zeros(N)
|
||||
L = np.zeros((N,3))
|
||||
for j in range(N):
|
||||
self.HPC(dt)
|
||||
|
||||
E[j] = self.E
|
||||
L[j] = self.L
|
||||
|
||||
if display and j%100==0:
|
||||
# display progression
|
||||
if len(self.bodylist) == 1:
|
||||
d.on_running(self, step=j, label="step {0:d}/{1:d}".format(j,N))
|
||||
else:
|
||||
d.on_running(self, step=j, label="step {0:d}/{1:d}".format(j,N))
|
||||
if display:
|
||||
d.close()
|
||||
if not savename is None:
|
||||
system("convert -delay 5 -loop 0 tmp/??????.png tmp/temp.gif && rm tmp/??????.png")
|
||||
system("convert tmp/temp.gif -fuzz 10% -layers Optimize plots/{0:s}_dynsyst.gif".format(savename))
|
||||
|
||||
if recover_param:
|
||||
return E, L
|
||||
|
||||
@property
|
||||
def mu(self):
|
||||
sum = 0
|
||||
prod = 1
|
||||
for body in self.bodylist:
|
||||
prod = prod * body.m
|
||||
mu = prod/self.M
|
||||
return mu
|
||||
|
||||
@property
|
||||
def ex(self): #exentricity of system (if composed of 2 bodies)
|
||||
if len(self.bodylist) != 2 :
|
||||
return np.nan
|
||||
else:
|
||||
k = (2.*self.E*(np.linalg.norm(self.L)**2))/((G**2)*(self.M**2)*(self.mu**3)) + 1.
|
||||
return k
|
||||
|
||||
@property
|
||||
def sma(self): #semi major axis of system (if composed of 2 bodies)
|
||||
if len(self.bodylist) != 2 :
|
||||
return np.nan
|
||||
else:
|
||||
sma = -G*self.M*self.mu/(2.*self.E)
|
||||
return sma
|
||||
|
||||
def __repr__(self): # Called upon "print(system)"
|
||||
return str([print(body) for body in self.bodylist])
|
||||
|
||||
def __str__(self): # Called upon "str(system)"
|
||||
return str([str(body) for body in self.bodylist])
|
||||
6
main.py
6
main.py
@@ -4,6 +4,8 @@ from sys import exit as sysexit
|
||||
import numpy as np
|
||||
import matplotlib.pyplot as plt
|
||||
from lib.objects import Body, System
|
||||
from lib.LeapFrog import leapfrog
|
||||
from lib.hermite import hermite
|
||||
from lib.plots import display_parameters
|
||||
from lib.units import *
|
||||
|
||||
@@ -41,9 +43,9 @@ def main():
|
||||
E, L = [], []
|
||||
for step0 in step:
|
||||
if integrator.lower() in ['leapfrog', 'frogleap', 'frog']:
|
||||
E0, L0 = dyn_syst.leapfrog(duration, step0, recover_param=True, display=display, savename=savename)
|
||||
E0, L0 = leapfrog(dyn_syst,duration, step0, recover_param=True, display=display, savename=savename)
|
||||
elif integrator.lower() in ['hermite','herm']:
|
||||
E0, L0 = dyn_syst.hermite(duration, step0, recover_param=True, display=display, savename=savename)
|
||||
E0, L0 = hermite(dyn_syst, duration, step0, recover_param=True, display=display, savename=savename)
|
||||
E.append(E0)
|
||||
L.append(L0)
|
||||
|
||||
|
||||
Binary file not shown.
|
Before Width: | Height: | Size: 118 KiB After Width: | Height: | Size: 120 KiB |
Binary file not shown.
|
Before Width: | Height: | Size: 104 KiB After Width: | Height: | Size: 78 KiB |
Reference in New Issue
Block a user