一篇文章教你用Python繪畫一個(gè)太陽(yáng)系

更新時(shí)間：2021年10月14日 12:05:09 作者：微小冷

這篇文章主要給大家介紹了關(guān)于如何利用Python繪畫一個(gè)太陽(yáng)系，文中通過示例代碼介紹的非常詳細(xì)，對(duì)大家的學(xué)習(xí)或者工作具有一定的參考學(xué)習(xí)價(jià)值，需要的朋友們下面隨著小編來一起學(xué)習(xí)學(xué)習(xí)吧

日地月三體

所謂三體，就是三個(gè)物體在重力作用下的運(yùn)動(dòng)。由于三點(diǎn)共面，所以三個(gè)質(zhì)點(diǎn)僅在重力作用下的運(yùn)動(dòng)軌跡也必然無(wú)法逃離平面。

三體運(yùn)動(dòng)所遵循的規(guī)律就是古老而經(jīng)典的萬(wàn)有引力

在這里插入圖片描述

則對(duì)于 m i 而言，

在這里插入圖片描述

且

在這里插入圖片描述

將其寫為差分形式

在這里插入圖片描述

由于我們希望觀察三體運(yùn)動(dòng)的復(fù)雜形式，而不關(guān)系其隨對(duì)應(yīng)的宇宙星體，所以不必考慮單位制，將其在二維平面坐標(biāo)系中拆分，則

在這里插入圖片描述

#后續(xù)代碼主要更改這里的參數(shù)
m = [1.33e20,3.98e14,4.9e12]
x = np.array([0,1.5e11,1.5e11+3.8e8])
y = np.array([0,0,0])
u = np.array([0,0,0])
v = np.array([0,2.88e4,1.02e3])

由于地月之間的距離相對(duì)于日地距離太近，所以在畫圖的時(shí)候?qū)⑵鋽U(kuò)大100倍，得到圖像

在這里插入圖片描述

盡管存在誤差，但最起碼看到了地球圍繞太陽(yáng)轉(zhuǎn)，月球圍繞地球轉(zhuǎn)。。。代碼為

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])
    #繪圖時(shí)的地月距離擴(kuò)大100倍，否則看不清
    tempX2S = 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")

日地火

在這里插入圖片描述

m = [1.33e20,3.98e14,4.28e13]
x = np.array([0,1.5e11,2.28e11])
y = np.array([0.0,0,0])
u = np.array([0.0,0,0])
v = np.array([0,2.88e4,2.4e4])
### 由于火星離地球很遠(yuǎn)，所以不必再改變尺度
def animate(n):
    trace0.set_data(xs[:n,0],ys[:n,0])
    trace1.set_data(xs[:n,1],ys[:n,1])
    trace2.set_data(xs[:n,2],ys[:n,2])
    pt0.set_data([xs[n,0]],[ys[n,0]])
    pt1.set_data([xs[n,1]],[ys[n,1]])
    pt2.set_data([xs[n,2]],[ys[n,2]])
    k_text.set_text(textTemplate % ts[n])
    return trace0, trace1, trace2, pt0, pt1, pt2, k_text

得到

在這里插入圖片描述

這個(gè)運(yùn)動(dòng)要比月球的運(yùn)動(dòng)簡(jiǎn)單得多——前提是開上帝視角，俯瞰太陽(yáng)系。如果站在地球上觀測(cè)火星的運(yùn)動(dòng)，那么這個(gè)運(yùn)動(dòng)可能相當(dāng)帶感

在這里插入圖片描述

所以這都能找到規(guī)律，托勒密那幫人也真夠有才的。

太陽(yáng)系

由于太陽(yáng)和其他星體之間的質(zhì)量相差懸殊，所以太陽(yáng)系內(nèi)的多體運(yùn)動(dòng)，都將退化為二體問題，甚至如果把太陽(yáng)當(dāng)作不動(dòng)點(diǎn)，那就成了單體問題了。

盡管如此，我們還是盡可能地模仿一下太陽(yáng)系的運(yùn)動(dòng)情況

	質(zhì)量	半長(zhǎng)軸(AU)	平均速度(km/s)
水星	0.055	0.387	47.89
金星	0.815	0.723	35.03
地球	1	1	29.79
火星	0.107	1.524	24.13
木星	317.8	5.203	13.06
土星	95.16	9.537	9.64
天王星	14.54	19.19	6.81
海王星	17.14	30.07	5.43
冥王星

除了水星偏心率為0.2，對(duì)黃道面傾斜為7°之外，其余行星的偏心率皆小于0.1，且對(duì)黃道面傾斜普遍小于4°。由于水星的軌道太小，偏不偏心其實(shí)都不太看得出來，所以就當(dāng)它是正圓也無(wú)所謂了，最后得圖

在這里插入圖片描述

au,G,RE,ME = 1.48e11,6.67e-11,1.48e11,5.965e24
m = np.array([3.32e5,0.055,0.815,1,
              0.107,317.8,95.16,14.54,17.14])*ME*6.67e-11
r = np.array([0,0.387,0.723,1,1.524,5.203,
              9.537,19.19,30.7])*RE
theta = np.random.rand(9)*np.pi*2
x = r*np.cos(theta)
y = r*np.sin(theta)
v = np.array([0,47.89,35.03,29.79,
              24.13,13.06,9.64,6.81,5.43])*1000
u = -v*np.sin(theta)
v = v*np.cos(theta)
name = "solar.gif"
fig = plt.figure(figsize=(10,10))
ax = fig.add_subplot(xlim=(-31*RE,31*RE),ylim=(-31*RE,31*RE))
ax.grid()
traces = [ax.plot([],[],'-', lw=0.5)[0] for _ in range(9)]
pts = [ax.plot([],[],marker='o')[0] for _ in range(9)]
k_text = ax.text(0.05,0.85,'',transform=ax.transAxes)
textTemplate = 't = %.3f days\n'
N = 500
dt = 3600*50
ts =  np.arange(0,N*dt,dt)
xs,ys = [],[]
for _ in ts:
    x_ij = (x-x.reshape(len(m),1))
    y_ij = (y-y.reshape(len(m),1))
    r_ij = np.sqrt(x_ij**2+y_ij**2)
    for i in range(len(m)):
        for j in range(len(m)):
            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):
    for i in range(9):
        traces[i].set_data(xs[:n,i],ys[:n,i])
        pts[i].set_data(xs[n,i],ys[n,i])
    k_text.set_text(textTemplate % (ts[n]/3600/24))
    return traces+pts+[k_text]
ani = animation.FuncAnimation(fig, animate, 
    range(N), interval=10, blit=True)
plt.show()
ani.save(name)

由于外圈的行星軌道又長(zhǎng)速度又慢，而內(nèi)層的剛好相反，所以這個(gè)圖很難兼顧，觀感上也不太好看。

如果只畫出木星之前的星體，順便加上小行星帶，可能會(huì)好一些。

在這里插入圖片描述