# 行星運動

> 作者：王一哲
> 日期：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)）

"""
 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 中寫法為

    scene.lights = [local_light(pos=vec(0, 0, 0), color=color.white)]

但是在 VPython 7 已經不支援這樣的寫法，改為在太陽的位置上放置一個光源

    lamp = local_light(pos=vec(0,0,0), color=color.white)

4. 產生表示速度、加速度的箭頭。

5. 開啟繪圖視窗，畫星球與太陽的距離、行星速度、加速度與時間關係圖。

<br />

### 物體運動

1. 計算行星的加速度

    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)）

"""
 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 的字典格式資料，語法為

    data = {"key1": value1, "key2": value2, ...}

引號的部分單引號或雙引號皆可，只要有成對使用即可。每個 key 會對應到一筆資料，資料可以是任意的格式，但是每個 key 都要不一樣，不能重複。如果想要呼叫 data 中 key2 對應的資料，語法為

    data["key2"]

以程式中的資料為例

    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)）

"""
 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. 產生行星時，先把要產生的行星名稱存到串列中

    names = ["Mercury", "Venus", "Earth", "Mars"]

第31 ~ 35行，開啟一個名為 planets 的空白串列，再用 for 迴圈將 names 當中的行星名稱取出來指定給變數 name，在使用 sphere 產生球體時丟到儲存資料的字典當中作為 key 值，就可以讀取對應的資料。接下來再用 **planets.append()** 將括號中的物件加到串列 planets 最後面。Python 當中有一個特殊的寫法，可以在一行程式碼中使用 for 迴圈産生串列，所以第31 ~ 35行可以改寫為

    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。用

    planet.a = - G*sun.m / planet.pos.mag2 * planet.pos.norm()

計算加速度，再用

    planet.v += planet.a*dt

更新速度，最後用

    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)）

"""
 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 則語法為：

class 名稱(來源):
    def 函式1(參數):
        return 回傳值1

    def 函式2(參數):
        return 回傳值2

<br />

例如程式 12-4 中自訂的 class 程式碼

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 的語法為

class 名稱:
    def __init__(self, 參數1, 參數2, ...):
        self.參數1 = 參數1
        self.參數2 = 參數2
        ...
        
    def 函式1(參數):
        return 回傳值1 # 若無回傳值則可省略此行

    def 函式2(參數):
        return 回傳值2 # 若無回傳值則可省略此行

<br />

例如程式 12-4 中自訂的 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

<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)）

"""
 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`