# 重力場及電場

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

<br />

若在空間中有一個質量為 $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|} $$

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

<img style="display: block; margin-left: auto; margin-right: auto" height="80%" width="80%" src="https://imgur.com/sCJld1d.png">
<div style="text-align:center">地球重力場示意圖</div>
<br />

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

<br />

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

$$ \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 就能畫出兩個帶電球體在空間中建立的電場。

<br />

## 程式 18-1. 用 for 迴圈產生箭頭 （[取得程式碼](https://github.com/YiZheWangTw/VPythonTutorial/blob/master/18.%E9%87%8D%E5%8A%9B%E5%A0%B4%E5%8F%8A%E9%9B%BB%E5%A0%B4/18-1_arrows.py)）　（[GlowScript 網站動畫連結](http://www.glowscript.org/#/user/yizhe/folder/Public/program/18-1arrows)）

"""
 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=vec(4, 0, 0), radius=1, color=color.green))

<br />

這個程式與之前的程式相比十分簡短，只要將箭頭畫完就可以了，因此連產生動畫的部分都沒有。比較特別之處在於使用了 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。程式執行的成果如下圖。
<br />

<img style="display: block; margin-left: auto; margin-right: auto" height="80%" width="80%" src="https://i.imgur.com/1NnZkKr.png">
<div style="text-align:center">程式18-1畫面截圖</div>
<br />

## 程式 18-2. 重力場, 可以改成畫不同的星球 （[取得程式碼](https://github.com/YiZheWangTw/VPythonTutorial/blob/master/18.%E9%87%8D%E5%8A%9B%E5%A0%B4%E5%8F%8A%E9%9B%BB%E5%A0%B4/18-2_gravitational_field.py)）　（[GlowScript 網站動畫連結](http://www.glowscript.org/#/user/yizhe/folder/Public/program/18-2gravitationalfield)）

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

<br />

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

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

<img style="display: block; margin-left: auto; margin-right: auto" height="80%" width="80%" src="https://imgur.com/sCJld1d.png">
<div style="text-align:center">程式18-2畫面截圖</div>
<br />

## 程式 18-3. 重力場, 可以改成畫不同的星球 （[取得程式碼](https://github.com/YiZheWangTw/VPythonTutorial/blob/master/18.%E9%87%8D%E5%8A%9B%E5%A0%B4%E5%8F%8A%E9%9B%BB%E5%A0%B4/18-3_gravitational_field_2.py)）　（[GlowScript 網站動畫連結](http://www.glowscript.org/#/user/yizhe/folder/Public/program/18-3gravitationalfield2)）

"""
 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(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(pos=vec(-d/2, 0, 0), radius=radius["Earth"], m=mass["Earth"], texture=textures.earth)
mars = Planet(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)

<br />

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

<img style="display: block; margin-left: auto; margin-right: auto" height="80%" width="80%" src="https://imgur.com/ztcJizC.png">
<div style="text-align:center">程式18-3畫面截圖</div>
<br />

## 程式 18-4. 電場, 2個球狀帶電體 （[取得程式碼](https://github.com/YiZheWangTw/VPythonTutorial/blob/master/18.%E9%87%8D%E5%8A%9B%E5%A0%B4%E5%8F%8A%E9%9B%BB%E5%A0%B4/18-4_electric_field.py)）　（[GlowScript 網站動畫連結](http://www.glowscript.org/#/user/yizhe/folder/Public/program/18-4electricfield)）

"""
 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.988E9             # 靜電力常數
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), radius=self.radius, 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))

# 更新箭頭的長度及方向, 長度乘以 1E-9 才不會太長, 記錄電場強度最大值, 量值接近最大值偏紅色, 量值接近 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)

<br />

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

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

<img style="display: block; margin-left: auto; margin-right: auto" height="80%" width="80%" src="https://imgur.com/GFQ8aD7.png">
<div style="text-align:center">q1 = 1, q2 = -5 模擬結果畫面截圖</div>
<br />

<img style="display: block; margin-left: auto; margin-right: auto" height="80%" width="80%" src="https://imgur.com/0d5oPc1.png">
<div style="text-align:center">q1 = 5, q2 = -5 模擬結果畫面截圖</div>
<br />

<img style="display: block; margin-left: auto; margin-right: auto" height="80%" width="80%" src="https://imgur.com/tvDP67Q.png">
<div style="text-align:center">q1 = 5, q2 = -1 模擬結果畫面截圖</div>
<br />

## 結語

這是我們第一次嘗試用 VPython 畫出向量場的示意圖，看起來效果還不錯，之後應該可以用類似的方法畫出載流導線或線圈產生的磁場示意圖，請參考〈[電流的磁效應](https://hackmd.io/@yizhewang/r1g8P0uf7)〉。

<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. **arrow**: http://www.glowscript.org/docs/VPythonDocs/arrow.html

<br />

---

## 補充 for 迴圈的使用方法

> 2020/6/10：3 層 for 迴圈
> 2021/12/23：for 迴圈基本語法

### 1 層 for 迴圈

若程式碼為

for i in range(5): print(i)

執行時會輸出

0
1
2
3
4

<br />

range是用來產生數列的函式，格式為

range(起始值, 結束值, 增量)

其中起始值預設為0，增量預設為1，執行時不包含結束值。
<br />

若程式碼為

for i in range(1, 10, 2): print(i)

執行時會輸出

1
3
5
7
9

<br />

若程式碼為

for i in range(10, 1, -2): print(i)

執行時會輸出

10
8
6
4
2

<br />

### 2 層 for 迴圈

若程式碼為

for i in range(4):
    for j in range(4):
        print("i = {:d}     j = {:d}".format(i, j))

執行時會輸出

i = 0     j = 0
i = 0     j = 1
i = 0     j = 2
i = 0     j = 3
i = 1     j = 0
i = 1     j = 1
i = 1     j = 2
i = 1     j = 3
i = 2     j = 0
i = 2     j = 1
i = 2     j = 2
i = 2     j = 3
i = 3     j = 0
i = 3     j = 1
i = 3     j = 2
i = 3     j = 3

<br />

如果想要印出九九乘法表可以這樣寫

for i in range(2, 10):
    for j in range(1, 10):
        print("{:d} * {:d} = {:d}".format(i, j, i*j), end="\t")
    print()

執行時會輸出

2 * 1 = 2	2 * 2 = 4	2 * 3 = 6	2 * 4 = 8	2 * 5 = 10	2 * 6 = 12	2 * 7 = 14	2 * 8 = 16	2 * 9 = 18	
3 * 1 = 3	3 * 2 = 6	3 * 3 = 9	3 * 4 = 12	3 * 5 = 15	3 * 6 = 18	3 * 7 = 21	3 * 8 = 24	3 * 9 = 27	
4 * 1 = 4	4 * 2 = 8	4 * 3 = 12	4 * 4 = 16	4 * 5 = 20	4 * 6 = 24	4 * 7 = 28	4 * 8 = 32	4 * 9 = 36	
5 * 1 = 5	5 * 2 = 10	5 * 3 = 15	5 * 4 = 20	5 * 5 = 25	5 * 6 = 30	5 * 7 = 35	5 * 8 = 40	5 * 9 = 45	
6 * 1 = 6	6 * 2 = 12	6 * 3 = 18	6 * 4 = 24	6 * 5 = 30	6 * 6 = 36	6 * 7 = 42	6 * 8 = 48	6 * 9 = 54	
7 * 1 = 7	7 * 2 = 14	7 * 3 = 21	7 * 4 = 28	7 * 5 = 35	7 * 6 = 42	7 * 7 = 49	7 * 8 = 56	7 * 9 = 63	
8 * 1 = 8	8 * 2 = 16	8 * 3 = 24	8 * 4 = 32	8 * 5 = 40	8 * 6 = 48	8 * 7 = 56	8 * 8 = 64	8 * 9 = 72	
9 * 1 = 9	9 * 2 = 18	9 * 3 = 27	9 * 4 = 36	9 * 5 = 45	9 * 6 = 54	9 * 7 = 63	9 * 8 = 72	9 * 9 = 81

<br />

### 3 層 for 迴圈

若程式碼為

N, counts = 4, 0

for i in range(N):
    for j in range(N):
        for k in range(N):
            counts += 1
            print("i = {:d}, j = {:d}, k = {:d}".format(i, j, k))

print("counts = {:d}".format(counts))

執行時會輸出

i = 0, j = 0, k = 0
i = 0, j = 0, k = 1
i = 0, j = 0, k = 2
i = 0, j = 0, k = 3
i = 0, j = 1, k = 0
i = 0, j = 1, k = 1
i = 0, j = 1, k = 2
i = 0, j = 1, k = 3
i = 0, j = 2, k = 0
i = 0, j = 2, k = 1
i = 0, j = 2, k = 2
i = 0, j = 2, k = 3
i = 0, j = 3, k = 0
i = 0, j = 3, k = 1
i = 0, j = 3, k = 2
i = 0, j = 3, k = 3
i = 1, j = 0, k = 0
i = 1, j = 0, k = 1
i = 1, j = 0, k = 2
i = 1, j = 0, k = 3
i = 1, j = 1, k = 0
i = 1, j = 1, k = 1
i = 1, j = 1, k = 2
i = 1, j = 1, k = 3
i = 1, j = 2, k = 0
i = 1, j = 2, k = 1
i = 1, j = 2, k = 2
i = 1, j = 2, k = 3
i = 1, j = 3, k = 0
i = 1, j = 3, k = 1
i = 1, j = 3, k = 2
i = 1, j = 3, k = 3
i = 2, j = 0, k = 0
i = 2, j = 0, k = 1
i = 2, j = 0, k = 2
i = 2, j = 0, k = 3
i = 2, j = 1, k = 0
i = 2, j = 1, k = 1
i = 2, j = 1, k = 2
i = 2, j = 1, k = 3
i = 2, j = 2, k = 0
i = 2, j = 2, k = 1
i = 2, j = 2, k = 2
i = 2, j = 2, k = 3
i = 2, j = 3, k = 0
i = 2, j = 3, k = 1
i = 2, j = 3, k = 2
i = 2, j = 3, k = 3
i = 3, j = 0, k = 0
i = 3, j = 0, k = 1
i = 3, j = 0, k = 2
i = 3, j = 0, k = 3
i = 3, j = 1, k = 0
i = 3, j = 1, k = 1
i = 3, j = 1, k = 2
i = 3, j = 1, k = 3
i = 3, j = 2, k = 0
i = 3, j = 2, k = 1
i = 3, j = 2, k = 2
i = 3, j = 2, k = 3
i = 3, j = 3, k = 0
i = 3, j = 3, k = 1
i = 3, j = 3, k = 2
i = 3, j = 3, k = 3
counts = 64

<br />

---

###### tags:`VPython`