三体系统

code

from pylab import *
from math import *

me=1.0
mj=316.7*1000
ms=333333.3
dt=0.001
t=40
LIM=t/dt
x=[0 for i in range(0,int(LIM))]
xe=x[:];ye=x[:];vex=x[:];vey=x[:]
xj=x[:];yj=x[:];vjx=x[:];vjy=x[:]
xs=x[:];ys=x[:];vsx=x[:];vsy=x[:]
t=x[:];res=x[:];rjs=x[:];rej=x[:]
def init(a,e):
    x0=a*(1+e)
    vy0=2*math.pi*sqrt((1-e)/(a*(1+e)))
    return [x0,vy0]
xe[0]=init(1,0.017)[0]
vey[0]=init(1,0.017)[1]
xj[0]=init(5.2,0.048)[0]
vjy[0]=init(5.2,0.048)[1]

#consider the motivation of the sun
vsx[0]=0
vsy[0]=-(me*vey[0]+mj*vjy[0])/ms


res[0]=sqrt(pow(xe[0]-xs[0],2)+pow(ye[0]-ys[0],2))
rjs[0]=sqrt(pow(xj[0]-xs[0],2)+pow(yj[0]-ys[0],2))  
rej[0]=sqrt(pow(xe[0]-xj[0],2)+pow(ye[0]-yj[0],2))

for i in range (0,int(LIM-1)):
    vex[i+1]=vex[i]-4*pow(math.pi,2)*(xe[i]-xs[i])/pow(res[i],3)*dt-4*pow(math.pi,2)*(mj/ms)*(xe[i]-xj[i])/pow(rej[i],3)*dt
    vey[i+1]=vey[i]-4*pow(math.pi,2)*(ye[i]-ys[i])/pow(res[i],3)*dt-4*pow(math.pi,2)*(mj/ms)*(ye[i]-yj[i])/pow(rej[i],3)*dt
    vjx[i+1]=vjx[i]-4*pow(math.pi,2)*(xj[i]-xs[i])/pow(rjs[i],3)*dt-4*pow(math.pi,2)*(me/ms)*(xe[i]-xj[i])/pow(rej[i],3)*dt
    vjy[i+1]=vjy[i]-4*pow(math.pi,2)*(yj[i]-ys[i])/pow(rjs[i],3)*dt-4*pow(math.pi,2)*(me/ms)*(ye[i]-yj[i])/pow(rej[i],3)*dt
    vsx[i+1]=vsx[i]+dt*4*pow(math.pi,2)*(me/ms)*(xe[i]-xs[i])/pow(res[i],3)+dt*4*pow(math.pi,2)*(mj/ms)*(xj[i]-xs[i])/pow(rjs[i],3)
    vsy[i+1]=vsy[i]+dt*4*pow(math.pi,2)*(me/ms)*(ye[i]-ys[i])/pow(res[i],3)+dt*4*pow(math.pi,2)*(mj/ms)*(yj[i]-ys[i])/pow(rjs[i],3)    
    xe[i+1]=xe[i]+vex[i+1]*dt
    ye[i+1]=ye[i]+vey[i+1]*dt
    xj[i+1]=xj[i]+vjx[i+1]*dt
    yj[i+1]=yj[i]+vjy[i+1]*dt
    xs[i+1]=xs[i]+vsx[i+1]*dt
    ys[i+1]=ys[i]+vsy[i+1]*dt
    res[i+1]=math.sqrt(pow(xe[i+1]-xs[i+1],2)+pow(ye[i+1]-ys[i+1],2))
    rjs[i+1]=math.sqrt(pow(xj[i+1]-xs[i+1],2)+pow(yj[i+1]-ys[i+1],2))
    rej[i+1]=math.sqrt(pow(xe[i+1]-xj[i+1],2)+pow(ye[i+1]-yj[i+1],2))
    t[i+1]=t[i]+dt

#plot
figure(figsize=[8,8])
plot(xe,ye,color='black')
plot(xj,yj,color='green')
plot(xs,ys,color='red')
legend(('Earth','Jupiter','Sun'),'upper left')
title('three-body trajectory M_j',fontsize=15)
xlabel('x/AU')
ylabel('y/AU')
savefig('three-body trajectory.png')
show()
#3D plot
fig = figure()
ax = fig.add_subplot(111, projection='3d')
ax.plot(xe,ye,t,color='black')
ax.plot(xj,yj,t,color='green')
ax.plot(xs,ys,t,color='red')
legend(('Earth','Jupiter','Sun'),'upper left')
ax.set_xlabel('x/Au')
ax.set_ylabel('y/AU')
ax.set_zlabel('t/yr')
savefig('three-body trajectory 3D.png')
show()