Tibeuleu/KozaiLidov
1
0
Fork 0
Code Activity

debug integrator

Browse Source
This commit is contained in:
Thibault Barnouin
2021-10-22 14:10:12 +02:00
parent b096021fc4
commit 71f00c26a8
3 changed files with 17 additions and 9 deletions
Download Patch File Download Diff File Expand all files Collapse all files

16
lib/integrator.py
View File

@@ -16,9 +16,11 @@ def dp_dt(m_array, q_array):
""" """
dp_array = np.zeros(q_array.shape) dp_array = np.zeros(q_array.shape)
for i in range(q_array.shape[0]): for i in range(q_array.shape[0]):
m_j = np.delete(m_array, i)
q_j = np.delete(q_array, i, 0) q_j = np.delete(q_array, i, 0)
dp_array = -m_array[i]*np.sum((m_j/np.sum((q_j-q_array[i])**3, axis=1)).reshape((q_j.shape[0],1))*(q_j-q_array[i]), axis=0) m_j = np.delete(m_array, i).reshape((q_j.shape[0],1))
dp_array[i] = -m_array[i]*np.sum(m_j/np.sum(np.sqrt(np.sum((q_j-q_array[i])**2, axis=0)))**3*(q_j-q_array[i]), axis=0)
dp_array[np.isnan(dp_array)] = 0.
print(dp_array)
return dp_array return dp_array
def frogleap(duration, step, m_array, q_array, p_array, display=False): def frogleap(duration, step, m_array, q_array, p_array, display=False):
@@ -38,11 +40,15 @@ def frogleap(duration, step, m_array, q_array, p_array, display=False):
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)
# In center of mass frame
q_cm = np.sum(m_array.reshape((q_array.shape[0],1))*q_array, axis=0)/m_array.sum()
q_array -= q_cm
if display: if display:
# display progression # display progression
q_cm = np.sum(m_array.reshape((q_array.shape[0],1))*q_array, axis=0)/m_array.sum() d.on_running(q_array[:,0], q_array[:,1])
d.on_running(q_array[:,0]-q_cm[0], q_array[:,1]-q_cm[1]) time.sleep(0.01)
time.sleep(0.1)
return q_array, p_array return q_array, p_array

2
lib/plots.py
View File

@@ -11,6 +11,8 @@ class DynamicUpdate():
min_x = -10 min_x = -10
max_x = 10 max_x = 10
plt.ion()
def on_launch(self): def on_launch(self):
#Set up plot #Set up plot
self.figure, self.ax = plt.subplots() self.figure, self.ax = plt.subplots()

8
main.py
View File

@@ -6,7 +6,7 @@ from lib.integrator import frogleap
def main(): def main():
#initialisation #initialisation
m = np.array([100, 2, 2]) m = np.array([1e5, 1, 1])
x1 = np.array([0, 0, 0]) x1 = np.array([0, 0, 0])
x2 = np.array([1, 0, 0]) x2 = np.array([1, 0, 0])
@@ -14,11 +14,11 @@ def main():
q = np.array([x1, x2, x3]) q = np.array([x1, x2, x3])
v1 = np.array([0, 0, 0]) v1 = np.array([0, 0, 0])
v2 = np.array([1, 0, 0]) v2 = np.array([0, 0, 0])
v3 = np.array([1, 0, 0]) v3 = np.array([0, 0, 0])
p = m*np.array([v1, v2, v3]) p = m*np.array([v1, v2, v3])
q, p = frogleap(10, 1, m, q, p, display=True) q, p = frogleap(10, 0.01, m, q, p, display=True)
return 0 return 0
if __name__ == '__main__': if __name__ == '__main__':