From e8d3e1c0beded5dc1f02938f2a3d241396ceed7b Mon Sep 17 00:00:00 2001 From: Thibault Barnouin Date: Fri, 22 Oct 2021 16:57:07 +0200 Subject: [PATCH] recover_param boolean allow to recover Mechanical Energy and Angular Momentum at each frogleap step --- lib/integrator.py | 22 +++++++++++++--------- 1 file changed, 13 insertions(+), 9 deletions(-) diff --git a/lib/integrator.py b/lib/integrator.py index eebfd00..097f592 100755 --- a/lib/integrator.py +++ b/lib/integrator.py @@ -33,20 +33,29 @@ def frogleap(duration, step, dyn_syst, recover_param=False, display=False): m_array = dyn_syst.get_masses() q_array = dyn_syst.get_positions() p_array = dyn_syst.get_momenta() - + E = np.zeros(N) + L = np.zeros(N) if display: d = DynamicUpdate() d.min_x, d.max_x = -1.5*np.abs(q_array).max(), +1.5*np.abs(q_array).max() d.on_launch() - for _ in range(N): + for j in range(N): # half-step drift q_array, p_array = q_array + step/2*p_array/m_array , p_array # full-step kick q_array, p_array = q_array , p_array - step*dp_dt(m_array, q_array) # half-step drift 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: # 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]) 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: - return dyn_syst + return E, L