Tibeuleu/KozaiLidov
1
0
Fork 0
Code Activity

Add E and L outputsfor hermite integrator

Browse Source
This commit is contained in:
Thibault Barnouin
2021-11-05 22:40:33 +01:00
parent 04a4dd6e6a
commit b2234afeca
3 changed files with 11 additions and 7 deletions
Download Patch File Download Diff File Expand all files Collapse all files

3
lib/integrator.py
View File

@@ -6,7 +6,6 @@ Implementation of the various integrators for numerical integration.
Comes from the assumption that the problem is analytically defined in position-momentum (q-p) space for a given hamiltonian H. Comes from the assumption that the problem is analytically defined in position-momentum (q-p) space for a given hamiltonian H.
""" """
from os import system from os import system
import time
import numpy as np import numpy as np
from lib.plots import DynamicUpdate from lib.plots import DynamicUpdate
@@ -74,7 +73,7 @@ def frogleap(duration, step, dyn_syst, recover_param=False, display=False):
d.on_running(q_array[0], q_array[1], q_array[2], step=j, label="step {0:d}/{1:d}".format(j,N)) d.on_running(q_array[0], q_array[1], q_array[2], step=j, label="step {0:d}/{1:d}".format(j,N))
else: else:
d.on_running(q_array[:,0], q_array[:,1], q_array[:,2], step=j, label="step {0:d}/{1:d}".format(j,N)) d.on_running(q_array[:,0], q_array[:,1], q_array[:,2], step=j, label="step {0:d}/{1:d}".format(j,N))
time.sleep(1e-5)
if display: if display:
system("convert -delay 5 -loop 0 tmp/?????.png tmp/temp.gif && rm tmp/?????.png") system("convert -delay 5 -loop 0 tmp/?????.png tmp/temp.gif && rm tmp/?????.png")
system("convert tmp/temp.gif -fuzz 30% -layers Optimize plots/dynsyst.gif && rm tmp/temp.gif") system("convert tmp/temp.gif -fuzz 30% -layers Optimize plots/dynsyst.gif && rm tmp/temp.gif")

13
lib/objects.py
View File

@@ -163,7 +163,7 @@ class System:
for body in self.bodylist: for body in self.bodylist:
body.p = body.v*body.m body.p = body.v*body.m
def hermite(self, duration, dt, display=False): def hermite(self, duration, dt, recover_param=False, display=False):
if display: if display:
try: try:
system("mkdir tmp") system("mkdir tmp")
@@ -173,9 +173,14 @@ class System:
d.on_launch() d.on_launch()
N = np.ceil(duration/dt).astype(int) N = np.ceil(duration/dt).astype(int)
E = np.zeros(N)
L = np.zeros((N,3))
for j in range(N): for j in range(N):
self.HPC(dt) self.HPC(dt)
E[j] = self.Eval()
L[j] = self.Lval()
if display: if display:
# display progression # display progression
q_array = self.get_positions() q_array = self.get_positions()
@@ -183,7 +188,7 @@ class System:
d.on_running(q_array[0], q_array[1], q_array[2], step=j, label="step {0:d}/{1:d}".format(j,N)) d.on_running(q_array[0], q_array[1], q_array[2], step=j, label="step {0:d}/{1:d}".format(j,N))
else: else:
d.on_running(q_array[:,0], q_array[:,1], q_array[:,2], step=j, label="step {0:d}/{1:d}".format(j,N)) d.on_running(q_array[:,0], q_array[:,1], q_array[:,2], step=j, label="step {0:d}/{1:d}".format(j,N))
#time.sleep(1e-5)
if recover_param:
return 1 return E, L

2
main.py
View File

@@ -35,7 +35,7 @@ def main():
duration, step = 100*3e7, 1e4 duration, step = 100*3e7, 1e4
E, L = frogleap(duration, step, dyn_syst, recover_param=True)#, display=True) E, L = frogleap(duration, step, dyn_syst, recover_param=True)#, display=True)
#dyn_syst.hermite(duration,step, display=True) E, L = dyn_syst.hermite(duration,step, recover_param=True, display=True)
fig1 = plt.figure(figsize=(30,15)) fig1 = plt.figure(figsize=(30,15))
ax1 = fig1.add_subplot(111) ax1 = fig1.add_subplot(111)