行星運動 - HackMD

# 行星運動 > 作者：王一哲 > 日期：2018/4/18 <br /> 克卜勒行星運動定侓(Kepler’s laws of planetary motion)共有以下3條： 1. 第一定律（軌道定律）：所有行星繞太陽公轉的穩定軌道為橢圓形，太陽在其中一個焦點上。 2. 第二定律（等面積速率定律）：行星與太陽連線於單位時間內掃過的面積相等。若太陽與行星連線為 $r$、行星速度為 $v$、$r$ 與 $v$ 的夾角為 $\theta$，則行星與太陽連線於單位時間內掃過的面積為 $$ \frac{\Delta A}{\Delta t} = \frac{1}{2}rv \sin \theta = 定值 $$ 3. 第三定律（週期定律）：所有繞太陽公轉的行星，公轉週期 $T$ 與平均軌道半徑 $a$ 的關係為 $$ \frac{a^3}{T^2} = 定值 $$ <br /> 我們知道第一定律是因為太陽與恆星之間只有重力作用，依照萬有引力定律可以證明只有橢圓是穩定軌道。第二定律則是因為重力通過太陽，所以行星相對於太陽的角動量守恆，因此行星與太陽連線於單位時間內掃過的面積相等。第三定律則可以用有引力定律當作向心力推導出來。如果配合太陽系星球的真實數據，應該可以將這三個定律畫出來。 <img style="display: block; margin-left: auto; margin-right: auto" height="80%" width="80%" src="https://i.imgur.com/8xTWd1U.png"> <div style="text-align:center">用 ejs 自製的行星運動定律動畫</div> <br /> ## 程式 12-1.行星運動, 自訂星球速度、距離, 可更改萬有引力定律中 r 的次方（[取得程式碼](https://github.com/YiZheWangTw/VPythonTutorial/blob/master/12.%E8%A1%8C%E6%98%9F%E9%81%8B%E5%8B%95/12-1_Kepler_1st.py)）　（[GlowScript 網站動畫連結](http://www.glowscript.org/#/user/yizhe/folder/Public/program/12-1Kepler1st)） ```python= """ VPython教學: 12-1.行星運動, 自訂星球速度、距離, 可更改萬有引力定律中 r 的次方 Ver. 1: 2018/2/25 Ver. 2: 2018/3/1 Ver. 3: 2019/9/8 作者: 王一哲 """ from vpython import * """ 1. 參數設定, 設定變數及初始值 """ G = 6.67408E-11 # 重力常數 size = 1E10 # 星球半徑, 放大約1000倍, 否則會看不見 sun_m = 1988500E24 # 太陽質量 d = 1.5E11 # 地球平均軌道半徑為 1.5E11 m v0 = 29290 # 地球於遠日點公轉速率為 29290 m/s theta = radians(90) # 速度與 +x 軸夾角, 用 radians 將單位換為 rad, 預設為 90 度 n = 2 # 萬有引力定律中 r 的次方 t = 0 # 時間 dt = 60*60 # 時間間隔 """ 2. 畫面設定 """ scene = canvas(title="Planetary Motion", width=600, height=600, x=0, y=0, background=color.black) # 產生太陽 sun, 行星 planet sun = sphere(pos=vec(0, 0, 0), radius=size, m=sun_m, color=color.orange, emissive=True) planet = sphere(pos=vec(d, 0, 0), radius=size, texture=textures.earth, make_trail=True, trail_color=color.blue, trail_radius=0.1*size, v=vec(v0*cos(theta), v0*sin(theta), 0)) lamp = local_light(pos=vec(0,0,0), color=color.white) # 產生表示速度、加速度用的箭頭 arrow_v = arrow(pos=planet.pos, axis=vec(0, 0, 0), shaftwidth=0.3*size, color=color.green) arrow_a = arrow(pos=planet.pos, axis=vec(0, 0, 0), shaftwidth=0.3*size, color=color.magenta) # 產生繪圖視窗 gdr = graph(title="r-t plot", width=600, height=450, x=0, y=600, xtitle="<i>t</i> (s)", ytitle="<i>r</i> (m)") gdv = graph(title="v-t plot", width=600, height=450, x=0, y=1050, xtitle="<i>t</i> (s)", ytitle="<i>v</i> (m/s)") gda = graph(title="a-t plot", width=600, height=450, x=0, y=1500, xtitle="<i>t</i> (s)", ytitle="<i>a</i> (m/s<sup>2</sup>") rt = gcurve(graph=gdr, color=color.blue) vt = gcurve(graph=gdv, color=color.green) at = gcurve(graph=gda, color=color.magenta) """ 3. 星球運動部分 """ while True: rate(60*24*2) # 更新行星加速度、速度、位置 planet.a = -G*sun.m / planet.pos.mag**n * planet.pos.norm() planet.v += planet.a*dt planet.pos += planet.v*dt # 更新代表速度、加速度的箭頭位置、方向、長度 arrow_v.pos = planet.pos arrow_a.pos = planet.pos arrow_v.axis = planet.v*2E6 # 乘以 2E6 放大箭頭, 否則會看不見 arrow_a.axis = planet.a*1E13 # 乘以 1E13 放大箭頭, 否則會看不見 # 畫出 r-t, v-t, a-t 圖 rt.plot(pos=(t, planet.pos.mag)) vt.plot(pos=(t, planet.v.mag)) at.plot(pos=(t, planet.a.mag)) # 更新時間 t += dt ``` <br /> ### 參數設定在此設定變數為 G、size、sun_m、d、v0、theta、n、t、dt，用途已寫在該行的註解當中。其中 E 代表乘以 10 的幾次方，例如 1E3 = 1000，也可以用 e 代替。由於我希望盡量使用真實數據，但是當距離為真實數據時，星球半徑就必須放大，否則會看不到星球。 <br /> ### 畫面設定 1. 產生太陽時最後一個選項為 **emissive=True** ，這個選項的預設值為 False，若設定為 True 則太陽會發光。 2. 產生地球時設定材質 **texture=textures.earth** ，這是因為 VPython 有內建地球的貼圖，不需要自己找圖片。但由於地球物件沒有設定顏色，軌跡的顏色預設為白色，若想改藍色就要加上 **trail_color=color.blue** 。若要改變軌跡的半徑為 0.1*size 則要加上 **trail_radius=0.1 \* size** 。 3. 若要使太陽發光，在舊版的 VPython 中寫法為 ```python scene.lights = [local_light(pos=vec(0, 0, 0), color=color.white)] ``` 但是在 VPython 7 已經不支援這樣的寫法，改為在太陽的位置上放置一個光源 ```python lamp = local_light(pos=vec(0,0,0), color=color.white) ``` 4. 產生表示速度、加速度的箭頭。 5. 開啟繪圖視窗，畫星球與太陽的距離、行星速度、加速度與時間關係圖。 <br /> ### 物體運動 1. 計算行星的加速度 ```python planet.a = - G*sun.m / planet.pos.mag**n * planet.pos.norm() ``` 為了要能夠改變萬有引力定律當中距離的次方，因此不使用 **planet.pos.mag2** 計算行星與太陽連線距離的平方，改用 **planet.pos.mag\*\*n** 計算行星與太陽連線距離的 n 次方。接著更新行星的速度、位置。 2. 更新代表速度、加速度的箭頭位置、方向、長度，由於行星與太陽間的距離很大，要箭頭放大才能看到。 3. 最後畫 r - t、v - t、a - t 圖、更新時間。 <br /> ### 模擬結果由於我是用地球的遠日距、位於遠日點時的速度代入程式中，因此預設的角度為90度，可以將角度改為其它的值，例如120度，地球的軌道就會變成很扁的橢圓形。也可以試著將 n 從 2 改成其它數值，只要不是 2 就沒有辦法形成穩定軌道。 <br /> <img style="display: block; margin-left: auto; margin-right: auto" height="80%" width="80%" src="https://i.imgur.com/skIubkQ.png"> <div style="text-align:center">角度為 120 度的模擬結果</div> <br /> <img style="display: block; margin-left: auto; margin-right: auto" height="80%" width="80%" src="https://i.imgur.com/bsrVdIA.png"> <div style="text-align:center">n = 2.015 的模擬結果</div> <br /> <img style="display: block; margin-left: auto; margin-right: auto" height="80%" width="80%" src="https://i.imgur.com/vJHpZkR.png"> <div style="text-align:center">n = 2.02</div> <br /> ## 程式 12-2.行星運動, 用 dictionary 儲存星球資料（[取得程式碼](https://github.com/YiZheWangTw/VPythonTutorial/blob/master/12.%E8%A1%8C%E6%98%9F%E9%81%8B%E5%8B%95/12-2_Kepler.py)）　（[GlowScript 網站動畫連結](http://www.glowscript.org/#/user/yizhe/folder/Public/program/12-2Kepler)） ```python= """ VPython教學: 12-2.行星運動, 用 dictionary 儲存星球資料 Ver. 1: 2018/2/25 Ver. 2: 2019/9/8 作者: 王一哲 """ from vpython import * """ 1. 參數設定, 設定變數及初始值, 太陽及行星半徑、質量、遠日距、遠日點速率, 資料來源 """ # 用 dictionary 儲存星球資料, 半徑 radius, 質量 mass, 遠日距 d_at_aphelion, 於遠日點的速率 v_at_aphelion radius = {"Mercury": 2439700, "Venus": 6051800, "Earth": 6371000, "Mars": 3389500, "Sun": 696392000} mass = {"Mercury": 0.33011E24, "Venus": 4.8675E24, "Earth": 5.9723E24, "Mars": 0.64171E24, "Sun": 1988500E24} d_at_aphelion = {"Mercury": 6982E7, "Venus": 10894E7, "Earth": 15210E7, "Mars": 24923E7} v_at_aphelion = {"Mercury": 38860, "Venus": 34790, "Earth": 29290, "Mars": 21970} G = 6.67408E-11 # 重力常數 eps = 10000 # 精準度 t = 0 # 時間 dt = 60*60 # 時間間隔 """ 2. 畫面設定 """ scene = canvas(title="Planetary Motion", width=600, height=600, x=0, y=0, background=color.black) # 產生太陽 sun, 地球 earth 及火星 mars sun = sphere(pos=vec(0,0,0), radius=radius["Sun"]*20, m=mass["Sun"], color=color.orange, emissive=True) earth = sphere(pos=vec(d_at_aphelion["Earth"], 0, 0), radius=radius["Earth"]*2E3, m=mass["Earth"], texture=textures.earth, make_trail=True, trail_color=color.blue, retain=365, v=vec(0, v_at_aphelion["Earth"], 0)) mars = sphere(pos=vec(d_at_aphelion["Mars"], 0, 0), radius=radius["Mars"]*2E3, m=mass["Mars"], color=color.red, make_trail=True, retain=365, v=vec(0, v_at_aphelion["Mars"], 0)) lamp = local_light(pos=vec(0, 0, 0), color=color.white) """ 3. 星球運動部分 """ while True: rate(60*24) # 更新行星加速度、速度、位置 earth.a = -G*sun.m / earth.pos.mag2 * earth.pos.norm() earth.v += earth.a*dt earth.pos += earth.v*dt mars.a = -G*sun.m / mars.pos.mag2 * mars.pos.norm() mars.v += mars.a*dt mars.pos += mars.v*dt # 判斷行星是否回到遠日點, 若回到出發點則顯示經過的時間 if abs(earth.pos.x - d_at_aphelion["Earth"]) <= eps: print("t_Earth =", t) if abs(mars.pos.x - d_at_aphelion["Mars"]) <= eps: print("t_Mars =", t) # 更新時間 t += dt ``` <br /> ### 參數設定 1. 字典(dictionary)是 Python 特殊的儲存資料格式，假設定義一個名稱為 data 的字典格式資料，語法為 ```python data = {"key1": value1, "key2": value2, ...} ``` 引號的部分單引號或雙引號皆可，只要有成對使用即可。每個 key 會對應到一筆資料，資料可以是任意的格式，但是每個 key 都要不一樣，不能重複。如果想要呼叫 data 中 key2 對應的資料，語法為 ```python data["key2"] ``` 以程式中的資料為例 ```python radius = {"Mercury": 2439700, "Venus": 6051800, "Earth": 6371000, "Mars": 3389500, "Sun": 696392000} ``` 若輸入 radius["Sun"] ，則系統會輸出 696392000。 2. 在此設定變數還有 G、eps、t、dt，用途已寫在該行的註解當中。 <br /> ### 畫面設定畫面設定的部分與程式 12-1 非常相似，但是在產生星球物件時有點不一樣，之前是直接輸入對應的變數名稱，在此則改用儲存資料用的字典 + 名稱的方式來處理，例如設定地球位於遠日點時是使用 **pos = vec(d_at_aphelion["Earth"], 0, 0)** ，這樣的好處在於產生每個星球時所用的程式碼幾乎一模一樣，只是修改了星球的名字，就能自動找到對應的資料。 <br /> ### 物體運動物體運動的部分與程式 12-1 非常相似，由於只畫出兩個行星，只要將計算加速度、更新速度、位移的程式碼針對兩個行星各寫一次即可。另外再用 if 判斷行星是否回到出發點附近，用來計算週期。雖然這樣的寫法不算精準，但至少還堪用。 <br /> ### 模擬結果由模擬程式得到的地球公轉週期為 31550400 s，約為 365.17 天，與回歸年 365.2421990741 天相差不到 0.02%。由模擬程式得到的火星公轉週期為 59335200 s，為 686.75 天，約為 1.88 年，火星的回歸週期實際值為 686.973 天，相差約 0.03%。由地球、火星的模擬結果來看，這個程式應該相當準確。 <br /> <img style="display: block; margin-left: auto; margin-right: auto" height="80%" width="80%" src="https://i.imgur.com/9f81GIB.png"> <div style="text-align:center">程式 12-2 畫面截圖</div> <br /> ## 程式 12-3.行星運動, 用 dictionary 儲存星球資料, 用 for 迴圈產生行星（[取得程式碼](https://github.com/YiZheWangTw/VPythonTutorial/blob/master/12.%E8%A1%8C%E6%98%9F%E9%81%8B%E5%8B%95/12-3_Kepler_for.py)）　（[GlowScript 網站動畫連結](http://www.glowscript.org/#/user/yizhe/folder/Public/program/12-3Keplerfor)） ```python= """ VPython教學: 12-3.行星運動, 用 dictionary 儲存星球資料, 用 for 迴圈產生行星 Ver. 1: 2018/2/26 Ver. 2: 2019/9/8 作者: 王一哲 """ from vpython import * """ 1. 參數設定, 設定變數及初始值, 太陽及行星半徑、質量、遠日距、遠日點速率, 資料來源 """ # 用 dictionary 儲存星球資料, 半徑 radius, 質量 mass, 遠日距 d_at_aphelion, 於遠日點的速率 v_at_aphelion radius = {"Mercury": 2439700, "Venus": 6051800, "Earth": 6371000, "Mars": 3389500, "Sun": 696392000} mass = {"Mercury": 0.33011E24, "Venus": 4.8675E24, "Earth": 5.9723E24, "Mars": 0.64171E24, "Sun": 1988500E24} material = {"Mercury": color.cyan, "Venus": color.yellow, "Earth": color.blue, "Mars": color.red, "Sun": color.orange} d_at_aphelion = {"Mercury": 6982E7, "Venus": 10894E7, "Earth": 15210E7, "Mars": 24923E7} v_at_aphelion = {"Mercury": 38860, "Venus": 34790, "Earth": 29290, "Mars": 21970} G = 6.67408E-11 # 重力常數 eps = 10000 # 精準度 t = 0 # 時間 dt = 60*60 # 時間間隔 """ 2. 畫面設定 """ scene = canvas(title="Planetary Motion", width=600, height=600, x=0, y=0, background=color.black) # 產生太陽 sun sun = sphere(pos=vec(0, 0, 0), radius=radius["Sun"]*20, m=mass["Sun"], color=color.orange, emissive=True) # 用 for 迴圈產生水星、金星、地球、火星 names = ["Mercury", "Venus", "Earth", "Mars"] planets = [] for name in names: planets.append(sphere(pos=vec(d_at_aphelion[name], 0, 0), radius=radius[name]*2E3, m=mass[name], color=material[name], make_trail=True, retain = 365, v=vec(0, v_at_aphelion[name], 0))) lamp = local_light(pos=vec(0, 0, 0), color=color.white) """ 3. 星球運動部分 """ while True: rate(60*24) # 用 for 迴圈自動跑完所有行星的資料 for planet in planets: planet.a = -G*sun.m / planet.pos.mag2 * planet.pos.norm() planet.v += planet.a*dt planet.pos += planet.v*dt # 更新時間 t += dt ``` <br /> ### 程式設計部分這是以程式 12-2 為基礎衍生出來的偷懶寫法，因為對每個行星而言，產生行星與計算加速度、更新速度、位置的程式碼幾乎一模一樣，但如果總共有4個行星的話，同樣的程式碼要寫4次，這樣實在太麻煩了，如果使用 for 迴圈就能少寫 3/4 的程式碼。 1. 產生行星時，先把要產生的行星名稱存到串列中 ```python names = ["Mercury", "Venus", "Earth", "Mars"] ``` 第31 ~ 35行，開啟一個名為 planets 的空白串列，再用 for 迴圈將 names 當中的行星名稱取出來指定給變數 name，在使用 sphere 產生球體時丟到儲存資料的字典當中作為 key 值，就可以讀取對應的資料。接下來再用 **planets.append()** 將括號中的物件加到串列 planets 最後面。Python 當中有一個特殊的寫法，可以在一行程式碼中使用 for 迴圈産生串列，所以第31 ~ 35行可以改寫為 ```python planets = [sphere(pos=vec(d_at_aphelion[name], 0, 0), radius=radius[name]*2E3, m=mass[name], color=material[name], make_trail=True, retain = 365, v=vec(0, v_at_aphelion[name], 0)) for name in names] ``` <br /> 2. 計算行星的加速度、更新速度、位置。用 for 迴圈將 planets 當中的行星物件取出來指定給變數 planet。用 ```python planet.a = - G*sun.m / planet.pos.mag2 * planet.pos.norm() ``` 計算加速度，再用 ```python planet.v += planet.a*dt ``` 更新速度，最後用 ```python planet.pos += planet.v*dt ``` 更新位置。 <br /> ### 模擬結果我們很成功地使用 for 迴圈減少程式碼，並且畫出水星、金星、地球、火星繞太陽公轉的動畫，如果只想畫出其中幾個行星，只要修改串列 names 即可。 <br /> <img style="display: block; margin-left: auto; margin-right: auto" height="80%" width="80%" src="https://i.imgur.com/Oc5mhe6.png"> <div style="text-align:center">程式 12-3 畫面截圖</div> <br /> ## 程式 12-4.行星運動, 用 dictionary 儲存星球資料, 用 class 產生行星（[取得程式碼](https://github.com/YiZheWangTw/VPythonTutorial/blob/master/12.%E8%A1%8C%E6%98%9F%E9%81%8B%E5%8B%95/12-3_Kepler_for.py)）　（[GlowScript 網站動畫連結](http://www.glowscript.org/#/user/yizhe/folder/Public/program/12-4Keplerclass)） ```python= """ VPython教學: 12-4.行星運動, 用 dictionary 儲存星球資料, 用 class 產生行星 Ver. 1: 2018/2/26 Ver. 2: 2018/10/25 改為不用繼承的 class Ver. 3: 2019/9/8 作者: 王一哲 """ from vpython import * """ 1. 參數設定, 設定變數及初始值, 太陽及行星半徑、質量、遠日距、遠日點速率, 資料來源 """ # 用 dictionary 儲存星球資料, 半徑 radius, 質量 mass, 遠日距 d_at_aphelion, 於遠日點的速率 v_at_aphelion radius = {"Mercury": 2439700, "Venus": 6051800, "Earth": 6371000, "Mars": 3389500, "Sun": 696392000} mass = {"Mercury": 0.33011E24, "Venus": 4.8675E24, "Earth": 5.9723E24, "Mars": 0.64171E24, "Sun": 1988500E24} material = {"Mercury": color.cyan, "Venus": color.yellow, "Earth": color.blue, "Mars": color.red, "Sun": color.orange} d_at_aphelion = {"Mercury": 6982E7, "Venus": 10894E7, "Earth": 15210E7, "Mars": 24923E7} v_at_aphelion = {"Mercury": 38860, "Venus": 34790, "Earth": 29290, "Mars": 21970} G = 6.67408E-11 # 重力常數 t = 0 # 時間 dt = 60*60 # 時間間隔 """ 2. 產生行星類別 """ # 使用繼承的 class 寫法 ''' class Planet(sphere): sun_m = mass["Sun"] def a(self): return - G*self.sun_m / self.pos.mag2 * self.pos.norm() ''' # 不用繼承的 class 寫法 class Planet: def __init__(self, pos, radius, mass, color, v): self.pos = pos self.radius = radius self.mass = mass self.color = color self.v = v self.a = 0 self.planet = sphere(pos=self.pos, radius=self.radius, mass=self.mass, color=self.color, make_trail=True, retain=365, v=self.v) def update(self, dt): self.dt = dt self.a = -G*mass["Sun"] / self.planet.pos.mag2 * self.planet.pos.norm() self.v += self.a * self.dt self.planet.pos += self.v * self.dt """ 3. 畫面設定 """ scene = canvas(title="Planetary Motion", width=600, height=600, x=0, y=0, background=color.black) sun = sphere(pos=vec(0,0,0), radius=radius["Sun"]*20, color=color.orange, emissive=True) lamp = local_light(pos=vec(0, 0, 0), color=color.white) # 用 for 迴圈產生水星、金星、地球、火星 names = ["Mercury", "Venus", "Earth", "Mars"] planets = [] for name in names: planets.append(Planet(pos=vec(d_at_aphelion[name], 0, 0), radius=radius[name]*2E3, mass=mass[name], color=material[name], v=vec(0, v_at_aphelion[name], 0))) # 第58 ~ 62行另一個寫法 #planets = [Planet(pos=vec(d_at_aphelion[name], 0, 0), radius=radius[name]*2E3, mass=mass[name], # color=material[name], v=vec(0, v_at_aphelion[name], 0)) for name in names] """ 4. 星球運動部分 """ while True: rate(60*24) # 用 for 迴圈自動跑完所有行星的資料 for planet in planets: planet.update(dt) # 更新時間 t += dt ``` <br /> ### 程式設計部分這是以程式 12-3 為基礎，再加上 **類別 (class)** 的寫法。如果想要在 Python 程式中產生許多相似的物件，而且我們還想要在這些物件中自訂函式，就可以使用 class 來達成目的，通常 class 的名稱開頭為大寫字母。class 的語法有兩種，如果繼承其它的 class 則語法為： ```python class 名稱(來源): def 函式1(參數): return 回傳值1 def 函式2(參數): return 回傳值2 ``` <br /> 例如程式 12-4 中自訂的 class 程式碼 ```python class Planet(sphere): sun_m = mass["Sun"] def a(self): return - G*self.sun_m / self.pos.mag2 * self.pos.norm() ``` <br /> class 名稱為 Planet，這個類別是由 sphere 衍生出來的，繼承了 sphere 所有的屬性，因此使用 Planet 產生的物件可以定義 pos、radius……等屬性。除此之外，還自行定義了函式 a，用來計算並回傳此行星物件的加速度。 <br /> 如果不使用繼承，則 class 的語法為 ```python class 名稱: def __init__(self, 參數1, 參數2, ...): self.參數1 = 參數1 self.參數2 = 參數2 ... def 函式1(參數): return 回傳值1 # 若無回傳值則可省略此行 def 函式2(參數): return 回傳值2 # 若無回傳值則可省略此行 ``` <br /> 例如程式 12-4 中自訂的 class 程式碼 ```python class Planet: def __init__(self, pos, radius, mass, color, v): self.pos = pos self.radius = radius self.mass = mass self.color = color self.v = v self.a = 0 self.planet = sphere(pos=self.pos, radius=self.radius, mass=self.mass, color=self.color, make_trail=True, retain=365, v=self.v) def update(self, dt): self.dt = dt self.a = -G*mass["Sun"] / self.planet.pos.mag2 * self.planet.pos.norm() self.v += self.a * self.dt self.planet.pos += self.v * self.dt ``` <br /> 使用 Planet 建立物件時，需要輸入 pos、radius、mass、color、v 等5個參數，並將參數值代入自訂函式 \_\_init\_\_，産生一個名為 planet 的 sphere 物件。使用另一個自訂函式 update 時，需要輸入參數 dt，函式可以計算 sphere 物件的加速度、速度，更新 planet 的位置。關於 class 更詳細的說明請參閱這篇文章：〈[自訂類別 Class](https://hackmd.io/@yizhewang/r1XyZogto)〉。 <br /> ### 模擬結果結果與程式 12-3 完全相同。 <br /> ## 程式 12-5.克卜勒第三行星運動定律（[取得程式碼](https://github.com/YiZheWangTw/VPythonTutorial/blob/master/12.%E8%A1%8C%E6%98%9F%E9%81%8B%E5%8B%95/12-5_Kepler_3rd.py)）　（[GlowScript 網站動畫連結](http://www.glowscript.org/#/user/yizhe/folder/Public/program/12-5Kepler3rd)） ```python= """ VPython教學: 12-5.克卜勒第三行星運動定律 Ver. 1: 2018/3/1 Ver. 2: 2018/4/17 Ver. 3: 3019/9/8 作者: 王一哲 """ from vpython import * """ 1. 參數設定, 設定變數及初始值, 太陽及行星半徑、質量、遠日距、遠日點速率, 資料來源太陽 https://nssdc.gsfc.nasa.gov/planetary/factsheet/sunfact.html 地球 https://nssdc.gsfc.nasa.gov/planetary/factsheet/earthfact.html """ size = 2E10 # 星球半徑, 放大約2000倍, 否則會看不見 sun_m = 1988500E24 # 太陽質量 d = 1.5E11 # 地球平均軌道半徑為 1.5E11 m v0 = 29780 # 地球公轉平均速率為 29780 m/s eps = 10000 # 計算週期用的精準度 n = 4 # 自訂行星 planet 軌道半徑為 n*d ec = 0.0 # 自訂行星軌道離心率(eccentricity), 若 n != 1 時設為 0 G = 6.67408E-11 # 重力常數 t = 0 # 時間 dt = 60*60 # 時間間隔 """ 2. 畫面設定 (1) 用 canvas 物件作為顯示動畫用的視窗 http://www.glowscript.org/docs/VPythonDocs/canvas.html (2) 用 sphere 物件產生星球 http://www.glowscript.org/docs/VPythonDocs/sphere.html (3) 星球的半徑要手動調整比例, 否則會看不到星球 """ scene = canvas(title="Kepler's Third Law of Planetary Motion", width=600, height=600, x=0, y=0, background=color.black) # 產生太陽 sun, 地球 earth 及自訂行星 planet sun = sphere(pos=vec(0,0,0), radius=size, m=sun_m, color=color.orange, emissive=True) earth = sphere(pos=vec(d, 0, 0), radius=size, texture=textures.earth, make_trail=True, trail_color=color.blue, retain=365, v=vec(0, v0, 0)) planet = sphere(pos=vec(n*(1+ec)*d, 0, 0), radius=size, color=color.red, make_trail=True, retain=365*n, v=vec(0, v0/sqrt(n)*sqrt((1-ec)/(1+ec)), 0)) line = cylinder(pos=vec(0, 0, 0), axis=vec(n*(1+ec)*d, 0, 0), radius=0.3*size, color=color.yellow) # 原來的寫法為 scene.lights = [local_light(pos = vec(0,0,0), color = color.white)] # 在 VPython 7 中 canvas.lights 無法設定為 local_light, 只能另外在太陽處放置另一個光源 lamp lamp = local_light(pos=vec(0,0,0), color=color.white) """ 3. 星球運動部分 """ while True: rate(60*24) # 計算行星加速度、更新速度、位置 earth.a = -G*sun.m / earth.pos.mag2 * earth.pos.norm() earth.v += earth.a*dt earth.pos += earth.v*dt planet.a = -G*sun.m / planet.pos.mag2 * planet.pos.norm() planet.v += planet.a*dt planet.pos += planet.v*dt # 判斷行星是否回到出發點 if abs(earth.pos.x - d) <= eps: print("t_Earth =", t) if abs(planet.pos.x - n*(1+ec)*d) <= eps: print("t_planet =", t) # 更新時間 t += dt ``` <br /> ### 程式設計部分這是以程式 12-1 為基礎，除了地球以外再加上一個自訂行星，用來說明克卜勒行三行星運動定律。假設自訂行星的平均軌道半徑為 **n * d**，離心率為 **ec**，則遠日距為 **n * (1 + ec) * d**。若地球於遠日點時的速率為 **v0**，則自訂行星於遠日點時的速率為 **v0 / sqrt(n) * sqrt((1 - ec) / (1 + ec))**。剩下的部分幾乎沒有改變。地球的速度推導過程如下：假設太陽的質量為 \$ M \$，地球的質量為 \$ m \$，地球的平均軌道半徑為 \$ d \$，地球地球的速率為 \$ v_0 \$，由太陽對地球的重力當作向心力可得 $$ \frac{GMm}{d^2} = m \cdot \frac{v_0^2}{d} ~\Rightarrow~ v_0 = \sqrt{\frac{GM}{d}} $$ <br /> 自訂行星位於遠日點的速率推導過程如下：假設太陽的質量為 \$ M \$，行星的質量為 \$ m \$，行星的平均軌道半徑為 $a = nd$，行星位於近日點的速率為 \$ v_1 \$ ，位於遠日點的速率為 \$ v_2 \$，由角動量守恆可得 $$ L = (1-e)amv_1 = (1+e)amv_2 ~\Rightarrow~ v_1 = \frac{1+e}{1-e}v_2 $$ 由力學能守恆可得 $$ \frac{1}{2}mv_1^2 + \left[ -\frac{GMm}{(1-e)a} \right ] = \frac{1}{2}mv_2^2 + \left[ -\frac{GMm}{(1+e)a} \right ] $$ 將 $v_1$、$v_2$ 的關係代入上式 $$ \left(\frac{1+e}{1-e} \right )^2 v_2^2 -\frac{2GM}{(1-e)a} = v_2^2 -\frac{2GM}{(1+e)a} $$ $$ \left[\left(\frac{1+e}{1-e} \right )^2 - 1 \right ] v_2^2 = \left( \frac{1}{1-e} - \frac{1}{1+e} \right) \frac{2GM}{a} $$ $$ \frac{4ev_2^2}{(1-e)^2} = \frac{4eGM}{(1-e)(1+e)a} $$ $$ v_2 = \sqrt{\frac{1+e}{1-e} \cdot \frac{GM}{a}} = \sqrt{\frac{1+e}{1-e} \cdot \frac{GM}{nd}} = \sqrt{\frac{1+e}{(1-e)n}} v_0 $$ <br /> ### 模擬結果 1. 假設 n = 4，則自訂行星的週期為 8 年。 2. 假設 n = 1，調整 0 < ec < 1，則自訂行星的週期為 1 年。 <img style="display: block; margin-left: auto; margin-right: auto" height="80%" width="80%" src="https://i.imgur.com/WWfvhLD.png"> <div style="text-align:center">n = 4、ec = 0，自訂行星的週期為 8 年</div> <br /> <img style="display: block; margin-left: auto; margin-right: auto" height="80%" width="80%" src="https://i.imgur.com/VCXKG4Q.png"> <div style="text-align:center">n = 1、ec = 0.5，自訂行星的週期為 1 年</div> <br /> ## 結語用 VPython 可以模擬出各種不同的狀況，例如修改萬有引力定律的型式，或是用實際數據代入模型中看看是否符合實際情形，同學們可以試著將程式修改成自己想要的樣子。 <br /> ## 太陽系天體資料來源 1. **太陽** https://nssdc.gsfc.nasa.gov/planetary/factsheet/sunfact.html 2. **水星** https://nssdc.gsfc.nasa.gov/planetary/factsheet/mercuryfact.html 3. **金星** https://nssdc.gsfc.nasa.gov/planetary/factsheet/venusfact.html 4. **地球** https://nssdc.gsfc.nasa.gov/planetary/factsheet/earthfact.html 5. **火星** https://nssdc.gsfc.nasa.gov/planetary/factsheet/marsfact.html <br /> ## VPython官方說明書 1. **canvas**: http://www.glowscript.org/docs/VPythonDocs/canvas.html 2. **sphere**: http://www.glowscript.org/docs/VPythonDocs/sphere.html 3. **cylinder**: http://www.glowscript.org/docs/VPythonDocs/cylinder.html 4. **graph**: http://www.glowscript.org/docs/VPythonDocs/graph.html --- ###### tags:`VPython`