# 斜向抛射

> 作者：王一哲
> 日期：2018/3/29

<br />

將一個小球由地面以初速 v₀、仰角 𝜃 抛出，小球受到重力作用向下加速，計算小球的飛行時間及水平射程。分為以下3種不同的狀況：

1. 只考慮重力的作用，小球撞到地板時停止運動。　（[GlowScript 網站動畫連結](http://www.glowscript.org/#/user/yizhe/folder/Public/program/6-1projection)）
2. 同時考慮重力及空氣阻力 \\( f = -bv \\)，同時畫出考慮空氣阻力與不考慮空氣阻力的小球。　（[GlowScript 網站動畫連結](http://www.glowscript.org/#/user/yizhe/folder/Public/program/6-2projectiona)）
3. 使用 for 迴圈，分別計算不同仰角 𝜃、不同的空氣阻力係數 b 對應的飛行時間和水平射程。　（[GlowScript 網站動畫連結](http://www.glowscript.org/#/user/yizhe/folder/Public/program/6-3projectionloop)）

成果如下：

<img style="display: block; margin-left: auto; margin-right: auto" height="100%" width="100%" src="https://i.imgur.com/Xcz900w.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/8QiHVSk.png">
<div style="text-align:center">同時考慮重力及空氣阻力，b = 0.1</div>
<br />

<img style="display: block; margin-left: auto; margin-right: auto" height="100%" width="100%" src="https://i.imgur.com/kudsA9l.png">
<div style="text-align:center">使用 for 迴圈改變仰角 𝜃 = 15º ~ 90º</div>
<br />

<img style="display: block; margin-left: auto; margin-right: auto" height="100%" width="100%" src="https://i.imgur.com/Jt67aCO.png">
<div style="text-align:center">使用 for 迴圈改變空氣阻力係數 b = 0 ~ 0.9</div>
<br />

## 程式 6-1：斜向抛射, 球落地停止運作 （[取得程式碼](https://github.com/YiZheWangTw/VPythonTutorial/blob/master/06.%E6%96%9C%E5%90%91%E6%8A%9B%E5%B0%84/6-1_projection.py)）

"""
 VPython教學: 6-1.斜向抛射, 球落地停止運作
 Ver. 1: 2018/2/21
 Ver. 2: 2019/9/7
 作者: 王一哲
"""
from vpython import *

"""
 1. 參數設定, 設定變數及初始值
"""
size = 1              # 小球半徑
v0 = 30               # 小球初速
theta = radians(30)   # 小球抛射仰角, 用 radians 將單位轉為弧度
L = 100               # 地板長度
g = 9.8               # 重力加速度 9.8 m/s^2
t = 0                 # 時間
dt = 0.001            # 時間間隔

"""
 2. 畫面設定
"""
scene = canvas(title="Projection", width=800, height=400, 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)
ball = sphere(pos=vec(-L/2, 0, 0), radius=size, color=color.red, make_trail=True,
              v=vec(v0*cos(theta), v0*sin(theta), 0), a=vec(0, -g, 0))

"""
 3. 物體運動部分
"""
while ball.pos.y - floor.pos.y >= size:
    rate(1000)
    ball.v += ball.a*dt
    ball.pos += ball.v*dt
    t += dt

print(t, ball.pos.x + L/2)

<br />

### 參數設定

在此定義的變數有 size、v0、theta、L、g、t、dt，用途都已經寫在該行的註解中。這裡出現了一個新的函式

radians(x)

這個函式可以將 x 的單位由角度轉換為弧度(rad)，radians 包含於 math 函式庫中，使用前需要引入函式庫，但是 vpython 已經包含了 math，不需要額外再引入 math。需要這樣做是因為 python 的三角函數輸入值皆以 rad 為單位。另一個類似的函式為

degrees(x)

這個函式可以將 x 的單位由 rad 轉換為角度。
<br />

### 畫面設定

畫面設定部分與程式 5-1 非常類以，不同之處在於小球的初位置及初速

ball = sphere(pos=vec(-L/2, 0, 0), radius=size, color=color.red, make_trail=True,
              v=vec(v0*cos(theta), v0*sin(theta), 0), a=vec(0, -g, 0))

小球出發時高度為0，初速度 \\( v_x = v_0 \cos \theta \\)、\\( v_y = v_0 \sin \theta \\)。
<br />

### 物體運動

物體運動部分與程式 5-1 非常類以，不同之處在於小球碰到地板時停止動畫，因此 while 迴圈中設定的條件為

ball.pos.y - floor.pos.y >= size

在 while 迴圈結束之後，印出飛行時間 t 及水平射程 R。可以手動修改 v0 或 theta 的數值，觀察 t 和 R 的變化。
<br />

## 程式 6-2.斜向抛射, 球落地停止運作, 有空氣阻力 （[取得程式碼](https://github.com/YiZheWangTw/VPythonTutorial/blob/master/06.%E6%96%9C%E5%90%91%E6%8A%9B%E5%B0%84/6-2_projection_a.py)）

"""
 VPython教學: 6-2.斜向抛射, 球落地停止運作, 有空氣阻力
 Ver. 1: 2018/2/21
 Ver. 2: 2019/11/13 改變儲存球落地時間的方法
 作者: 王一哲
"""
from vpython import *

"""
 1. 參數設定, 設定變數及初始值
"""
size = 1              # 小球半徑
m = 1                 # 小球質量
v0 = 30               # 小球初速
theta = radians(30)   # 小球抛射仰角, 用 radians 將單位轉為弧度
L = 100               # 地板長度
b = 0.1               # 空氣阻力 f=-bv
g = 9.8               # 重力加速度 9.8 m/s^2
t = 0                 # 時間
dt = 0.001            # 時間間隔
t1, t2 = 0, 0         # 小球落地時間
s1, s2 = False, False # 小球是否落地

"""
 2. 畫面設定
"""
scene = canvas(title="Projection with Air Resistance", width=800, height=400,
               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)
# ball: 有空氣阻力; ball2: 沒有空氣阻力
ball = sphere(pos=vec(-L/2, 0, 0), radius=size, color=color.blue, make_trail=True,
              v=vec(v0*cos(theta), v0*sin(theta), 0))
ball2 = sphere(pos=vec(-L/2, 0, 0), radius=size, color=color.red, make_trail=True,
               v=vec(v0*cos(theta), v0*sin(theta), 0), a=vec(0, -g, 0))

"""
 3. 物體運動部分, 小球觸地停止運動, 記錄落地時間, 畫出兩者軌跡
    印出落地時間及水平射程
"""
while ball.pos.y - floor.pos.y >= size or ball2.pos.y - floor.pos.y >= size:
    rate(1000)
    # ball 還沒有落地才更新加速度、速度、位置
    if s1 == False:
        f = -b*ball.v
        ball.a = vec(0, -g, 0) + f/m
        ball.v += ball.a*dt
        ball.pos += ball.v*dt
    # ball2 還沒有落地才更新加速度、速度、位置
    if s2 == False:
        ball2.v += ball2.a*dt
        ball2.pos += ball2.v*dt
    # ball 第一次碰到地板時將狀態改成 True 並記錄時間
    if ball.pos.y - floor.pos.y <= size and s1 == False:
        ball.v = vec(0, 0, 0)
        s1, t1 = True, t
    # ball2 第一次碰到地板時將狀態改成 True 並記錄時間
    if ball2.pos.y - floor.pos.y <= size and s2 == False:
        ball2.v = vec(0, 0, 0)
        s2, t2 = True, t
    # 更新時間
    t += dt

print("t1, r1, t2, r2")
print("{:f}, {:f}, {:f}, {:f}".format(t1, ball.pos.x + L/2, t2, ball2.pos.x + L/2))

<br />

### 參數設定

程式 6-2 與 6-1 很像，但是增加了小球質量 m 及空氣阻力係數 b。為了記錄小球飛行時間，新增了 s1、s2、t1、t2 等4個變數。如果對於加上空氣阻力的寫法已經很熟悉，其實可以跳過 6-1 直接寫 6-2。
<br />

### 畫面設定

為了與不考慮空氣阻力的理想狀況對照，需要畫出兩個球，ball 是考慮空氣阻力的狀況、ball2 是理想狀況。由於 ball 的加速度與速度有關，放在 while 迴圈裡再設定即可。
<br />

### 物體運動

1. 為了讓兩個小球都碰到地板後才停止動畫，因此 while 迴圈中設定的條件為

    ball.pos.y - floor.pos.y >= size or ball2.pos.y - floor.pos.y >= size

2. 由於 ball 需要考慮空氣阻力，需要加上

    f = -b*ball.v
    ball.a = vector(0, -g, 0) + f/m

用原來的速度計算空氣阻力，再代入 \\( F = ma \\) 更新加速度。

3. 為了記錄飛行時間，先將狀態 s1、s2 設定成 False，當小球第一次碰到地板時將狀態改為 True 並將時間 t 儲存到 t1 及 t2。
<br />

## 程式 6-3.斜向抛射, 使用for 迴圈改變仰角 theta, 空氣阻力係數 b （[取得程式碼](https://github.com/YiZheWangTw/VPythonTutorial/blob/master/06.%E6%96%9C%E5%90%91%E6%8A%9B%E5%B0%84/6-3_projection_loop.py)）

"""
 VPython教學: 6-3.斜向抛射, 使用 for 迴圈改變仰角 theta, 空氣阻力係數 b
 Ver. 3: 2019/10/11 改用 with 及 def 簡化程式
 作者: 王一哲
"""
from vpython import *

"""
 1. 參數設定, 設定變數及初始值
"""
size = 1    # 小球半徑
m = 1       # 小球質量
L = 100     # 地板長度
g = 9.8     # 重力加速度 9.8 m/s^2
dt = 0.001  # 時間間隔

"""
 2. 畫面設定
"""
scene = canvas(title="Projection with for loop", width=800, height=400, 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)

# 開啟檔案 data.csv, 屬性為寫入, 先寫入欄位的標題
with open("data.csv", "w", encoding="UTF-8") as file:
    file.write("theta(degree), b, t(s), R(m)\n")

"""
 3. 物體運動部分
"""
def main(i=0, n=16, v0=30, degree=30, b=0.0):
    t = 0
    theta = radians(degree)
    ball = sphere(pos=vec(-L/2, 0, 0), radius=size, color=vec((n-i)/n, 0, i/n),
                  make_trail=True, v=vec(v0*cos(theta), v0*sin(theta), 0))
    while ball.pos.y - floor.pos.y >= size:
        rate(1000)
        f = -b*ball.v
        ball.a = vec(0, -g, 0) + f/m
        ball.v += ball.a*dt
        ball.pos += ball.v*dt
        t += dt
    return degree, b, t, ball.pos.x + L/2

for i in range(0, 10): # 改變 theta 時設為16, 改變 b 時設為10
    #degree = 15 + 5*i  # 改變 theta
    #degree, b, t, r = main(i=i, degree=degree)
    b = 0.1*i          # 改變 b
    degree, b, t, r = main(i=i, n=10, b=b)
    print("{:f}, {:f}, {:f}, {:f}".format(degree, b, t, r))
    with open("data.csv", "a", encoding="UTF-8") as file:
        file.write(str(degree) + "," + str(b) + "," + str(t) + "," + str(r) + "\n")

<br />

### 程式設計部分

程式 6-3 是將 5-5、6-1、6-2 合為一體的產物，有兩個使用方法：

1. 主程式 main 中需要輸入的參數有 i、n、v0、degree、b，如果引用 main 時不輸入數值則採用等號後的預設值。

    main(i=0, n=16, v0=30, degree=30, b=0.0)

2. 改變仰角 theta，範圍為 15º ~ 90º，每次增加 5º，總共執行 16 次。為了記錄代入的角度數值，需要增加變數 degree，將 i 及 degree 代入主程式 main 並回傳計算結果。

    degree = 15 + 5 * i
    degree, b, t, r = main(i=i, degree=degree)

3. 改變空氣阻力係數 b，範圍為 0 ~ 0.9，每次增加 0.1，總共執行 10 次，將 i、n 及 degree 代入主程式 main 並回傳計算結果。

    b = 0.1 * i
    degree, b, t, r = main(i=i, n=10, b=b)

4. 為了方便使用者從圖中看出執行的順序，我將小球的顏色加上一點變化

    color = vec((n-i)/n, 0, i/n)

一開始小球是紅色的，隨著執行次數增加會慢慢變為藍色。

<br />

### 數據處理部分

若以 v0 = 30 m/s、b = 0、改變 theta 產生的數據如下

v0(m/s), theta(degree), b, t(s), R(m)
30,15,0.0,1.5839999999999363,45.900795265255965
30,20,0.0,3.677999999999706,59.03149043776981
30,25,0.0,6.265000000000427,70.3385473519133
30,30,0.0,9.32600000000027,79.5271128295263
30,35,0.0,12.836999999998325,86.28128482496055
30,40,0.0,16.771999999997515,90.43154651019775
30,45,0.0,21.101000000002806,91.83195767269348
30,50,0.0,25.791000000008538,90.44021668289648
30,55,0.0,30.806000000014667,86.29457484901386
30,60,0.0,36.10800000000655,79.53000000000226
30,65,0.0,41.65599999999362,70.34058348411925
30,70,0.0,47.40899999998021,59.02925653658275
30,75,0.0,53.32199999996643,45.9119104107409
30,80,0.0,59.35099999995238,31.407745894621137
30,85,0.0,65.44999999994846,15.946886250519292
30,90,0.0,71.57199999997769,0.0

<br />
將數據匯入 SciDAVis 作圖得到的成果如下

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

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

若以 v0 = 30 m/s、theta = 30、改變 b 產生的數據如下

theta(degree), b, t(s), R(m)
30,0.0,3.060999999999774,79.5271128295263
30,0.1,2.9189999999997895,65.7683163073318
30,0.2,2.7999999999998026,55.694591159308345
30,0.30000000000000004,2.6989999999998138,48.05540251489901
30,0.4,2.6129999999998232,42.10134734984786
30,0.5,2.5369999999998316,37.33313950127536
30,0.6000000000000001,2.469999999999839,33.44825688170242
30,0.7000000000000001,2.4109999999998455,30.2339703048744
30,0.8,2.3579999999998513,27.533728914252674
30,0.9,2.3099999999998566,25.238123689064427

<br />

將數據匯入 SciDAVis 作圖得到的成果如下

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

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

## 結語

使用 VPython 可以避開很多不容易用手算的數學，例如在空氣阻力不可忽略的前提下，計算斜向抛射的飛行時間和水平射程。雖然是一次將一小段時間代入程式中運動，難免會有一點誤差，但只要將代入的時間縮短一點，計算結果就不致於偏差太多。
<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. **math.radians**: https://docs.python.org/3/library/math.html

---

###### tags:`VPython`