用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")