diff --git a/lib/integrator.py b/lib/integrator.py index 1cd5c6c..d148160 100755 --- a/lib/integrator.py +++ b/lib/integrator.py @@ -16,9 +16,11 @@ def dp_dt(m_array, q_array): """ dp_array = np.zeros(q_array.shape) for i in range(q_array.shape[0]): - m_j = np.delete(m_array, i) 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 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) # half-step drift 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: # 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_cm[0], q_array[:,1]-q_cm[1]) - time.sleep(0.1) + d.on_running(q_array[:,0], q_array[:,1]) + time.sleep(0.01) return q_array, p_array diff --git a/lib/plots.py b/lib/plots.py index ee1e6d9..a1e99d5 100755 --- a/lib/plots.py +++ b/lib/plots.py @@ -11,6 +11,8 @@ class DynamicUpdate(): min_x = -10 max_x = 10 + plt.ion() + def on_launch(self): #Set up plot self.figure, self.ax = plt.subplots() diff --git a/main.py b/main.py index 9d56173..2de7042 100755 --- a/main.py +++ b/main.py @@ -6,7 +6,7 @@ from lib.integrator import frogleap def main(): #initialisation - m = np.array([100, 2, 2]) + m = np.array([1e5, 1, 1]) x1 = np.array([0, 0, 0]) x2 = np.array([1, 0, 0]) @@ -14,11 +14,11 @@ def main(): q = np.array([x1, x2, x3]) v1 = np.array([0, 0, 0]) - v2 = np.array([1, 0, 0]) - v3 = np.array([1, 0, 0]) + v2 = np.array([0, 0, 0]) + v3 = np.array([0, 0, 0]) 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 if __name__ == '__main__':