<style>
.markdown-body table{
display: unset;
}
</style>

# 圓周運動

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

<br />

如果一個小球的質量為 $m$，速度為 $v$，角速度為 $\omega$，轉彎時的曲率半徑為 $R$，則轉彎需要的向心加速度為
$$
a_c = \frac{v^2}{R} = R \omega^2
$$

由牛頓第二運動定律可得向心力為

$$
F_c = ma_c = m \cdot \frac{v^2}{R} = m \cdot R \omega^2
$$

如果想要用 VPython 畫出小球在水平面上做等速率圓周運動，我們需要想辦法計算向心加速度的大小與方向。如果成功地畫出水平面上的等速率圓周運動，也許就可以進一步挑戰鉛直面圓周運動。以下共有3個程式：

1. 只畫出小球的等速率圓周運動。　（[GlowScript 網站動畫連結](http://www.glowscript.org/#/user/yizhe/folder/Public/program/7-1circle)）
2. 畫出小球、繩子、轉軸的等速率圓周運動。　（[GlowScript 網站動畫連結](http://www.glowscript.org/#/user/yizhe/folder/Public/program/7-2circlerope)）
3. 鉛直面圓周運動，畫出速率、切線加速度、法線加速度與時間的關係圖，計算週期。　（[GlowScript 網站動畫連結](http://www.glowscript.org/#/user/yizhe/folder/Public/program/7-3verticalcircle)）

<br />

## 程式 7-1：圓周運動 （[取得程式碼](https://github.com/YiZheWangTw/VPythonTutorial/blob/master/07.%E5%9C%93%E5%91%A8%E9%81%8B%E5%8B%95/7-1_circle.py)）

"""
 VPython教學: 7-1.圓周運動
 Ver. 1: 2018/2/22
 Ver. 2: 2019/9/7
 作者: 王一哲
"""
from vpython import *

"""
 1. 參數設定, 設定變數及初始值
"""
size = 0.5            # 小球半徑
v0 = 10               # 小球初速
R = 5                 # 圓周運動半徑
L = 4*R               # 地板長度
t = 0                 # 時間
dt = 0.001            # 時間間隔

"""
 2. 畫面設定
"""
scene = canvas(title="Circle", width=800, height=400, x=0, y=0, background=vec(0, 0.6, 0.6))
scene.camera.pos = vec(0, L/2, L/2)
scene.camera.axis = vec(0, -L/2, -L/2)
floor = box(pos=vec(0, 0, 0), size=vec(L, 0.01, L), texture=textures.metal)
ball = sphere(pos=vec(R, size, 0), radius=size, color=color.red, make_trail=True, retain=100, v=vec(0, 0, -v0))
arrow_v = arrow(pos=ball.pos, axis=ball.v, radius=0.2*size, shaftwidth=0.4*size, color=color.green)
arrow_a = arrow(pos=ball.pos, axis=vec(0, 0, 0), radius=0.2*size, shaftwidth=0.4*size, color=color.blue)

"""
 3. 物體運動部分
"""
while True:
    rate(1000)
    axis = ball.pos - vec(0, size, 0)
    ball.a = -(ball.v.mag2/R)*axis.norm()
    ball.v += ball.a*dt
    ball.pos += ball.v*dt
    arrow_v.pos = ball.pos
    arrow_v.axis = ball.v
    arrow_a.pos = ball.pos
    arrow_a.axis = ball.a
    t += dt

<br />

<img style="display: block; margin-left: auto; margin-right: auto" height="100%" width="100%" src="https://i.imgur.com/U2PaMd4.png">
<div style="text-align:center">水平面圓周運動</div>
<br />

### 參數設定

在此定義的變數有 size、v0、R、L、t、dt，用途都已經寫在該行的註解中。
<br />

### 畫面設定

1. 這裡用到了一個新的功能

    scene.camera.pos = vec(0, L/2, L/2)
    scene.camera.axis = vec(0, -L/2, -L/2)

camera.pos 用來觀察者所在的位置，camera.axis 則是設定觀察者看畫面的方向。

2. 小球是在 xz 平面上運動，出發點在畫面正右方距離 R 處，初速度方向朝 -z 軸、量值為 v0。

3. arrow_v 和 arrow_a 是用來表示小球速度和加速度的箭頭。
<br />

### 物體運動

1. 用 **axis = ball.pos - vec(0, size, 0)** 找出小球相對於轉軸的位置向量。

2. 用 **ball.a = -(ball.v.mag2 / R) * axis.norm()** 計算小球的向心加速度。其中 mag2 是用來計算向量的量值平方，假設向量的名稱為 A ，語法有以下兩種

    A.mag2 = mag2(A)

其中 **norm()** 是用來計算單位向量，假設向量的名稱為 A ，語法有以下兩種

    A.norm() = norm(A)

3. 最後更新小球的速度、位置，更新箭頭的起點位置、方向及長度，更新時間。
<br />

### 2023年4月22日補充：attach_arrow

VPython 提供了 attach_arrow 功能 [[5](https://www.glowscript.org/docs/VPythonDocs/trail.html)]，可以將箭頭加在某個物件上，箭頭位置會隨著物件移動，長度及方向也會隨著物件性質而改變。如果使用 attach_arrow 功能，可以將程式碼第26 ~ 28行改成

ball = sphere(pos=vec(R, size, 0), radius=size, color=color.red, make_trail=True, retain=100, v=vec(0, 0, -v0), a=vec(-v0**2/R, 0, 0))
arrow_v = attach_arrow(ball, "v", radius=0.2*size, shaftwidth=0.4*size, color=color.green)
arrow_a = attach_arrow(ball, "a", radius=0.2*size, shaftwidth=0.4*size, color=color.blue)

同時刪除程式碼第39 ~ 42行，也可以做到相同的動畫效果，但是不需要在 while 迴圈中更新箭頭的位置、長度、方向，程式碼會更簡潔。
<br /><br />

## 程式 7-2.圓周運動, 加上轉軸及繩子 （[取得程式碼](https://github.com/YiZheWangTw/VPythonTutorial/blob/master/07.%E5%9C%93%E5%91%A8%E9%81%8B%E5%8B%95/7-2_circle_rope.py)）

"""
 VPython教學: 7-2.圓周運動, 加上轉軸及繩子
 Ver. 1: 2018/2/22
 Ver. 2: 2019/9/7
 作者: 王一哲
"""
from vpython import *

"""
 1. 參數設定, 設定變數及初始值
"""
size = 0.5            # 小球半徑
v0 = 10               # 小球初速
R = 5                 # 圓周運動半徑
L = 4*R               # 地板長度
t = 0                 # 時間
dt = 0.001            # 時間間隔

"""
 2. 畫面設定
"""
scene = canvas(title="Circle with Rope", width=800, height=400, x=0, y=0, background=vec(0, 0.6, 0.6))
scene.camera.pos = vec(0, L/2, L/2)
scene.camera.axis = vec(0, -L/2, -L/2)
floor = box(pos=vec(0, -size, 0), size=vec(L, 0.01, L), texture=textures.metal)
ball = sphere(pos=vec(R, 0, 0), radius=size, color=color.red, make_trail=True, retain=100, v=vec(0, 0, -v0))
center = cylinder(pos=vec(0, -size, 0), axis=vec(0, 2*size, 0), radius=0.1*size, color=color.white)
rope = cylinder(pos=vec(0, 0, 0), axis=ball.pos, radius=0.1*size, color=color.yellow)
arrow_v = arrow(pos=ball.pos, axis=ball.v, radius=0.2*size, shaftwidth=0.4*size, color=color.green)
arrow_a = arrow(pos=ball.pos, axis=vec(0, 0, 0), radius=0.2*size, shaftwidth=0.4*size, color=color.blue)

"""
 3. 物體運動部分
"""
while True:
    rate(1000)
    axis = ball.pos
    ball.a = -(ball.v.mag2/R)*axis.norm()
    ball.v += ball.a*dt
    ball.pos += ball.v*dt
    rope.axis = axis
    arrow_v.pos = ball.pos
    arrow_v.axis = ball.v
    arrow_a.pos = ball.pos
    arrow_a.axis = ball.a
    t += dt

<br />

<img style="display: block; margin-left: auto; margin-right: auto" height="100%" width="100%" src="https://i.imgur.com/dFbQ3tX.png">
<div style="text-align:center">水平圓周運動（加上繩子和轉軸）</div>
<br />

<img style="display: block; margin-left: auto; margin-right: auto" height="100%" width="100%" src="https://i.imgur.com/ZxVdzHQ.png">
<div style="text-align:center">水平圓周運動（繩子被表示加速度的箭頭擋住）</div>
<br />

程式 7-2 和 7-1 非常像，只是增加了繩子和轉軸，新增的程式碼為27、28、41行

center = cylinder(pos = vec(0, -size, 0), axis = vec(0, 2*size, 0), radius = 0.1*size, color = color.white)
rope = cylinder(pos = vec(0, 0, 0), axis = ball.pos, radius = 0.1*size, color = color.yellow)
rope.axis = axis

分別用來畫出轉軸、畫出繩子、更新繩子的長度及方向。

如果想要隱藏顯示速度、加速度用的箭頭，可以註解第 29、30、42 ~ 45 行，這樣就會只畫出球、軸心、繩子。
<br /><br />

VPython 當中有另一個畫曲線的物件，名稱為 curve，如果只將兩個點連線，看起來就會是直線，可以用來畫軸心及繩子，第27、28、41行程式碼分別改為

center = curve(pos=[vec(0, -size, 0), vec(0, size, 0)], radius=0.1*size, color=color.white)
rope = curve(pos=[ball.pos, vec(0, 0, 0)], radius=0.1*size, color=color.yellow)
rope.modify(0, pos=ball.pos)

<br /><br />

## 程式7-3.鉛直面圓周運動 （[取得程式碼](https://github.com/YiZheWangTw/VPythonTutorial/blob/master/07.%E5%9C%93%E5%91%A8%E9%81%8B%E5%8B%95/7-3_vertical_circle.py)）

from vpython import *

"""
 1. 參數設定, 設定變數及初始值
"""
size = 0.5            # 小球半徑
R = 5                 # 圓周運動半徑
g = 9.8               # 重力加速度 9.8 m/s^2
v0 = 1*sqrt(g*R)      # 小球初速, 1 ~ 7 sqrt(g*R)
ratio = 0.1           # 速度, 加速度箭頭長度與實際的比例
i = 0                 # 小球回到出發點的次數
t = 0                 # 時間
dt = 0.0001           # 時間間隔, 取 0.0001 以降低誤差

"""
 2. 畫面設定
"""
scene = canvas(title="Vertical Circle", width=600, height=600, x=0, y=0, background=color.black)
ball = sphere(pos=vec(0, R, 0), radius=size, color=color.red, make_trail=True, retain=200, v=vec(-v0, 0, 0))
center = cylinder(pos=vec(0, 0, -size), axis=vec(0, 0, 2*size), radius=0.1*size, color=color.white)
rope = cylinder(pos=vec(0, 0, 0), axis=ball.pos, radius=0.1*size, color=color.yellow)
arrow_v = arrow(pos=ball.pos, axis=vec(0, 0, 0), radius=0.2*size, shaftwidth=0.4*size, color=color.green)
arrow_a = arrow(pos=ball.pos, axis=vec(0, 0, 0), radius=0.2*size, shaftwidth=0.4*size, color=color.blue)
gd = graph(title="plot", width=600, height=450, x=0, y=600, xtitle="<i>t</i>(s)",
           ytitle="green: <i>v</i>(m/s), red: <i>a<sub>t</sub></i>(m/s<sup>2</sup>), blue: <i>a<sub>n</sub></i>(m/s<sup>2</sup>)")
vt_plot = gcurve(graph=gd, color=color.green)
at_plot = gcurve(graph=gd, color=color.red)
an_plot = gcurve(graph=gd, color=color.blue)

"""
 3. 自訂函式, findan 計算法線加速度, findat 計算切線加速度
"""
def findan(v, pos):
    an = -v.mag2 / R * pos.norm()
    return an

def findat(pos):
    x = pos.x
    y = pos.y
    r = sqrt(x**2 + y**2)
    sintheta = abs(x)/r
    costheta = abs(y)/r
    absat = g*sintheta
    aty = -absat*sintheta
    if (x <= 0 and y <= 0) or (x >=0 and y>= 0):
        atx = +absat*costheta
    elif (x <= 0 and y >= 0) or (x >= 0 and y <= 0):
        atx = -absat*costheta
    at = vec(atx, aty, 0)
    return at

"""
 4. 物體運動部分, 小球回到出發點 5 次停止運作
"""
while i < 5:
# 由於 dt 較小，每秒計算 5000 次使動畫速度加快
    rate(5000)
# xp 是小球原來的位置, xc 是小球現在的位置, 用來判斷小球是否回到出發點    
# 計算小球 an, at, 更新加速度, 速度, 位置
    xp = ball.pos.x
    an = findan(ball.v, ball.pos)
    at = findat(ball.pos)
    ball.a = an + at
    ball.v += ball.a*dt
    ball.pos += ball.v*dt
    xc = ball.pos.x
    rope.axis = ball.pos
# 若小球回到出發點, 將 i 加 1, 印出時間 t, 由於誤差會累積, 取第一次回到出發點的時間作為週期
    if xp > 0 and xc < 0:
        i += 1
        print(i, t)
# 更新代表速度, 加速度的箭頭
    arrow_v.pos = ball.pos
    arrow_v.axis = ball.v * ratio
    arrow_a.pos = ball.pos
    arrow_a.axis = ball.a * ratio
# 更新 v-t, at-t, an-t 圖
    vt_plot.plot(t, ball.v.mag)
    at_plot.plot(t, at.mag)
    an_plot.plot(t, an.mag)
# 更新時間
    t += dt

<br />

<img style="display: block; margin-left: auto; margin-right: auto" height="80%" width="80%" src="https://i.imgur.com/o15CTqK.png">
<div style="text-align:center">鉛直圓周運動</div>
<br />

<img style="display: block; margin-left: auto; margin-right: auto" height="80%" width="80%" src="https://imgur.com/fPVSVaB.png">
<div style="text-align:center">鉛直圓周運動速度、切線加速度、法線加速度與時間關係圖</div>
<br />

### 理論分析

為了使計算過程較為簡單，我將圓心設定在原點，小球由 \\( (0, R, 0) \\) 處以初速度 \\( (-v_0, 0, 0) \\) 出發。加速度可以分為改變速度方向用的法線加速度 \\( a_n \\)，也就是向心加速度 \\( a_c \\)，以及改變速度量值用的切線加速度 \\( a_t \\)。\\( a_n \\) 的計算方式與程式 7-1 相同，不再贅述。$a_t$ 的來源為重力加速度的切線方向分量，假設小球與鉛垂線的夾角為$\theta$，則

$$ a_t = g \sin \theta $$

$a_t$ 的 y 方向分量方向下向，量值為

$$ a_{t,y} = a_t \sin \theta = g \sin^2 \theta $$

$a_t$ 的 $x$ 方向分量方向與位置有關，在第1和第3象限向右，在第2和第4象限向左，量值為

$$ a_{t,x} = a_t \cos \theta = g \sin \theta \cos \theta $$

如果只想找出週期，可以改由力學能守恆計算。假設小球在最高點的速度量值為

$$ v_0 = n \sqrt{gR} $$

再由力學能守恆可寫出任意點與最高點的關係式 \[9\]

$$ \frac{1}{2}mv^2 + mgR \cos \theta = \frac{1}{2}mn^2 gR + mgR $$

$$ v = R \frac{d \theta}{dt} = \sqrt{(n^2 + 2 - 2\cos \theta)gR} $$

$$ T = \int_0^T dt = 2 \sqrt{\frac{R}{g}} \int_0^{\pi} \frac{d \theta}{\sqrt{n^2 + 2 - 2\cos \theta}} $$

<br />

### 參數設定

在此定義的變數有 size、R、g、v0、ratio、i、t、dt，用途都已經寫在該行的註解中。為了減少代入的時間長度造成的誤差，將 dt 的值調整為 0.0001。
<br />

### 畫面設定

1. 小球是在 xy 平面上運動，出發點在畫面正上方距離 R 處，初速度方向朝 -x 軸、量值為 v0，不需要調整觀察者的位置及方向。
2. 轉軸的方向改為指向 +z 軸方向。
3. arrow_v 和 arrow_a 是用來表示小球速度和加速度的箭頭。
4. 用 graph 開啟繪圖視窗，用 gcurve 畫出 v - t、at - t、an - t 關係圖。

<br />

### 自訂函式

1. 自訂函式 findan，輸入的參數為 v 和 pos，格式皆為向量，用來計算小球的法線加速度an。
2. 自訂函式 findat，輸入的參數為 pos，格式為向量，用來計算小球的切線加速度 \\( a_t \\)。繩子和鉛垂線的夾角為 \\( \theta \\)，先計算 \\( \sin \theta \\) 和 \\( \cos \theta \\) ，由於 \\( \theta \\)只 取銳角，因此兩者皆為正值。再計算 \\( a_{t,y} \\) 和 \\( a_{t,x} \\) ，並由 pos 判斷 \\( a_{t,x} \\) 的正、負號。

<br />

### 物體運動

1. 由小球的位置判斷是否回到出發點。由於小球是由最高點開始逆時鐘方向旋轉，小球回到出發點的條件為：小球原來的位置 xp 在 +x 區域、現在的位置 xc 在 -x 區域。如果回到出發點，印出經過的時間 t、將次數 i 加 1，當小球回到出發點 5 次時停止動畫。
2. 由於 dt 較小，將每秒計算的次數改為 5000，使動畫速度加快。
3. 使用自訂函式 findan、findat 計算小球的法線加速度及切線加速度。
4. 更新小球的速度、位置，更新箭頭的起點位置、方向及長度，更新 v - t、at - t、an - t 關係圖，更新時間。

<br />

### 數據處理部分

若小球在最高點的速度量值為 \\( v_0 = n \sqrt{gR} \\)，將 n 用 1 ~ 7 代入，分別由模擬程式及理論計算找出週期，數據如下

<div style="text-align:center">

| n | v0 | T理論值 | T模擬值 |
| - | -- | :----------------: | :----------------: |
| 1 | 7 | 2.8841511713977500 | 2.8838000000016613 |
| 2 | 14 | 1.8728982530658000 | 1.8728999999998102 |
| 3 | 21 | 1.3617223964215400 | 1.3617999999999999 |
| 4 | 28 | 1.0602946238731300 | 1.0602999999998997 |
| 5 | 35 | 0.8646052540953560 | 0.8645999999999211 |
| 6 | 42 | 0.7284267507699600 | 0.7283999999999361 |
| 7 | 49 | 0.6286255969937600 | 0.6285999999999471 |

</div>

<br />
理論值與模擬值幾乎一樣。將數據匯入 SciDAVis 畫出 T - v0 關係圖如下

<img style="display: block; margin-left: auto; margin-right: auto" height="100%" width="100%" src="https://i.imgur.com/E3ihBGt.jpg">
<div style="text-align:center"><i>T</i> - <i>v</i>₀ 關係圖</div>
<br /><br />

如果想要用 Python 計算數值積分，需要用到另一個套件 [SymPy](https://www.sympy.org/en/index.html)，程式碼如下：

from __future__ import division
from sympy import *

x = symbols('x', communtative = True)

def f(n, x):
    return 1/sqrt(n**2 + 2 - 2*cos(x))

for n in range(1, 8, 1):
    print("n = {:d} \tT = {:.16f}".format(n, integrate(f(n, x), (x, 0, pi)).evalf() * 2 * sqrt(5/9.8)))

<br />

## 結語

模擬鉛直圓周運動時，如果代入的 dt 較長，誤差就會很快累積並放大，接下來的動畫及數據就會與原來預設的條件相差很多，這是需要特別注意的地方。
<br />

## 參考資料
1. **canvas**: http://www.glowscript.org/docs/VPythonDocs/canvas.html
2. **box**: http://www.glowscript.org/docs/VPythonDocs/box.html
3. **sphere**: http://www.glowscript.org/docs/VPythonDocs/sphere.html
4. **arrow**: https://www.glowscript.org/docs/VPythonDocs/arrow.html
5. **attach_arrow**: https://www.glowscript.org/docs/VPythonDocs/trail.html
6. **vec**: https://www.glowscript.org/docs/VPythonDocs/vector.html
7. **cylinder**: https://www.glowscript.org/docs/VPythonDocs/cylinder.html
8. **curve**: https://www.glowscript.org/docs/VPythonDocs/curve.html
9. 徐國誠(2013)。鉛直面圓周運動的週期。**翰林自然科學天地**，**42**，10-15。

---

###### tags:`VPython`