# 使用 For 迴圈計算水平抛射資料
> 作者：王一哲
> 日期：2018/3/28
<br />
在前一篇文章〈[水平抛射](https://hackmd.io/s/Hy4UoZfMm)〉中，我們寫了一個模擬水平抛射運動的程式，再手動修高度 h 計算飛行時間 t 及水平射程 R。但是手動修改是很麻煩的，如果能多寫一個 for 迴圈，就能讓程式幫我們自動代入不同的 h 並將資料存成文字檔。以下共有兩個程式：
- 程式5-4：水平抛射, 改變h, 記錄R （[GlowScript 網站動畫連結](http://www.glowscript.org/#/user/yizhe/folder/Public/program/5-4projectileh)）
- 程式5-5：由程式5-4修改, 用for迴圈改變h, 記錄t、R
<br />
## 程式 5-4：水平抛射, 改變h, 記錄R （[取得程式碼](https://github.com/YiZheWangTw/VPythonTutorial/blob/master/05.%E6%B0%B4%E5%B9%B3%E6%8A%9B%E5%B0%84/5-4_projectile_h.py)）
```python=
"""
VPython教學: 5-4.水平抛射, 改變h, 記錄R
Ver. 1: 2018/2/19
Ver. 2: 2019/9/6
作者: 王一哲
"""
from vpython import *
"""
1. 參數設定, 設定變數及初始值
"""
size = 1 # 小球半徑
v0 = 5 # 小球水平初速
h = 15 # 小球離地高度5, 10, 15, 20, 25, 30, 35, 40, 45, 50
L = 50 # 地板長度
g = 9.8 # 重力加速度 9.8 m/s^2
t = 0 # 時間
dt = 0.001 # 時間間隔
"""
2. 畫面設定
"""
scene = canvas(title="Projectile", width=800, height=600, x=0, y=0, center=vec(0, h/2, 0), background=vec(0, 0.6, 0.6))
floor = box(pos=vec(0, -size, 0), size=vec(L, 0.01, 10), texture=textures.metal)
ball = sphere(pos=vec(-L/2, h, 0), radius=size, texture=textures.wood, make_trail=True, v=vec(v0, 0, 0), a=vec(0, -g, 0))
"""
3. 物體運動部分, 小球觸地時停止
"""
while ball.pos.y - floor.pos.y > size + 0.5*floor.height:
rate(1000)
ball.v += ball.a*dt
ball.pos += ball.v*dt
t += dt
print(t, ball.pos.x + L/2)
```
<br />
這個程式是程式 5-1 的簡化版，由於我只想算出飛行時間 t 及水平射程 R，小球碰到地板後就可以停止動畫，因此 while 迴圈中設定的條件為
```python
ball.pos.y - floor.pos.y > size + 0.5*floor.height
```
手動修改 h，分別以 5、10、15、……50代入，再將計算得到的 t 和 R 手動輸入 SciDAVis 並作圖，這個作法雖然可行，但是相當費時。
<br />
### 數據處理原則
通常我們只用 XY 散布圖，而且只畫數據點，不會將數據點連線，畫完數據點之後再加上擬合直線式曲線。如果應變數與自變數的關係圖為曲線，還要想辦法化為直線，最後的成果一定要用線性擬合，還要解釋最接近直線的斜率、y軸截距的物理意義，要求相當嚴謹。使用 SciDAVis 畫 XY 散布圖的方法請參考另一篇文章〈[SciDAVis繪圖：XY散布圖](https://hackmd.io/s/HkNHsKFzN)〉。
先畫 t - h 關係圖，由於圖形看起來像是二次曲線，就試著用多項式擬合，最高次方為二次，結果為
$$ y = 0.656 + 0.070 x - 0.00058 x^2$$
但這畢竟是猜測的結果，不一定正確。比較好的作法是取 \\( \log t \\) 、\\( \log h \\)，畫 \\( \log t - \log h \\) 關係圖，畫出來的圖為斜直線，最接近直線為
$$ y = 0.500 x - 0.345 $$
由斜率可以猜測 \\( t \propto \sqrt h \\)。理論上
$$ h = \frac{1}{2}gt^2 ~\Rightarrow~ \log h = \log \left( \frac{1}{2}gt^2 \right) = \log \left( \frac{1}{2}g \right) + \log t^2 = 2 \log t + c$$ $$ \log t = 0.5 \log h + c' $$
為了驗證這個猜測，新增一欄數據計算 \\( \sqrt h \\)，畫 \\( t - \sqrt h \\) 關係圖，畫出來的圖為斜直線，最接近直線為
$$ y = 0.45178 x - 0.00034 $$
理論上
$$ h = \frac{1}{2}gt^2 ~\Rightarrow~ t = \sqrt{\frac{2h}{g}} ~\Rightarrow~ slope = \sqrt{\frac{2}{g}} \approx 0.45175 $$
理論值與模擬得到的值很接近，可以確定程式應該沒有問題。
<br />
<img style="display: block; margin-left: auto; margin-right: auto" height="100%" width="100%" src="https://i.imgur.com/SkGwcNE.png">
<div style="text-align:center">於 SciDAVis 手動輸入數據</div>
<br />
<img style="display: block; margin-left: auto; margin-right: auto" height="100%" width="100%" src="https://imgur.com/U1zNT9t.jpg">
<div style="text-align:center">t - h 關係圖</div>
<br />
<img style="display: block; margin-left: auto; margin-right: auto" height="100%" width="100%" src="https://i.imgur.com/Emntfzv.jpg">
<div style="text-align:center">log(t) - log(h) 關係圖</div>
<br />
<img style="display: block; margin-left: auto; margin-right: auto" height="100%" width="100%" src="https://i.imgur.com/xR0iB99.jpg">
<div style="text-align:center">t - h<sup>1/2</sup> 關係圖</div>
<br /><br />
## 程式 5-5：水平抛射, 用for迴圈改變h, 記錄t、R （[取得程式碼](https://github.com/YiZheWangTw/VPythonTutorial/blob/master/05.%E6%B0%B4%E5%B9%B3%E6%8A%9B%E5%B0%84/5-5_projectile_loop.py)）
```python=
"""
VPython教學: 5-5.水平抛射, 用for迴圈改變h, 記錄t、R
Ver. 1: 2018/3/21
Ver. 2: 2018/6/13 加上用 matplotlib.pyplot 繪圖的部分
Ver. 3: 2019/10/10 改用 with 開啟檔案, 上課用的版本
作者: 王一哲
"""
from vpython import *
import matplotlib.pyplot as plt
"""
1. 參數設定, 設定變數及初始值
"""
size = 1 # 小球半徑
v0 = 5 # 小球水平初速
L = 50 # 地板長度
g = 9.8 # 重力加速度 9.8 m/s^2
dt = 0.001 # 時間間隔
data_y, data_x = [], [] # 儲存繪圖資料的串列
"""
2. 畫面設定
"""
scene = canvas(title="Projectile", width=800, height=600, x=0, y=0, center=vec(0, 5, 0), background=vec(0, 0.6, 0.6))
floor = box(pos=vec(0, -size, 0), size=vec(L, 0.01, 10), texture=textures.metal)
# 開啟檔案 5-5_data.csv, 屬性為寫入, 先寫入欄位的標題
with open("5-5_data.csv", "w", encoding="UTF-8") as file:
file.write("h(m), t(s), R(m)\n")
"""
3. 主程式
"""
def main(h):
t = 0
ball = sphere(pos=vec(-L/2, h, 0), radius=size, texture=textures.wood, make_trail=True,
v=vec(v0, 0, 0), a=vec(0, -g, 0))
while ball.pos.y - floor.pos.y > size + 0.5*floor.height:
rate(1000)
ball.v += ball.a*dt
ball.pos += ball.v*dt
t += dt
return t, ball.pos.x + L/2
"""
4. 用for迴圈改變h, 計算t、R, 寫入檔案
"""
for h in range(5, 51, 1):
t, r = main(h)
with open("5-5_data.csv", "a", encoding="UTF-8") as file:
file.write(str(h) + "," + str(t) + "," + str(r) + "\n")
data_y.append(t)
data_x.append(h)
# 用 matplotlib.pyplot 繪圖
plt.figure(figsize=(8, 6), dpi=100) # 設定圖片尺寸
plt.xlabel(r'$h ~\mathrm{(m)}$', fontsize=16) # 設定坐標軸標籤
plt.ylabel(r'$t ~\mathrm{(s)}$', fontsize=16)
plt.xticks(fontsize=12) # 設定坐標軸數字格式
plt.yticks(fontsize=12)
plt.grid(color='red', linestyle='--', linewidth=1) # 設定格線顏色、種類、寬度
plt.plot(data_x, data_y, marker='o', markerfacecolor='blue', markersize=8, color='skyblue', linewidth=3) # 繪圖並設定線條顏色、寬度
plt.savefig('t-h_plot.svg') # 儲存圖片
plt.savefig('t-h_plot.png')
plt.show() # 顯示圖片
```
<br />
<img style="display: block; margin-left: auto; margin-right: auto" height="100%" width="100%" src="https://i.imgur.com/OyULerj.png">
<div style="text-align:center">程式 5-5 模擬結果畫面截圖</div>
<br />
**程式 5-5 當中運用了串列 (list) 及自訂函式 (def) 兩個新的物件，詳細的說明請參考另外兩篇講義〈[串列](https://hackmd.io/@yizhewang/H1rSL1W8K)〉、〈[自訂函式](https://hackmd.io/@yizhewang/S1awsJ-IY)〉** ，以下先說明與程式 5-4 的不同之處。
1. 為了將資料儲存到文字檔中，需要加上以下的程式碼：
```python
with open("5-5_data.csv", "w", encoding="UTF-8") as file
file.write("h(m), t(s), R(m)\n")
```
第一行程式碼用來開啟名稱為 5-5_data.csv 的文字檔，如果檔案已存在則覆蓋檔案內容，文字編碼為 UTF-8。第二行程式碼會在 5-5_data.csv 中寫入文字 h(m), t(s), R(m) 並換行。
2. 第 34 ~ 42 行，用 def 自訂名為 main 的函式，輸入變數 h 數值，自動跑完小球落下的過程，回傳落下需要的時間 t 以及水平射程 ball.pos.x + L/2 。如果要讓動畫不要執行地太快，要在 while 迴圈當中加上 rate(1000)。
3. 為了讓程式自動改變 h 的數值並代入後續的程式碼中運算，需要用到 for 迴圈。
```python
for h in range(5, 51, 1):
t, r = main(h)
with open("5-5_data.csv", "a", encoding="UTF-8") as file:
file.write(str(h) + "," + str(t) + "," + str(r) + "\n")
data_y.append(t)
data_x.append(h)
```
h 的數值由 5 開始代入，每次增加 1，直到數量為 50 為止。將每個小球的飛行時間及水平射程回傳為 t 及 r，再用
```python
file.write(str(h) + "," + str(t) + "," + str(r) + "\n")
```
將 h、t、r 轉為字串、用逗號分隔、寫入檔案。
4. 由於是用 with 開啟檔案，不需要加上 file.close() 關閉檔案。
5. 儲存到文字檔中的資料格式如下（[取得檔案](https://raw.githubusercontent.com/YiZheWangTw/VPythonTutorial/master/05.%E6%B0%B4%E5%B9%B3%E6%8A%9B%E5%B0%84/data.csv)）：
```python
h(m), t(s), R(m)
5,1.0099999999999996,5.049999999998995
6,1.105999999999989,5.5299999999989
7,1.1949999999999792,5.974999999998811
8,1.2769999999999702,6.38499999999873
9,1.3549999999999616,6.774999999998652
10,1.4279999999999535,7.1399999999985795
```
6. 第 55 ~ 64 行的程式碼是利用 matplotlib 套件繪圖，如果想要用其它軟體讀取 5-5_data.csv 的資料並繪圖，可以刪除這幾行。
7. 由於 [GlowScript](https://www.glowscript.org/#/) 網站上不能將資料儲存成文字檔，也不支援 matplotlib，可以使用以下的程式碼，將不支援的部分刪掉，並用 print 輸出需要的資料。程式碼第1行 **Web VPython 3.2** 是在網站上建立程式時自動産生的，會自動採用最新的版本，不需要手動輸入，也不需要再寫 **from vpython import \***，這是 [GlowScript 動畫連結](https://www.glowscript.org/#/user/yizhe/folder/Public/program/5-5projectliLoop)。
```python=
Web VPython 3.2
"""
1. 參數設定, 設定變數及初始值
"""
size = 1 # 小球半徑
v0 = 5 # 小球水平初速
L = 50 # 地板長度
g = 9.8 # 重力加速度 9.8 m/s^2
dt = 0.001 # 時間間隔
"""
2. 畫面設定
"""
scene = canvas(title="Projectile", width=800, height=600, x=0, y=0, center=vec(0, 5, 0), background=vec(0, 0.6, 0.6))
floor = box(pos=vec(0, -size, 0), size=vec(L, 0.01, 10), texture=textures.metal)
"""
3. 主程式
"""
def main(h):
t = 0
ball = sphere(pos=vec(-L/2, h, 0), radius=size, texture=textures.wood, make_trail=True,
v=vec(v0, 0, 0), a=vec(0, -g, 0))
while ball.pos.y - floor.pos.y > size + 0.5*floor.height:
rate(1000)
ball.v += ball.a*dt
ball.pos += ball.v*dt
t += dt
return t, ball.pos.x + L/2
"""
4. 用for迴圈改變h, 計算t、R, 寫入檔案
"""
for h in range(5, 51, 1):
t, r = main(h)
print("{}, {}, {}".format(h, t, r))
```
<img style="display: block; margin-left: auto; margin-right: auto" height="100%" width="100%" src="https://imgur.com/sHfjAis.png">
<img style="display: block; margin-left: auto; margin-right: auto" height="100%" width="100%" src="https://imgur.com/jjsIYXs.png">
<div style="text-align:center">於 GlowScript 網站執行的結果</div>
<br />
<br />
### 數據處理原則
仿照程式 5-4 的方法，先畫 t - h 關係圖，多項式擬合的結果為
$$ y = 0.728 + 0.074 x - 0.0005 x^2 $$
再畫 \\( \log t - \log h \\) 關係圖，畫出來的圖為斜直線，最接近直線為
$$ y = 0.500 x - 0.345 $$
由斜率可以猜測 \\( t \propto \sqrt h \\)。接著新增一欄數據計算 \\( \sqrt h \\)，畫 \\( t - \sqrt h \\) 關係圖，畫出來的圖為斜直線，最接近直線為
$$ y = 0.45179 x - 0.00043 $$
斜率的理論值為 0.45175，理論值與模擬得到的值很接近，可以確定程式應該沒有問題。
<img style="display: block; margin-left: auto; margin-right: auto" height="100%" width="100%" src="https://imgur.com/oQp10Ta.jpg">
<div style="text-align:center">t - h 關係圖</div>
<br />
<img style="display: block; margin-left: auto; margin-right: auto" height="100%" width="100%" src="https://i.imgur.com/1TTOupe.jpg">
<div style="text-align:center">log(t) - log(h) 關係圖</div>
<br />
<img style="display: block; margin-left: auto; margin-right: auto" height="100%" width="100%" src="https://i.imgur.com/0xWTrqt.jpg">
<div style="text-align:center">t - h<sup>1/2</sup> 關係圖</div>
<br />
<img style="display: block; margin-left: auto; margin-right: auto" height="100%" width="100%" src="https://imgur.com/BnQsnlH.png">
<div style="text-align:center">使用 Matplotlib 繪製的 t - h 關係圖</div>
<br /><br />
## 結語
本文是以改變 h 為例，可以試著用相同的方法改變其它的變數值，例如改變空氣阻力的係數，找出飛行時間與空氣阻力係數的關係。而且用 for 迴圈自動改變數值，再將運算結果存入文字檔中，比手動改變數值所需要的時間少很多，這樣才能發揮寫程式的好處。
<br /><br />
## VPython官方說明書
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. **texture**: http://www.glowscript.org/docs/VPythonDocs/textures.html
<br />
---
###### tags: `VPython`
Sign in with Wallet
Connect another wallet