Tibeuleu/KozaiLidov
1
0
Fork 0
Code Activity

recover_param boolean allow to recover Mechanical Energy and Angular Momentum at each frogleap step

Browse Source
This commit is contained in:
Thibault Barnouin
2021-10-22 16:57:07 +02:00
parent ec0cf2fd19
commit e8d3e1c0be
1 changed files with 13 additions and 9 deletions
Download Patch File Download Diff File Expand all files Collapse all files

22
lib/integrator.py
View File

@@ -33,20 +33,29 @@ def frogleap(duration, step, dyn_syst, recover_param=False, display=False):
m_array = dyn_syst.get_masses() m_array = dyn_syst.get_masses()
q_array = dyn_syst.get_positions() q_array = dyn_syst.get_positions()
p_array = dyn_syst.get_momenta() p_array = dyn_syst.get_momenta()
E = np.zeros(N)
L = np.zeros(N)
if display: if display:
d = DynamicUpdate() d = DynamicUpdate()
d.min_x, d.max_x = -1.5*np.abs(q_array).max(), +1.5*np.abs(q_array).max() d.min_x, d.max_x = -1.5*np.abs(q_array).max(), +1.5*np.abs(q_array).max()
d.on_launch() d.on_launch()
for _ in range(N): for j in range(N):
# half-step drift # half-step drift
q_array, p_array = q_array + step/2*p_array/m_array , p_array q_array, p_array = q_array + step/2*p_array/m_array , p_array
# full-step kick # full-step kick
q_array, p_array = q_array , p_array - step*dp_dt(m_array, q_array) q_array, p_array = q_array , p_array - step*dp_dt(m_array, q_array)
# half-step drift # half-step drift
q_array, p_array = q_array + step/2*p_array/m_array , p_array q_array, p_array = q_array + step/2*p_array/m_array , p_array
#print(p_array)
for i, body in enumerate(dyn_syst.bodylist):
body.q = q_array[i]
body.p = p_array[i]
if body.m != 0.:
body.v = body.p/body.m
E[j] = dyn_syst.Eval()
L[j] = dyn_syst.Lval()
if display: if display:
# In center of mass frame # In center of mass frame
@@ -55,10 +64,5 @@ def frogleap(duration, step, dyn_syst, recover_param=False, display=False):
d.on_running(q_array[:,0]-q_cm[0], q_array[:,1]-q_cm[1]) d.on_running(q_array[:,0]-q_cm[0], q_array[:,1]-q_cm[1])
time.sleep(0.01) time.sleep(0.01)
for i, body in enumerate(dyn_syst.bodylist):
body.q = q_array[i]
body.p = p_array[i]
if body.m != 0.:
body.v = body.p/body.m
if recover_param: if recover_param:
return dyn_syst return E, L