一個地球人在臺北: 重力場及電場

作者：王一哲
日期：2018/5/7

若在空間中有一個質量為 $M$ 且質量均勻分布的球體，若以 $M$ 的球心為原點，在空間中位置向量為 $M$ 處的重力場為

$$ \vec g = - \frac{GM}{r^2} \hat r $$

上式中的負號帶表重力場方向指向球心。若有多個球體，第 $i$ 個球體質量為 $M_i$、球心位置為 $r_i$，則重力場為

$$ \vec g = \sum \vec g_i = - G \sum \frac{M_i}{(\vec r - \vec r_i)^2} \frac{\vec r - \vec r_i}{|\vec r - \vec r_i|} $$

如果要在黑板上畫出各個位置的重力場強度及方向，這幾乎是不可能的任務，下圖是我畫出來的地球重力場示意圖

地球重力場示意圖

在以下的課程中，我們想要在空間每隔一段距離取一個點，以箭頭的長度及方向來表示該點的重力場。因此我們在程式 18-1 中先練習使用 for 迴圈，在空間中每隔一段距離畫一個箭頭。在程式 18-2 中除了畫出箭頭還要計算此處的重力場，再更新箭頭的長度及方向。如果能夠成功畫出一個球體的重力場，在程式 18-3 中則更進一步地畫出兩個球體，甚至是多個球體建立的重力場。

在畫完重力場之後，我們可以用相同的方法畫出帶電球體在空間中建立的電場，其數學型式為

$$ \vec E = \sum \vec E_i = k \sum \frac{Q_i}{(\vec r - \vec r_i)^2} \frac{\vec r - \vec r_i}{|\vec r - \vec r_i|} $$

我們只要稍微修改程式 18-3 就能畫出兩個帶電球體在空間中建立的電場。

程式 18-1. 用 for 迴圈產生箭頭

取得程式碼
GlowScript 網站動畫連結

"""
 VPython教學: 18-1. 用 for 迴圈產生箭頭
 日期: 2018/3/2
 作者: 王一哲
"""
from vpython import *

N, L = 4, 20

scene = canvas(title="Arrows", width=600, height=600, x=0, y=0, background=color.black)

fields = []
for i in range(N):
    for j in range(N):
        for k in range(N):
            fields.append(arrow(pos=vec(L/N*i - L/2, L/N*j - L/2, L/N*k - L/2),
                                axis=vector(4, 0, 0), radius=1, color=color.green))

這個程式與之前的程式相比十分簡短，只要將箭頭畫完就可以了，因此連產生動畫的部分都沒有。比較特別之處在於使用了 3 層的 for 迴圈，第 1 層迴圈的變數是 i ，用來產生位置的 x 坐標；第 2 層迴圈的變數是 j ，用來產生位置的 y 坐標；第 3 層迴圈的變數是 k ，用來產生位置的 z 坐標。由於使用了 range(N)、N = 4 來產生數值，也就是只會產生 0、1、2、3，所以畫箭頭的位置為 -L/2、-L/4、0、L/4。程式執行的成果如下圖。

程式18-1畫面截圖

程式 18-2. 重力場, 可以改成畫不同的星球

取得程式碼
GlowScript 網站動畫連結

"""
 VPython教學: 18-2. 重力場, 可以改成畫不同的星球
 Ver. 1: 2018/3/2
 Ver. 2: 2019/9/14 箭頭的顏色隨著量值改變
 作者: 王一哲
"""
from vpython import *

"""
 1. 參數設定, 設定變數及初始值, 太陽及行星半徑、質量、遠日距、遠日點速率, 資料來源
"""
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}
G = 6.67408E-11   # 重力常數
N = 5             # 將顯示的空間每邊切成 N 等份

"""
 2. 產生行星類別, 回傳行星產生的重力場
"""
class planet_c(sphere):
    def g(self, pos):
        return -G*self.m / mag2(pos-self.pos) * norm(pos-self.pos)

"""
 3. 畫面設定
"""
# 產生動畫視窗
L = radius["Earth"]*4
scene = canvas(title="Gravitatinoal Field", width=600, height=600, x=0, y=0, background=color.black, range=L)

# 產生地球
earth = planet_c(pos=vec(0, 0, 0), radius=radius["Earth"], m=mass["Earth"], texture=textures.earth)

# 計算畫箭頭的位置, 如果不在地球內則加到串列 locations 當中
locations = []
for i in range(N+1):
    for j in range(N+1):
        for k in range(N+1):
            location = vec(L/N*i - L/2, L/N*j - L/2, L/N*k - L/2)
            if(mag(location - earth.pos) > earth.radius):
                locations.append(location)

# 依序讀取串列 locations 的元素, 在對應的位置產生箭頭
fields = []
for location in locations:
    fields.append(arrow(pos=location, axis=vec(0, 0, 0), color=color.green))

# 更新箭頭的長度及方向, 長度乘以 1E6 才能看見, 記錄重力場強度最大值, 量值接近最大值偏紅色, 量值接近 0 偏綠色
fmax = 0
for field in fields:
    field.axis = earth.g(field.pos)*1E6
    if(field.axis.mag >= fmax): fmax = field.axis.mag

for field in fields:
    field.color = vec(field.axis.mag/fmax, 1 - field.axis.mag/fmax, 0)

由於我們在之前的課程〈行星運動〉當中已經蒐集了不少行星資料，並且在程式 12-4 中學到如何利用 class 產生自訂類別的物件，因此我們可以將程式 12-4 和 18-1 合併，寫出程式 18-2。主要的步驟為：

參數設定
自訂類別 planet_c，並在當中新增自訂的方法(method) g，用來回傳此物件在輸入的位置 pos 處產生的重力場。
產生動畫視窗、地球。
計算畫箭頭的位置，如果不在地球內則加到串列 locations 當中。
依序讀取串列 locations 的元素，在對應的位置產生箭頭。
第 49 ~ 55 行，依序讀取箭頭的位置，計算此處的重力場，更新箭頭的長度及方向。由於重力場的量值與畫面的寬度相比太短，因此需要將重力場乘以 1 × 10⁶ 才能在畫面上看到箭頭。將重力場量值當中的最大值記錄為 fmax，量值越大箭頭的顏色越接近紅色，量值越小箭頭的顏色越接近綠色。

程式18-2畫面截圖

程式 18-3. 重力場, 可以改成畫不同的星球

取得程式碼
GlowScript 網站動畫連結

"""
 VPython教學: 18-3. 重力場, 2個星球, 自行調整距離
 Ver. 1: 2018/3/2
 Ver. 2: 2019/9/14 箭頭的顏色隨著量值改變
 作者: 王一哲
"""
from vpython import *

"""
 1. 參數設定, 設定變數及初始值, 太陽及行星半徑、質量、遠日距、遠日點速率, 資料來源
"""
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}
G = 6.67408E-11   # 重力常數
N = 5             # 將顯示的空間每邊切成 N 等份

"""
 2. 產生行星類別, 回傳行星產生的重力場
"""
class planet_c(sphere):
    def g(self, pos):
        return -G*self.m / mag2(pos-self.pos) * norm(pos-self.pos)

"""
 3. 畫面設定
"""
# 產生動畫視窗
d = 2E7
L = 2*d
scene = canvas(title="Gravitatinoal Field", width=600, height=600, x=0, y=0, background=color.black, range=L)

# 產生地球及火星
earth = planet_c(pos=vec(-d/2, 0, 0), radius=radius["Earth"], m=mass["Earth"], texture=textures.earth)
mars = planet_c(pos=vec(d/2, 0, 0), radius=radius["Mars"], m=mass["Mars"], color=color.red)

# 計算畫箭頭的位置, 如果不在地球或火星內則加到串列 locations 當中
locations = []
for i in range(N+1):
    for j in range(N+1):
        for k in range(N+1):
            location = vec(L/N*i - L/2, L/N*j - L/2, L/N*k - L/2)
            if(mag(location-earth.pos) > earth.radius and mag(location-mars.pos) > mars.radius):
                locations.append(location)

# 依序讀取串列 locations 的元素, 在對應的位置產生箭頭
fields = []
for location in locations:
    fields.append(arrow(pos=location, axis=vec(0, 0, 0), color=color.green))

# 更新箭頭的長度及方向, 長度乘以 1E6 才能看見, 記錄重力場強度最大值, 量值接近最大值偏紅色, 量值接近 0 偏綠色
fmax = 0
for field in fields:
    field.axis = (earth.g(field.pos) + mars.g(field.pos))*1E6
    if(field.axis.mag >= fmax): fmax = field.axis.mag

for field in fields:
    field.color = vec(field.axis.mag/fmax, 1 - field.axis.mag/fmax, 0)

程式 18-3 與 18-2 幾乎一模一樣，只是多畫了火星。由於要產生兩個行星，在此應該就能看出自訂類別較為方便之處。需要特別注意一點，兩個行星之間的距離並未按照真實數據繪製，否則畫面中什麼東西都看不到。

程式18-3畫面截圖

程式 18-4. 電場, 2個球狀帶電體

取得程式碼
GlowScript 網站動畫連結

"""
 VPython教學: 18-4. 電場, 2個球狀帶電體
 Ver. 1: 2018/3/2
 Ver. 2: 2019/9/14 箭頭的顏色隨著量值改變, 改為不用繼承的類別
 作者: 王一哲
"""
from vpython import *

"""
 1. 參數設定, 設定變數及初始值
"""
size = 1                 # 帶電球體半徑
d = 10                   # 帶電球體連心線距離
L = 1.8*d                # 畫面寬度
q1, c1 = 5, color.blue   # 帶電球體1電量、顏色
q2, c2 = -5, color.red   # 帶電球體2電量、顏色
ke = 8.989E9             # 靜電力常數
N = 6                    # 將顯示的空間每邊切成 N 等份

"""
 2. 產生帶電球體類別, 回傳帶電球體產生的電場
"""
class Ball:
    def __init__(self, pos, radius, color, charge):
        self.pos = pos
        self.radius = radius
        self.color = color
        self.charge = charge
        self.ball = sphere(pos=vec(self.pos), charge=self.charge, color=self.color)
    def electric(self, pos2):
        return ke*self.charge / mag2(pos2-self.pos) * norm(pos2-self.pos)

"""
 3. 畫面設定
"""

# 產生動畫視窗
scene = canvas(title="Electric Field", width=600, height=600, x=0, y=0, background=color.black, range=L)

# 產生帶電球體1、2
b1 = Ball(pos=vec(-d/2, 0, 0), radius=size, color=c1, charge=q1)
b2 = Ball(pos=vec(d/2, 0, 0), radius=size, color=c2, charge=q2)

# 計算畫箭頭的位置, 如果不在帶電球體內則加到串列 locations 當中
locations = []
for i in range(N+1):
    for j in range(N+1):
        for k in range(N+1):
            location = vec(L/N*i - L/2, L/N*j - L/2, L/N*k - L/2)
            if(mag(location-b1.pos) > 2*size and mag(location-b2.pos) > 2*size):
                locations.append(location)

# 依序讀取串列 locations 的元素, 在對應的位置產生箭頭
fields = []
for location in locations:
    fields.append(arrow(pos=location, axis=vec(0, 0, 0), color=color.green))

# 更新箭頭的長度及方向, 長度乘以 1E6 才能看見, 記錄重力場強度最大值, 量值接近最大值偏紅色, 量值接近 0 偏綠色
fmax = 0
for field in fields:
    field.axis = (b1.electric(field.pos) + b2.electric(field.pos))*1E-9
    if(field.axis.mag >= fmax): fmax = field.axis.mag

for field in fields:
    field.color = vec(field.axis.mag/fmax, 1 - field.axis.mag/fmax, 0)

程式 18-4 與 18-3 幾乎一模一樣，以下只解釋不同之處：

將球體半徑、畫面寬度大幅縮小，設定帶電球體的電量、靜電力常數。
第 23 ~ 31 行，自訂類別 Ball 時不採用繼承的方式，產生此類別的物件時要輸入位置 pos、半徑 radius、顏色 color、電量 charge，產生一個球體，並且自訂計算電場用的函式 electric。
若左側小球電量為 q1、右側小球電量為 q2，以下是 3 種不同組合的模擬結果。