行星運動

作者：王一哲
日期：2018/4/18

克卜勒行星運動定侓(Kepler’s laws of planetary motion)共有以下3條：

第一定律（軌道定律）：所有行星繞太陽公轉的穩定軌道為橢圓形，太陽在其中一個焦點上。
第二定律（等面積速率定律）：行星與太陽連線於單位時間內掃過的面積相等。若太陽與行星連線為
$r$ 、行星速度為
$v$ 、
$r$ 與
$v$ 的夾角為
$θ$ ，則行星與太陽連線於單位時間內掃過的面積為

$\frac{Δ A}{Δ t} = \frac{1}{2} r v \sin θ = 定值$
第三定律（週期定律）：所有繞太陽公轉的行星，公轉週期
$T$ 與平均軌道半徑
$a$ 的關係為

$\frac{a^{3}}{T^{2}} = 定值$

我們知道第一定律是因為太陽與恆星之間只有重力作用，依照萬有引力定律可以證明只有橢圓是穩定軌道。第二定律則是因為重力通過太陽，所以行星相對於太陽的角動量守恆，因此行星與太陽連線於單位時間內掃過的面積相等。第三定律則可以用有引力定律當作向心力推導出來。如果配合太陽系星球的真實數據，應該可以將這三個定律畫出來。

Image Not Showing Possible Reasons

The image file may be corrupted
The server hosting the image is unavailable
The image path is incorrect
The image format is not supported

Learn More →

用 ejs 自製的行星運動定律動畫

程式 12-1.行星運動, 自訂星球速度、距離, 可更改萬有引力定律中 r 的次方（取得程式碼）　（GlowScript 網站動畫連結）






























































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

參數設定

在此設定變數為 G、size、sun_m、d、v0、theta、n、t、dt，用途已寫在該行的註解當中。其中 E 代表乘以 10 的幾次方，例如 1E3 = 1000，也可以用 e 代替。由於我希望盡量使用真實數據，但是當距離為真實數據時，星球半徑就必須放大，否則會看不到星球。

畫面設定

產生太陽時最後一個選項為 emissive=True ，這個選項的預設值為 False，若設定為 True 則太陽會發光。
產生地球時設定材質 texture=textures.earth ，這是因為 VPython 有內建地球的貼圖，不需要自己找圖片。但由於地球物件沒有設定顏色，軌跡的顏色預設為白色，若想改藍色就要加上 trail_color=color.blue 。若要改變軌跡的半徑為 0.1*size 則要加上 trail_radius=0.1 * size 。
若要使太陽發光，在舊版的 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)
```
產生表示速度、加速度的箭頭。
開啟繪圖視窗，畫星球與太陽的距離、行星速度、加速度與時間關係圖。

物體運動

計算行星的加速度
```
planet.a = - G*sun.m / planet.pos.mag**n * planet.pos.norm()
```
為了要能夠改變萬有引力定律當中距離的次方，因此不使用 planet.pos.mag2 計算行星與太陽連線距離的平方，改用 planet.pos.mag**n 計算行星與太陽連線距離的 n 次方。接著更新行星的速度、位置。
更新代表速度、加速度的箭頭位置、方向、長度，由於行星與太陽間的距離很大，要箭頭放大才能看到。
最後畫 r - t、v - t、a - t 圖、更新時間。

模擬結果

由於我是用地球的遠日距、位於遠日點時的速度代入程式中，因此預設的角度為90度，可以將角度改為其它的值，例如120度，地球的軌道就會變成很扁的橢圓形。也可以試著將 n 從 2 改成其它數值，只要不是 2 就沒有辦法形成穩定軌道。

Image Not Showing Possible Reasons

The image file may be corrupted
The server hosting the image is unavailable
The image path is incorrect
The image format is not supported

Learn More →

角度為 120 度的模擬結果

Image Not Showing Possible Reasons

The image file may be corrupted
The server hosting the image is unavailable
The image path is incorrect
The image format is not supported

Learn More →

n = 2.015 的模擬結果

Image Not Showing Possible Reasons

The image file may be corrupted
The server hosting the image is unavailable
The image path is incorrect
The image format is not supported

Learn More →

n = 2.02

程式 12-2.行星運動, 用 dictionary 儲存星球資料（取得程式碼）　（GlowScript 網站動畫連結）





















































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

參數設定

字典(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。
在此設定變數還有 G、eps、t、dt，用途已寫在該行的註解當中。

畫面設定

畫面設定的部分與程式 12-1 非常相似，但是在產生星球物件時有點不一樣，之前是直接輸入對應的變數名稱，在此則改用儲存資料用的字典 + 名稱的方式來處理，例如設定地球位於遠日點時是使用 pos = vec(d_at_aphelion["Earth"], 0, 0) ，這樣的好處在於產生每個星球時所用的程式碼幾乎一模一樣，只是修改了星球的名字，就能自動找到對應的資料。

物體運動

物體運動的部分與程式 12-1 非常相似，由於只畫出兩個行星，只要將計算加速度、更新速度、位移的程式碼針對兩個行星各寫一次即可。另外再用 if 判斷行星是否回到出發點附近，用來計算週期。雖然這樣的寫法不算精準，但至少還堪用。

模擬結果

由模擬程式得到的地球公轉週期為 31550400 s，約為 365.17 天，與回歸年 365.2421990741 天相差不到 0.02%。由模擬程式得到的火星公轉週期為 59335200 s，為 686.75 天，約為 1.88 年，火星的回歸週期實際值為 686.973 天，相差約 0.03%。由地球、火星的模擬結果來看，這個程式應該相當準確。

Image Not Showing Possible Reasons

The image file may be corrupted
The server hosting the image is unavailable
The image path is incorrect
The image format is not supported

Learn More →

程式 12-2 畫面截圖

程式 12-3.行星運動, 用 dictionary 儲存星球資料, 用 for 迴圈產生行星（取得程式碼）　（GlowScript 網站動畫連結）


















































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

程式設計部分

這是以程式 12-2 為基礎衍生出來的偷懶寫法，因為對每個行星而言，產生行星與計算加速度、更新速度、位置的程式碼幾乎一模一樣，但如果總共有4個行星的話，同樣的程式碼要寫4次，這樣實在太麻煩了，如果使用 for 迴圈就能少寫 3/4 的程式碼。

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

計算行星的加速度、更新速度、位置。用 for 迴圈將 planets 當中的行星物件取出來指定給變數 planet。用
```
planet.a = - G*sun.m / planet.pos.mag2 * planet.pos.norm()
```
計算加速度，再用
```
planet.v += planet.a*dt
```
更新速度，最後用
```
planet.pos += planet.v*dt
```
更新位置。

模擬結果

我們很成功地使用 for 迴圈減少程式碼，並且畫出水星、金星、地球、火星繞太陽公轉的動畫，如果只想畫出其中幾個行星，只要修改串列 names 即可。

Image Not Showing Possible Reasons

The image file may be corrupted
The server hosting the image is unavailable
The image path is incorrect
The image format is not supported

Learn More →

程式 12-3 畫面截圖

程式 12-4.行星運動, 用 dictionary 儲存星球資料, 用 class 產生行星（取得程式碼）　（GlowScript 網站動畫連結）













































































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

程式設計部分

這是以程式 12-3 為基礎，再加上 類別 (class) 的寫法。如果想要在 Python 程式中產生許多相似的物件，而且我們還想要在這些物件中自訂函式，就可以使用 class 來達成目的，通常 class 的名稱開頭為大寫字母。class 的語法有兩種，如果繼承其它的 class 則語法為：

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

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

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

class Planet(sphere):
    sun_m = mass["Sun"]
    def a(self):
        return - G*self.sun_m / self.pos.mag2 * self.pos.norm()

class 名稱為 Planet，這個類別是由 sphere 衍生出來的，繼承了 sphere 所有的屬性，因此使用 Planet 產生的物件可以定義 pos、radius……等屬性。除此之外，還自行定義了函式 a，用來計算並回傳此行星物件的加速度。

如果不使用繼承，則 class 的語法為

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

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

例如程式 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

使用 Planet 建立物件時，需要輸入 pos、radius、mass、color、v 等5個參數，並將參數值代入自訂函式 __init__，産生一個名為 planet 的 sphere 物件。使用另一個自訂函式 update 時，需要輸入參數 dt，函式可以計算 sphere 物件的加速度、速度，更新 planet 的位置。關於 class 更詳細的說明請參閱這篇文章：〈自訂類別 Class〉。

模擬結果

結果與程式 12-3 完全相同。

程式 12-5.克卜勒第三行星運動定律（取得程式碼）　（GlowScript 網站動畫連結）































































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

程式設計部分

這是以程式 12-1 為基礎，除了地球以外再加上一個自訂行星，用來說明克卜勒行三行星運動定律。假設自訂行星的平均軌道半徑為 n * d，離心率為 ec，則遠日距為 n * (1 + ec) * d。若地球於遠日點時的速率為 v0，則自訂行星於遠日點時的速率為 v0 / sqrt(n) * sqrt((1 - ec) / (1 + ec))。剩下的部分幾乎沒有改變。

地球的速度推導過程如下：

假設太陽的質量為

M

，地球的質量為

m

，地球的平均軌道半徑為

d

，地球地球的速率為

v_{0}

，由太陽對地球的重力當作向心力可得

\frac{G M m}{d^{2}} = m \cdot \frac{v_{0}^{2}}{d} \Rightarrow v_{0} = \sqrt{\frac{G M}{d}}

自訂行星位於遠日點的速率推導過程如下：

假設太陽的質量為

M

，行星的質量為

m

，行星的平均軌道半徑為

a = n d

，行星位於近日點的速率為

v_{1}

，位於遠日點的速率為

v_{2}

，由角動量守恆可得

L = (1 - e) a m v_{1} = (1 + e) a m v_{2} \Rightarrow v_{1} = \frac{1 + e}{1 - e} v_{2}

由力學能守恆可得

\frac{1}{2} m v_{1}^{2} + [- \frac{G M m}{(1 - e) a}] = \frac{1}{2} m v_{2}^{2} + [- \frac{G M m}{(1 + e) a}]

將

v_{1}

、

v_{2}

的關係代入上式

{(\frac{1 + e}{1 - e})}^{2} v_{2}^{2} - \frac{2 G M}{(1 - e) a} = v_{2}^{2} - \frac{2 G M}{(1 + e) a}

[{(\frac{1 + e}{1 - e})}^{2} - 1] v_{2}^{2} = (\frac{1}{1 - e} - \frac{1}{1 + e}) \frac{2 G M}{a}

\frac{4 e v_{2}^{2}}{(1 - e)^{2}} = \frac{4 e G M}{(1 - e) (1 + e) a}

v_{2} = \sqrt{\frac{1 + e}{1 - e} \cdot \frac{G M}{a}} = \sqrt{\frac{1 + e}{1 - e} \cdot \frac{G M}{n d}} = \sqrt{\frac{1 + e}{(1 - e) n}} v_{0}

模擬結果

假設 n = 4，則自訂行星的週期為 8 年。
假設 n = 1，調整 0 < ec < 1，則自訂行星的週期為 1 年。

Image Not Showing Possible Reasons

The image file may be corrupted
The server hosting the image is unavailable
The image path is incorrect
The image format is not supported

Learn More →

n = 4、ec = 0，自訂行星的週期為 8 年

Image Not Showing Possible Reasons

The image file may be corrupted
The server hosting the image is unavailable
The image path is incorrect
The image format is not supported

Learn More →

n = 1、ec = 0.5，自訂行星的週期為 1 年

結語

用 VPython 可以模擬出各種不同的狀況，例如修改萬有引力定律的型式，或是用實際數據代入模型中看看是否符合實際情形，同學們可以試著將程式修改成自己想要的樣子。

太陽系天體資料來源

VPython官方說明書

canvas: http://www.glowscript.org/docs/VPythonDocs/canvas.html
sphere: http://www.glowscript.org/docs/VPythonDocs/sphere.html
cylinder: http://www.glowscript.org/docs/VPythonDocs/cylinder.html
graph: http://www.glowscript.org/docs/VPythonDocs/graph.html

行星運動

程式 12-1.行星運動, 自訂星球速度、距離, 可更改萬有引力定律中 r 的次方（取得程式碼） （GlowScript 網站動畫連結）

參數設定

畫面設定

物體運動

模擬結果

程式 12-2.行星運動, 用 dictionary 儲存星球資料（取得程式碼） （GlowScript 網站動畫連結）

參數設定

畫面設定

物體運動

模擬結果

程式 12-3.行星運動, 用 dictionary 儲存星球資料, 用 for 迴圈產生行星（取得程式碼） （GlowScript 網站動畫連結）

程式設計部分

模擬結果

程式 12-4.行星運動, 用 dictionary 儲存星球資料, 用 class 產生行星（取得程式碼） （GlowScript 網站動畫連結）

程式設計部分

模擬結果

程式 12-5.克卜勒第三行星運動定律（取得程式碼） （GlowScript 網站動畫連結）

程式設計部分

模擬結果

結語

太陽系天體資料來源

VPython官方說明書

tags:VPython

Read more

APCS實作題2018年：置物櫃分配

Python 自訂函式

終端速度

Python基本語法

程式 12-1.行星運動, 自訂星球速度、距離, 可更改萬有引力定律中 r 的次方（取得程式碼）　（GlowScript 網站動畫連結）

程式 12-2.行星運動, 用 dictionary 儲存星球資料（取得程式碼）　（GlowScript 網站動畫連結）

程式 12-3.行星運動, 用 dictionary 儲存星球資料, 用 for 迴圈產生行星（取得程式碼）　（GlowScript 網站動畫連結）

程式 12-4.行星運動, 用 dictionary 儲存星球資料, 用 class 產生行星（取得程式碼）　（GlowScript 網站動畫連結）

程式 12-5.克卜勒第三行星運動定律（取得程式碼）　（GlowScript 網站動畫連結）

tags:`VPython`