# Lecture 2 Vpython動力學模擬 - F=ma的數值解
>建國高中特色選修課程 - 物理現象的程式設計與模擬
>作者：曾靖夫、賴奕帆
>日期：2018/7/5
---
## 一、牛頓第二運動定律與迴圈
在牛頓力學中，牛頓第二運動定律是我們處理所有運動都會使用的定律。原因是你只要知道初始條件，當受力決定了速度就決定了（$\Delta \vec{v} = \vec{a}\ \Delta t$），而速度決定了位置就決定了（$\Delta \vec{r} = \vec{v}\ \Delta t$）。因此你只要知道物體受了什麼力，就可以預測未來每一瞬間的運動狀態。你也可以說這是力學問題的暴力解法，因為即使作用力再複雜都一定可以算出正確的物體運動狀態，只是不是用手算而是讓電腦幫我們算，即是在科學計算所說的「數值解」。
這裡簡單介紹一下牛頓第二運動定律的基本內容$$\vec{F}_{合} = m\vec{a}$$式中的 $\vec{F}_{合}$ 為作用在物體的所有外力合，$m$ 為所選定討論的受力物，$\vec{a}$ 為物體的加速度。這裡只要 $\vec{F}_{合}$ 確定了，加速度就決定了，加速度就會造成速度變化，有了速度就會造成位移。所以我們就可以在程式中將作用力寫入，電腦就會自動模擬出之後每一刻的運動狀態。
在這個基礎下，使用迴圈去模擬每個時間間隔所發生的變化是最適當的，只要設定的時間間隔夠短，模擬出來的結果就近似為每一時刻的瞬時狀態。以下簡單分析每個迴圈產生的變化：
<style>
.markdown-body table{
display: unset;
}
</style>
<center>
| loop | 時間 t = t + dt | 速度 v = v + a*dt | 位置 x = x + v*dt |
|:----:|:---------------:|:-----------------:|:-----------------:|
| 1 | t1 = 0 + dt | v1 = v0 + a0*dt | x1 = x0 + v1*dt |
| 2 | t2 = t1 + dt | v2 = v1 + a1*dt | x2 = x1 + v2*dt |
| 3 | t3 = t2 + dt | v3 = v2 + a2*dt | x3 = x2 + v3*dt |
| ... | ... | ... | ... |
>
</center>
每個迴圈所計算的東西，會累加到下一個迴圈中，產生一連串的連續動態的變化，這就是我們進行動力學模擬的重要基礎。
---
## 二、自由落體運動的模擬（均勻重力場）
這裡我們馬上來分析最熟悉的變速度運動－等加速度運動，這裡以均勻重力場中不受阻力的拋體運動為例，物體只受到固定的重力作用，因此加速度也必為定值。首先我們分析初速度為零的自由落體運動，參考程式如下：
### Exapmle 1 : 等加速度運動－自由落體
```python=
from vpython import *
g = 9.8 #重力加速度 9.8 m/s^2
size = 0.5 #球半徑 0.5 m
height = 15.0 #球初始高度 15 m
m = 1.0 #球質量1kg
Fg = vector(0, -m*g, 0) #重力
scene = canvas(width=600, height=600,x=0, y=0, center = vector(0,height/2,0)) #設定畫面
floor = box(length=20, height=0.01, width=10, color=color.green) #畫地板
ball = sphere(radius = size, color=color.yellow, make_trail= True, trail_type="points", interval=100) #畫球
ball.pos = vector(0, height, 0) #球初始位置
ball.v = vector(0, 0, 0) #球初速
dt = 0.001 #時間間隔 0.001 秒
t = 0.0 #模擬初始時間為0秒
while ball.pos.y > size: #模擬直到球落地 即y=球半徑
rate(1/dt) #每一秒跑 1000 次
t = t + dt #計時器
ball.a = Fg/m #球的加速度
ball.v = ball.v + ball.a*dt #球的末速度 = 前一刻速度 + 加速度*時間間隔
ball.pos = ball.pos + ball.v * dt #球的末位置 = 前一刻位置 + 速度*時間間隔
print (t, ball.v)
```
執行後物件視窗的畫面如下：
<img style="display: block; margin-left: auto; margin-right: auto" height="60%" width="60%"
src="https://i.imgur.com/4VKvSuz.png">
:::info
* **物件的軌跡控制**
1. make_trail=True：這是Vpython畫出物件運動軌跡的指令，預設的軌跡為一條連續的細線，若設定為False則不會顯現出軌跡。
2. trail_type="points"：這是將軌跡線設定為「點」的形式。
3. interval=100：這是設定每多少單位打一個「點」，這裡的100指的是每100個迴圈打一個點。在這個例子中100個迴圈代表時間經過100*0.001秒，即每0.1秒打一個點。同學可將其理解為立體的打點計時器，就可以透過觀察點和點的間隔有沒有變化來判斷物體運動型態的改變。
:::
接下來如果我們要使小球碰到地面之後反彈，則需要引入程式撰寫中很常用到的條件判斷結構－<font color = "red">if和else</font>。這裡我們假設球碰到地面的時候，以相同度速率彈起，如果整個過程沒有受到任何阻力，小球應該會達到與釋放時相同的高度，在程式中我們很容易驗證這個觀念。在程式中只須增加一點點while迴圈的內容，程式範例如下：
```python=
while True: #無窮迴圈
rate(1000)
ball.v.y += -ball.a*dt #此寫法意義上與ball.v = ball.v + ball.a*dt完全相同
ball.pos += ball.v*dt #此寫法意義上與ball.pos = ball.pos + ball.v*dt完全相同
if ball.pos.y <= size and ball.v.y < 0: #條件：球心高度小於球半徑且速度沿-y軸
ball.v.y = - ball.v.y #條件成立則球的速度加一負號表示反彈
```
到目前為止，相信同學已經知道如何寫出在均勻重力場中的等加速度運動，我們只要再加入各種初始條件就可以模擬出各種拋體運動。比如說，你只需要將while迴圈上方的初始速度設定為ball.v = vector(5,5,0)，就會發現球將以仰角 $45^{\circ}$、速率 $5\sqrt{2}$ 斜向拋出。你也可以用三角函數來描述拋射方向，將球度初始速度設定為 $\left ( v_0\ {\rm sin}\theta , v_0\ {\rm cos}\theta, 0 \right )$，將發現球以初速 $v_0$、仰角 $\theta$ 沿 $x-y$ 平面斜向拋出，你也將發現拋體運動是一種很典型的平面運動。
以上我們討論了無空氣阻力的拋體運動，但在現實生活中，這樣的運動太過理想。其實我們也知道空氣阻力對物體運動的影響其實不小，這裡將試著引入空氣阻力。一般而言，在物體運動速度不大的情況下，物體受空氣阻力的大小與其速率成正比，阻力方向與速度方向相反，數學式可以表示為： $$\vec{f}_r=-k\vec{v}$$ 在物體所處環境不變的情況下，式中 $k$ 可視為一常數。
### 檢核作業 2-1
請完成拋體運動的模擬，使球碰到地板會發生反彈，並使用指令描繪出球的運動軌跡。接著請再新增一個球並加入空氣阻力的影響，保留沒有空氣阻力拋體的軌跡與之對照，空氣阻力常數 $k$ 可以自行設定。
<img style="display: block; margin-left: auto; margin-right: auto" height="90%" width="90%"
src="https://i.imgur.com/w5sJci2.png">
執行成果如上圖所示，黃色球為無空氣阻力，紅色為加入空氣阻力的結果。
:::info
* 空氣阻力在程式中的寫法為：AirResistance = - damp * ball.v
damp即是空氣阻力公式中的 $k$。
:::
---
## 三、拋體運動的運動方程式
在這個段落我們還是簡單列出拋體運動的數學方程式，但並不仔細證明推導，目的在於讓同學理解程式工具的重要性。利用程式我們只要用很簡單的觀念將「力」寫入程式中，而不需要經過繁雜的數學推導，就可將物體的運動狀態完整模擬。
如果拋體運動沒有受到空氣阻力，只受到均勻向下的重力作用，我們可以將物體的運動分解為兩個方向，一為沿著水平面的方向，二為沿著鉛直方向。在水平方向上物體並沒有受力，因此在這個方向物體的運動應為等速度運動；而在鉛直方向物體受到固定的重力，所以在這個方向為等加速度運動。
<img style="display: block; margin-left: auto; margin-right: auto" width="70%"
src="https://i.imgur.com/U0dCD9h.png">
如上圖所示，我們分別由上方和側面照光，就可以透過影子很清楚的呈現它在兩個方向上的運動狀態，按這個概念速度速度向量隨時間的關係為：$$\begin{cases} v_x(t)=v_{0x} \\ v_y(t)=v_{0y}-gt \\ v_z(t)=0 \end{cases}$$ 位置向量隨時間的關係為：$$\begin{cases} x(t)=v_{0x}\ t \\ y(t)=v_{0y}\ t - \dfrac{1}{2}gt^2 \\ z(t)=0 \end{cases}$$ 如果考慮了空氣阻力之後，則運動方程式如下所示：$$\begin{cases} x(t)=\dfrac{v_{0x}}{k}\left ( 1-e^{-kt} \right ) \\ y(t)=-\dfrac{gt}{k}+\dfrac{kv_{0y}+g}{k^{2}}\left ( 1-e^{-kt} \right ) \\ z(t)=0 \end{cases}$$這是以正統的力學方法，從牛頓第二運動定律解出來的運動方程式，即「理論解」，光看答案你都可以想像其推導過程的複雜程度。因此用程式數值解最大的優點是，我們只要設定好的初始條件，不必進行冗長的數學推導，電腦就會自動幫我們模擬出物體每瞬間的運動狀態。
### 指定作業 2-1
請利用程式驗證今天在程式中加入空氣阻力的拋體運動模擬與使用牛頓力學解出的運動方程式完全一致。請在視窗中留下三個球與其運動軌跡，如下圖所示，黃色為無阻力的理想拋體，紅色為加入空氣阻力的拋體，綠色為上方方程式所畫出來的軌跡。這邊我們要求執行流程為，黃球和紅球同時拋出，當紅球和黃球都落地就讓它們直接停在地面上，之後再拋出綠球，將發現綠球會貼著紅色軌跡走，確認「數值解」和「理論解」是吻合的。
<img style="display: block; margin-left: auto; margin-right: auto" width="70%"
src="https://i.imgur.com/jVJVCJu.png">
:::success
* 提示：牛頓力學解出的考慮空氣阻力之方程式，在程式中寫法如下
時間 t 以計時器 t = t + dt 累加後代入以下三個方向的式子中運算。
$x$ 方向：$x(t)$
```python
ball.pos.x = ball.v.x * (1 - exp(-k*t))/k
```
$y$ 方向：$y(t)$
```python
ball.pos.y = - g*t/k + (k*ball.v.y + g) * (1 - exp(-k*t))/k**2
```
$z$ 方向：$z(t)$
```python
ball.pos.z = 0
```
:::
---
## 四、三角函數簡介
在物理學中很常利用三角函數來表示空間向量在各個方向上的分量。而在數學上三角函數的定義非常單純，只是直角三角形的邊角關係，以下我們對照直角座標系，給出三種三角函數的定義：
<img style="display: block; margin-left: auto; margin-right: auto" width="35%"
src="https://i.imgur.com/lqhNADC.png">
1. 餘弦函數：$\ {\rm cos} \theta = \dfrac{x}{r}$ （投影到 $x$ 軸）
2. 正弦函數：$\ {\rm sin} \theta = \dfrac{y}{r}$ （投影到 $y$ 軸）
3. 正切函數：$\ {\rm tan} \theta = \dfrac{y}{x}$ （斜率）
<img style="display: block; margin-left: auto; margin-right: auto" width="30%"
src="https://i.imgur.com/zCVajBa.png">
<font size=2> <center>圖片來源：https://zh.wikipedia.org/wiki/%E5%BC%A7%E5%BA%A6 </center></font>
這裡的角度 $\theta$ 一般可以用度（$\rm deg$）來作為單位，但是在Python中使用更有意義的單位稱為弧度（$\rm rad$）。弧度的定義為，圓心角所對應張開的弧長除以半徑，上圖是 $1 \rm rad$ 的定義因此一個圓周角即為 弧度，因此一般熟知的圓周角 $360\ {\rm deg} = 2\pi\ {\rm rad}$。
同學們知道了三角函數的定義之後，就可以了解在平面上如果有一向量 $\vec{r}$，則此向量可以透過其投影在 $x$ 軸和 $y$ 軸的分量來描述，一般使用直角座標表示法寫為$\left ( r\ {\rm sin}\theta , r\ {\rm cos}\theta \right )$，而在程式中我們就是用相同的方法來描述空間向量，只是更多了 $z$ 方向的空間維度。
從以上的座標圖，你也將發現 $\sqrt{x^2+y^2}=r$，這就是圓方程式，從圓心到圓弧上一點的距離是固定的。所以我們可以很方便的用三角函數去描述一些固定大小，但是不同方向的物理量。例如：在拋體運動中，如果固定拋射速率，只改變不同的拋射角度，則我們就會將初速度向量寫為 $\vec{v_{0}}=\left ( v_0\ {\rm sin}\theta , v_0\ {\rm cos}\theta, 0 \right )$，意思是將物體以初速 $v_{0}$，在 $x-y$ 平面上以仰角為 $\theta$ 拋出。或許同學可以透過迴圈控制計算出仰角為何的拋體，會具有最大的水平射程。同學們可以試著練習看看，利用Python算算看是否和你的經驗吻合。
### 指定作業 2-2
固定初速 $v_{0}$，拋射角度由 $\theta=3\ {\rm deg}$ 開始，每次完成拋射增加 $3\ {\rm deg}$ ：
1. 請在物件視窗，做出多個拋體運動並留下軌跡，同時也使用label指令在畫面上產生一個標籤顯示每次拋射時的角度，並將水平射程R、飛行時間T對上角度theta的關係，以graph指令呈現函數圖，如下所示：
:::info
label指令請參考 [vpython7 - label](http://www.glowscript.org/docs/VPythonDocs/label.html) 或 [vpython6 - label](http://vpython.org/contents/docs/label.html)。
:::
<img style="display: block; margin-left: auto; margin-right: auto" width="90%"
src="https://i.imgur.com/B4bGh4x.gif">
2. 請在執行視窗留下水平射程數據，並驗證 $\theta=45\ {\rm deg}$ 時具有最大水平射程。
<img style="display: block; margin-left: auto; margin-right: auto" width="90%"
src="https://i.imgur.com/g8kUgZk.png">
:::success
* 提示：
1. 初速改以三角函數表示其 $x$ 和 $y$ 方向的分量。
2. 利用條件判斷，在球落地後重新再以原速度大小發射一次。
3. 每次重新發射，善用迴圈累加技巧，使每次發射仰角增加3度。
```python=
while True:
ball.v = ball.v + ball.a*dt
...
if ball.pos.y <= size:
theta = theta + 3 * pi/180 #每次迴圈增加仰角3度
ball.pos = vector(0, size, 0) #位置重設
ball.v = vector(v0*cos(theta), v0*sin(theta), 0) #速度重設
...
```
4. label 文字標籤設定：
在 while 迴圈上面先產生物件：
```python
show_angle = label(pos=vector(0,-7*size,0), box = False, height = 20, color=color.yellow)
```
label 指令中的 box 是用來設定文字是否需要外框，False 為無框、True 為有框；opacity 是設定文字標籤的透明度，數值設定在 0~1 之間，0 表示完全透明、1 表示完全不透明；height為文字大小；text 為顯示在執行畫面上的文字，需用單引號設定為字串。
在 while 迴圈中將數值輸出到 label 中：
```python
show_angle.text = 'theta = %1.0f deg' %(theta/pi*180)
```
label 字串內容中的「%1.0f」表示將要顯示的數值取到整數位，如果是「%1.2f」表示取到小數點第二位，而這個標籤要引入的數值是寫在最後的「%(theta/pi*180))」，這是將theta換單位為deg後輸出。
:::
---
## 五、等速率圓周運動
同學們了解三角函數的定義後，可以馬上利用 $x$ 和 $y$ 方向的分量合成出圓軌跡，我們只要在迴圈中每一次增加一點小小的角度，迴圈持續的跑下去就會變成是動態的圓周運動。
現在我們做一個回轉半徑 $R=1$ 公尺，角速率大小 $\omega = 2\pi \ {\rm rad/s}$的等速率圓周運動。這裡的角速度大小指的是每單位時間掃過的角度，$2\pi \ {\rm rad/s}$代表每秒鐘繞圓周一圈。
### Exapmle 2 : 等速率圓周運動
```python=
from vpython import *
size = 0.1 #球的大小
theta = 0.0 #初始角度
R = 1.0 #圓周運動半徑
omega = 2*pi #角速度大小=單位時間繞過的角度
t = 0.0 #初始時間
scene = canvas(width=500, height=500, center=vector(0,0,0), background=vector(148.0/225,228.0/225,204.0/225))
ball = sphere(radius=size, color=color.blue, make_trail=True, interval = 1, retain = 900)
ball.pos = vector(R,0,0) #球的初始位置
t = 0.0 #初始時間
dt = 0.001 #時間間隔
while True:
rate(1/dt)
theta += omega*dt #t時刻的角度
ball.pos = vector(R*cos(theta), R*sin(theta), 0)
t += dt #計時器
```
:::info
* 為使物件視窗畫出來的運動軌跡較美觀不會重疊，這裡可以使用retain指令來固定軌跡點的數量，參考寫法如下：
ball = sphere(radius=size, color=color.blue, make_trail=True, interval = 1, <font color = "red">retain = 900</font>)
retain = 900：表示共留下900個軌跡點。
:::
到這邊我們已經能做出圓周運動的模擬了，再進一步可以利用這個圓周運動來設計計時功能，以下提供一個簡單的視覺化計時器範例。
### Example 3 : 三點記錄法 - 讓程式自動錄下依順序三個相鄰時間點
```python=
pre_theta = theta
while True:
rate(1/dt)
t = t + dt
pre_pre_theta = pre_theta #前前時刻的角度
pre_theta = theta #前一時刻的角度
theta = theta + omega*dt #現在時刻的角度
if pre_theta%(2*pi) > pre_pre_theta%(2*pi) and pre_theta%(2*pi) > theta%(2*pi):
print ('period = ', t)
#print pre_pre_theta%(2*pi), pre_theta%(2*pi), theta%(2*pi)
t = 0 #這行不寫t就不會歸零，而會印出t = 1, 2, 3, 4, ...
```
請同學「有腦的」複製 Example 3 的 code 貼到 Example 2 中，使物體回到原點的時候，都在執行視窗印出記錄數值。
#### 三點記錄法簡化版
同學可以將下面這小段程式執行看看，就能清楚電腦是怎麼錄下依序的三個時刻。
```python=
import time #python內建時鐘
t = 0 #初始時刻
dt = 1 #時間間隔
pre_t = t #只是為了防出錯，可以隨便設，第2個迴圈起就正常了
while True:
time.sleep(1) #每隔1秒跑一次迴圈
pre_pre_t = pre_t #前前時刻
pre_t = t #前一時刻
t = t + dt #現在時刻
print ('pre_pre_t =', pre_pre_t, ', pre_t =', pre_t, ', t =', t)
```
執行結果如下：
<img style="display: block; margin-left: 20px; margin-right: auto" width="38%"
src="https://i.imgur.com/2YCfre1.png">
:::info
* **視覺化計時器程式說明**
1. **如何紀錄現在時刻、前一時刻和前前時刻**：
同學一定要記得程式的執行是從上而下一行一行執行的：
前前時刻的角度： pre_pre_theta = pre_theta
前一時刻的角度： pre_theta = theta
現在時刻的角度： theta = theta + omega * dt
**++迴圈 1++：**
pre_theta被紀錄為theta，theta在迴圈結束時紀錄為**theta + omega * dt**
**++迴圈 2++：**
pre_pre_theta被紀錄為pre_theta，pre_theta被紀錄為前一迴圈的theta，即theta + omega * dt，theta在迴圈結束時被紀錄為**theta + omega * <font color = "red">2</font> * dt**
...
以下迴圈請同學依此類推，你將可以看出我們依時間順序紀錄三個時刻。
2. **餘數除法** %
在程式執行過程中，theta在超過 $2\pi\ {\rm rad}$ 或 $360 {\rm deg}$ 度之後，數值仍然會持續累加，所以我們沒法透過角度數值來判斷此圓周運動是否回到原點。
**「%」的功能說明**：如果圓周運動已經繞了$1$又$1/4$圈，繞過的角度累加為 $450 {\rm deg}$，但這角度對我們是無感的，因為看起來根本就是 $90$ 度。而迴圈並不會察覺這件事，而會一直累加下去，故而需要用到餘數除法，讓theta超過一圈後歸零。
同學們可以在執行視窗打入以下算式練習看看：
```python
450 % 360 = 90 # 1又1/4圈 和 1/4圈 意義是一樣的
540 % 360 = 180 # 1又1/2圈 和 1/2圈 意義是一樣的
```
3. **判斷回到原點的原理**
請參考下圖，利用剛剛紀錄下來的三個數據點，配合餘數除法，挑出物體回到點的時刻。
前前時刻： 359.7 % 360 = 359.7
前一時刻： **359.9 % 360 = 359.9** <font color = "red">...**三個記錄的角度最大值** </font>
現在時刻： 360.1 % 360 = 0.1
照理說，現在時刻的角度應該都要大於以前的時刻，但到達折返點時會發生<font color = "red">**前一時刻的角度是三個時刻中最大**</font>的現象，只要此條件判斷成立，即是此圓周運動回到原點的時刻。
<img style="display: block; margin-left: auto; margin-right: auto" width="60%"
src="https://i.imgur.com/Enj66Ph.png">
>
4. **回到原點的條件判斷指令**
```python
if pre_theta%(2*pi) > pre_pre_theta%(2*pi) and pre_theta%(2*pi) > theta%(2*pi):
# 物體跨過原點的條件判斷
```
:::
### 檢核作業 2-2：
請完成視覺化的圓周圓周運動計時器，並在執行視窗紀錄每繞一圈所花的時間。執行畫面如下圖所示，物件視窗中當物體跨過原點時，執行視窗就會打印出時間。
<img style="display: block; margin-left: auto; margin-right: auto" width="95%"
src="https://i.imgur.com/7jwrOK6.png">
:::success
* 提示：Example 3 貼入 Example 2 中。
:::
### 指定作業 2-3
延續++課堂作業 2-2++，請將此圓周運動用 cylinder 指令每隔一小段時間就畫出細線，將此圓周分割成 N 等分。執行畫面如下：
<img style="display: block; margin-left: auto; margin-right: auto" width="80%"
src="https://i.imgur.com/5Whrp1I.png">
上左圖為圓周運動繞第一圈時，這一圈用來計算週期不劃線，等到回到出發點瞬間，電腦將記錄下繞一圈的週期 period_t，然後再將此時間除以 N，於第二圈開始每隔 period_t/N 用 cylinder 指令畫一條線，如上右圖所示為切成 20 等分的樣子。
你可能會疑惑，週期不是一開始就知道，為什麼要大費周章繞完一圈讓電腦求出週期？在 Example 2 中我們預先設定了角速率為 $2 \pi$，也就是一秒繞一圈，那週期早就知道了不是嗎？
事實上，有很多週期性運動我們一開始並不知道週期，例如：行星受萬有引力繞太陽運行，我們沒法從初速和受力就心算出繞行週期，因此必須要等它繞回原點利用程式裡的計時器把週期測出來。這個技巧會在未來 Lecture 5 克普勒第二定律用上，事實上行星運動是橢圓軌道，你更難一開始就猜到週期究竟是多少。
以下我們提供建議的寫法：
:::success
* 提示：等時距分割週期性運動
1. 設定布林值判斷小工具：
在 while 迴圈上面設定如下指令：
```python
back = False
...
while True:
...
```
back 是我們取的「事件名稱」，用來描述「回到原點」這件事情，程式一開始運行時先把它設定為 False，因為球剛出發還沒有回到出發點。
2. 在迴圈中跨過原出發點的條件判斷裡更改 back 的布林值為 True：
```python
if pre_theta%(2*pi) > pre_pre_theta%(2*pi) and pre_theta%(2*pi) > theta%(2*pi):
print ('period = ', t)
back = True #如果跨過原點代表back這個事件發生了，所以把False改成True
period_t = t #用一個新的變數period_t把週期存下來
t = 0
```
這裡為什麼要多設一個新的變數 period_t？是因為這一組條件判斷結束時 t 會被歸零，故要在 t 被歸零前把計時器的數據存出來當作週期。
3. 在迴圈中 2. 跨過原出發點條件判斷的**上面**加入以下程式：
```python
if back:
plot_t = t % (period_t/N) #將週期切成N等分，並用餘數除法設定畫圖時間點
if plot_t + dt >= ... and plot_t < ... : #畫圖時間判斷點
cylinder(radius=size/50, color=color.black, pos= ... , axis= ... ) #畫細線
```
當 back 為 True 時，程式就會進入這一段條件結構然後畫圖，至於畫圖的技巧與 Lecture 1 進階作業 1-2 相同。
4. 第三圈開始不畫線：你可以再想想看，在第三圈之後怎麼讓程式不再劃線，因為你會發現之後每一圈畫的線都會跟前面的重疊，也沒必要浪費電腦空間繼續畫下去。甚至你會發現跑了很多圈之後，因為誤差的累加，每一圈畫的線會有細微的錯開，時間久了線條看起來會有點暈開的感覺。這就留給大家想一想，暫不提供解法給大家了！
:::
---
## 六、自由落體運動的等時距取點（輸出及載入MS excel檔）-選讀教材
為了方便使模擬的數據與實驗數據對照，我們可以透過 xlwt 套件將模擬數據輸出成 MS excel 檔。範例程式如下：
```python=
from vpython import *
import xlwt
g = 9.8 #重力加速度 9.8 m/s^2
size = 0.5 #球半徑 0.5 m
height = 15.0 #球初始高度 15 m
m = 1.0 #球質量1kg
Fg = vector(0, -m*g, 0) #重力
book = xlwt.Workbook(encoding="utf-8")
sheet1 = book.add_sheet("Sheet 1", cell_overwrite_ok=True)
scene = canvas(width=600, height=600,x=0, y=0, center = vector(0,height/2,0)) #設定畫面
floor = box(length=20, height=0.01, width=10, color=color.green) #畫地板
ball = sphere(radius = size, color=color.yellow, make_trail= True, trail_type="points", interval=100) #畫球
ball.pos = vector(0, height, 0) #球初始位置
ball.v = vector(0, 0, 0) #球初速
dt = 0.001 #時間間隔0.001秒
t = 0.0 #模擬初始時間為0秒
T = 0.1 #每隔T秒紀錄一組數據至xls檔
q = 1 #xls檔的列數計數
while ball.pos.y > size: #模擬直到球落地 即y=球半徑
rate(1/dt) #每一秒跑 1000 次
t = t + dt #計時器
ball.a = Fg/m #球的加速度
ball.v = ball.v + ball.a*dt #球的末速度 = 前一刻速度 + 加速度*時間間隔
ball.pos = ball.pos + ball.v * dt #球的末位置 = 前一刻位置 + 速度*時間間隔
#寫入表頭名稱
sheet1.write(0,1,'時間')
sheet1.write(0,2,'y方向位置')
sheet1.write(0,3,'y方向速度')
sheet1.write(0,4,'y方向加速度')
#將數據寫入xls檔
if t%T < T and t%T + dt >= T: #設定寫入數據時機
sheet1.write(q,1,t)
sheet1.write(q,2,ball.pos.y)
sheet1.write(q,3,ball.v.y)
sheet1.write(q,4,ball.a.y)
q = q + 1 #計數器：每寫完一次數據q增加1，使下次數據在xls檔中往下一列的欄位寫入
book.save("free fall data.xls")
print (t, ball.v)
```
輸出之後的 MS excel 表格如下：
<img style="display: block; margin-left: auto; margin-right: auto" width="50%"
src="https://i.imgur.com/KT278dY.png">
接著可以利用 MS excel 中的「插入圖表」功能，這些數據繪製成「位置-時間」、「速度-時間」、「加速度-時間」的運動學分析圖，如下所示：
<img style="display: block; margin-left: auto; margin-right: auto" width="85%"
src="https://i.imgur.com/cjgSyAt.png">
若有實驗數據，則可以直接在 MS excel 紀錄數據之後，一樣輸出為如上圖形，即可立即對比實驗數據與模擬數據的吻合度，讓實驗與理論的驗證更為直接、方便！
另外，如果要將獲得的實驗數據載入 python 中進行運算與繪圖，則需要使用到 xlrd 這個套件來讀取 MS excel 檔。以下我們以方才輸出的檔案（free fall data.xls）進行載入的練習，範例程式如下：
```python=
from vpython import *
import xlrd
size = 0.5
height = 15.0
filename = 'free fall data.xls' #xls檔案名稱
book = xlrd.open_workbook(filename) #以xlrd打開表格給book
sheel_1 = book.sheet_by_index(0) #獲取index=0的工作表，即第一張工作表
scene = canvas(width=600, height=600,x=0, y=0, center = vector(0,height/2,0)) #設定畫面
floor = box(length=20, height=0.01, width=10, color=color.green) #畫地板
ball = sphere(radius = size, color=color.yellow, make_trail= True, trail_type="points", interval=1) #畫球
i = 1 #xls檔中的列數記數
while i < sheel_1.nrows:
time.sleep(0.1) #每0.1秒載入一次數據
t = sheel_1.cell_value(i,1) #將時間欄位指定給時間變數t
ball.pos.y = sheel_1.cell_value(i,2) #將位置欄位指定給球的y方向位置
i = i + 1 #迴圈結束前計數器+1，以讀取下一個欄位
```
執行後你將會發現，xls 檔中的數據完全被還原到 vpython 的繪圖環境中。也就是說，未來你可以將真實實驗測得的數據，在 vpython 中還原，做到「虛實整合」的效果。不過我們剛剛採點數據點較少，因此還原之後的模擬過程並不細緻，所以要將實驗數據較真實還原的話，就必須要密集的採點。
實驗上，較方便的方法可以使用攝影機錄下實驗過程，再使用類似像 tracker 之類的影像分析軟體，將時間與座標點的數據取出，如果每個影片畫格都採點的話，就可以非常細緻的在 vpython 的環境中擬真呈現實驗。但是採點過度密集可能造成讀取速度降低，而產生模擬畫面的延遲，所以「適量」就好！
以上簡單提供如何利用 python 寫入與讀取 MS excel 檔，若想知道更多可以參考以下連結： https://www.itread01.com/articles/1498941489.html 。
---
# 作業繳交：
平台連結：Joy of Code : http://203.64.158.237/student/lessons/D/
此平台由南港高中-高慧君老師設計建置
Sign in with Wallet
Connect another wallet