2

I'm looking to animate my graph (below) and I'm not sure where or how to start since I have no experience animating. I'm not sure how it works or what the structure of the code should be, so if someone can offer a pseudo-code or an algorithm, I would greatly appreciate it. I have provided the code I used to graph the plot below, too.

enter code here
from scipy.integrate import odeint
import scipy as sci
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
import matplotlib.animation as ani

# Universal Gravitational Const.

G = 6.674e-11

# Defining Mass

m1 = 2
m2 = 3.5
m3 = 2.3



# Init positions in graph (array)

pos1 = [-5,0,1]
pos2 = [5,0,10]
pos3 = [0,1,3]

p01 = np.array(pos1)
p02 = np.array(pos2)
p03 = np.array(pos3)


# Init velocities (array)

vi1 = [1,0.01,0]
vi2 = [-5,0,1]
vi3 = [0,-1,0]


v01 = np.array(vi1)
v02 = np.array(vi2)
v03 = np.array(vi3)

#Function
def derivs_func(y,t,G,m1,m2,m3):
    d1 = y[:3]
    d2 = y[3:6]   
    d3 = y[6:9]
    v1 = y[9:12]
    v2 = y[12:15]
    v3 = y[15:18]


    dist12 = np.sqrt((pos2[0]-pos1[0])**2 + (pos2[1]-pos1[1])**2 + (pos2[2]-pos1[2])**2)
    dist13 = np.sqrt((pos3[0]-pos1[0])**2 + (pos3[1]-pos1[1])**2 + (pos3[2]-pos1[2])**2)
    dist23 = np.sqrt((pos3[0]-pos2[0])**2 + (pos3[1]-pos2[1])**2 + (pos3[2]-pos2[2])**2)

    dv1dt = m2 * (d2-d1)/dist12**3 + m3 * (d3-d1)/dist13**3 
    dv2dt = m1 * (d1-d2)/dist12**3 + m3 * (d3-d2)/dist23**3 
    dv3dt = m1 * (d1-d3)/dist13**3 + m2 * (d2-d3)/dist23**3 
    dd1dt = v1 
    dd2dt = v2
    dd3dt = v3

    derivs = np.array([dd1dt,dd2dt,dd3dt,dv1dt,dv2dt,dv3dt]) 
    derivs3 = derivs.flatten() 

    return derivs3 

yo = np.array([p01, p02, p03, v01, v02, v03]) 
y0 = yo.flatten()  

time = np.linspace(0,200,500) 

sol = odeint(derivs_func, y0, time, args = (G,m1,m2,m3)) 

x1 = sol[:,:3]
x2 = sol[:,3:6]
x3 = sol[:,6:9]


fig = plt.figure(figsize = (15,15))
ax = fig.add_subplot(111,projection = '3d')

ax.plot(x1[:,0],x1[:,1],x1[:,2],color = 'b')
ax.plot(x2[:,0],x2[:,1],x2[:,2],color = 'm')
ax.plot(x3[:,0],x3[:,1],x3[:,2],color = 'g')

ax.scatter(x1[-1,0],x1[-1,1],x1[-1,2],color = 'b', marker = 'o', s=30, label = 'Mass 1')
ax.scatter(x2[-1,0],x2[-1,1],x2[-1,2],color = 'm', marker = 'o',s=90, label = 'Mass 2')
ax.scatter(x3[-1,0],x3[-1,1],x3[-1,2],color = 'g', marker = 'o',s=60, label = 'Mass 3')
ax.legend()

enter image description here

theheretic
  • 57
  • 6
  • Hi, I think you've got some of the correct ideas. The problem is that you're treating the 3 dimensions as independent which you can't do. At each step the derivatives need to change and these derivatives depend on the distance between the objects, for which you need all three dimensions and you need to update the distance every time step rather than use the initial distance like you do. Try following this guide (https://towardsdatascience.com/modelling-the-three-body-problem-in-classical-mechanics-using-python-9dc270ad7767) and that should show you how to use all 3 dims at once. – Ianhi Nov 24 '19 at 22:40
  • Additionally your equation for the change in velocity is not quite correct. It should depend on 1/dist^3, and the change in position should not depend on the relative positions, only on the velocity and the time step. – Ianhi Nov 24 '19 at 22:42
  • I took a look at the website that was given, and I followed the procedures correctly but my graph still doesn't show complete orbits. I'm not sure if there is something else wrong with the equations I'm using or the time or the initial conditions. I'm not sure where my error is. I've also updated my code and the graph from the previous question. – theheretic Nov 27 '19 at 02:02

0 Answers0