用Python搓一個太陽系

import numpy as np

import matplotlib.pyplot as plt

from matplotlib import animation

m = [1.33e20,3.98e14,4.9e12]

x = np.array([0,1.5e11,1.5e11+3.8e8])

y = np.array([0.0,0,0])

u = np.array([0.0,0,0])

v = np.array([0,2.88e4,2.88e4+1.02e3])

fig = plt.figure(figsize=(12,12))

ax = fig.add_subplot(xlim=(-2e11,2e11),ylim=(-2e11,2e11))

ax.grid()

trace0, = ax.plot([],[],'-', lw=0.5)

trace1, = ax.plot([],[],'-', lw=0.5)

trace2, = ax.plot([],[],'-', lw=0.5)

pt0, = ax.plot([x[0]],[y[0]] ,marker='o')

pt1, = ax.plot([x[0]],[y[0]] ,marker='o')

pt2, = ax.plot([x[0]],[y[0]] ,marker='o')

k_text = ax.text(0.05,0.85,'',transform=ax.transAxes)

textTemplate = 't = %.3f days\n'

N = 1000

dt = 36000

ts =  np.arange(0,N*dt,dt)/3600/24

xs,ys = [],[]

for _ in ts:

     x_ij = (x-x.reshape(3,1))

     y_ij = (y-y.reshape(3,1))

     r_ij = np.sqrt(x_ij**2+y_ij**2)

     for i in range(3):

         for j in range(3):

             if i!=j :

                 u[i] += (m[j]*x_ij[i,j]*dt/r_ij[i,j]**3)

                 v[i] += (m[j]*y_ij[i,j]*dt/r_ij[i,j]**3)

     x += u*dt

     y += v*dt

     xs.append(x.tolist())

     ys.append(y.tolist())

xs = np.array(xs)

ys = np.array(ys)

def animate(n):

     trace0.set_data(xs[:n,0],ys[:n,0])

     trace1.set_data(xs[:n,1],ys[:n,1])

     # 繪圖時的地月距離擴大 100 倍，否則看不清

     tempX2S =外匯跟單gendan5.com xs[:n,1]+100*(xs[:n,2]-xs[:n,1])

     tempY2S = ys[:n,1]+100*(ys[:n,2]-ys[:n,1])

     trace2.set_data(tempX2S,tempY2S)

     pt0.set_data([xs[n,0]],[ys[n,0]])

     pt1.set_data([xs[n,1]],[ys[n,1]])

     tempX = xs[n,1]+100*(xs[n,2]-xs[n,1])

     tempY = ys[n,1]+100*(ys[n,2]-ys[n,1])

     pt2.set_data([tempX],[tempY])

     k_text.set_text(textTemplate % ts[n])

     return trace0, trace1, trace2, pt0, pt1, pt2, k_text

ani = animation.FuncAnimation(fig, animate,

     range(N), interval=10, blit=True)

plt.show()

ani.save("3.gif")