<font color=#008000>

>作者: 吳姿瑩
>更新:2022.12.29
</font>[color=#008000]

# Lesson15: 彈簧運動（簡諧運動）
###### tags: `運動科學模擬` `多元選修` `vpython` `簡諧運動` `彈簧`

## :memo:什麼是簡諧運動Simple Harmonic Motion?
簡諧運動即是簡單又和諧的運動，位置對時間的函數圖，可用簡單的數學正弦$sin(t)$函數來表達的運動，皆是簡諧運動，簡稱為S.H.M.其中位置、速度、加速度的時間函數可表示為

- $x(t)=A\times sin(\omega t) ，其中A表示振幅，\omega 表示角速率（轉動的快慢）$
- $v(t)=A\omega \times cos(\omega t)$
- $a(t)=-A\omega^2 \times sin(\omega t)$

其函數圖如下：
<a href="https://imgur.com/La4KJfz"><img src="https://i.imgur.com/La4KJfz.png?2" title="source: imgur.com" /></a>

其中彈簧運動亦是簡諧運動其中一種，彈簧運動符合虎克定律$F=-kx$，其中$k$為彈力常數，$x$為形變量($x = 總長度－原長$)，彈簧運動如下：
<a href="https://imgur.com/ych64gG"><img src="https://i.imgur.com/ych64gG.gif" title="source: imgur.com" /></a>

<a href="https://imgur.com/McZKDXB"><img src="https://i.imgur.com/McZKDXB.gif" title="source: imgur.com" /></a>

## :memo:簡諧運動物理觀念
我們考慮較簡單的水平光滑面的彈簧運動，木塊受到３個力，即重力、正向力、彈力，其中合力即為彈力。將木塊向右拉一段距離振幅$R$ 且靜止釋放，所受合力為
$F=ma=-kx$
根據數學、簡諧運動的概念可解出
- $x(t)=R\times sin(\omega t) ，其中A表示振幅，\omega 表示角速率（轉動的快慢）$
- $v(t)=R\omega \times cos(\omega t)$
- $a(t)=-R\omega^2 \times sin(\omega t)$

其中角速率$\omega =\sqrt\frac{k}{m}$
根據週期$T =\frac{2\pi}{\omega} = 2\pi \sqrt\frac{m}{k}$

## :memo:設置一個物體box彈簧運動
### 1.參數設定
|**意涵** |**設定**|
|-----|--------|
|木塊質量 1kg|m = 4 |
|木塊邊長 0.8m |size = 0.8 |
|振幅 5 m|R = 5|
|彈性力常數 k |k = 1|
|彈簧原長 L0 |L0 = R + size|
|木塊回到初位置的次數|i = 0|
|時間間隔|dt = 0.01|
|地板厚度 0.01 m| thin = 0.01|

並在執行程式時，先輸出木塊簡諧運動的週期理論值

print("週期理論值為 {:.6f} s".format(2*pi*sqrt(m/k)))

其中==f==代表可結合==字串==輸出一同輸出，亦即結合運算式陳述(expression statements)與print()函式的運用，運用字串的 str.format()寫法可給予更詳細的輸出格式指令。
在==f前方加上6==表示，輸出的數字以小數點後6位來呈現。

### 2.畫面設定
- 有關畫布==canvas==的設定，請設定title="Simple Harmonic Motion", caption='spring motion'的畫布

畫面的 ***正中央為起點(0, 0, 0)*** 的位置

- 利用==box設定一個地板==，為了讓木塊可以平放在地板上(木塊的起始高度位置為０)，調整地板高度下降一個木塊高及地板的一半高
- 利用==box設定一個牆壁==，牆壁的起始點在最左方，故須退縮一個彈簧長度，讓彈簧可連接至牆壁
- 利用==box設定一個木塊==，木塊的起始位置在畫面做右方，故木塊中央位置拉至一個振幅R 加上木塊寬度之一半

其中這次使用到的顏色不再是單一顏色，使用vpython中的材質顏色textures。

# 畫面設定，視窗、地板、木塊
scene = canvas(title="Simple Harmonic Motion", caption='spring motion', 
               width=640, height=400, x=0, y=0,center=vec(0, 0.3, 0), 
               background=vec(204/255,204/255,255/255))
floor = box(pos=vec(0, -(size+thin)/2, 0), 
            size=vec(2*L0, thin, R), texture=textures.wood)
wall = box(pos=vec(-L0, 0, 0), 
           size=vec(0.1, size, R), texture=textures.metal)
block = box(pos=vec(R+size/2, 0, 0), 
            size=vec(size, size, size), texture=textures.rough, 
            v=vec(0, 0, 0))

- 設置==彈簧資訊，helix表螺線館之意==，其中螺線管的位置由左端開始，亦即牆面的位置(-L0, 0, 0)，radius半徑即是螺線管旋轉半徑大小，thickness厚度即是螺線管粗細
- 特別注意彈簧需連接至木塊上，其中==axis軸的資訊即是連接至何方==，而axis指向某物體（木塊）時，是連接至木塊中心，故將彈簧指向axis長度設置為，==$木塊位置 - 木塊一半寬 - 彈簧左端位置$==

#彈簧相關資訊
spring = helix(pos=vec(-L0, 0, 0), radius=0.2*size, 
               thickness=0.05*size, color=color.green)
spring.axis = block.pos - vec(size/2, 0, 0) - spring.pos

- 根據前一堂課的箭頭練習，==箭頭視覺化，將速度、加速度大小畫在木塊正上方==，故箭頭的位置即是木塊位置再加上高度為size的向量大小，讓箭頭隨著木塊一同移動，黑色表速度大小、藍色表加速度大小。

# 產生表示速度、加速度的箭頭
arrow_v = arrow(pos=block.pos + vec(0, size, 0), 
                axis=vec(0, 0, 0), shaftwidth=0.3*size, color=color.black)
arrow_a = arrow(pos=block.pos + vec(0, 2*size, 0), 
                axis=vec(0, 0, 0), shaftwidth=0.3*size, color=color.blue)

- 繪圖部分總共畫了4張圖，分別為==x-t圖、v-t圖、a-t圖、以及將位置、速度、加速度一同畫在同一張圖中==做比較

# 繪圖部分
gd = graph(title="x-t plot", width=600, height=450, x=0, y=850, 
           xtitle="<i>t</i>(s)", ytitle="<i>x</i>(m)")
gd_v = graph(title="v-t plot", width=600, height=450, x=0, y=1300, 
             xtitle="<i>t</i>(s)", ytitle="<i>v</i>(m/s)")
gd_a = graph(title="a-t plot", width=600, height=450, x=0, y=1750, 
             xtitle="<i>t</i>(s)", ytitle="<i>a</i>(m/s<sup>2</sup>)")
gd_tot = graph(title="plot", width=600, height=450, x=0, y=2200, 
               xtitle="<i>t</i>(s)",  
               ytitle="red: <i>x</i>(m), black: <i>v</i>(m/s), blue: <i>a</i>(m/s<sup>2</sup>)")

xt = gcurve(graph=gd, color=color.red)
vt = gcurve(graph=gd_v, color=color.black)
at = gcurve(graph=gd_a, color=color.blue)

xt_tot = gcurve(graph=gd_tot, color=color.red, label="red x-t plot")
vt_tot = gcurve(graph=gd_tot, color=color.black, label="black v-t plot")
at_tot = gcurve(graph=gd_tot, color=color.blue, label="blue a-t plot")

### 3.物體的運動
- 前一堂課使用了click以後才開始執行程式運動的部分

#等待click才會開始執行
T = text(text='Please click to run! ',pos=vec(-3, 5, 0), align='left', height = 0.5, color=color.black)
scene.waitfor('click')

- 利用==變數i來計算木塊歷經多少週期==，當達到3個週期的時候則結束迴圈
- 在以下程式第29行，將計算木塊完成的次數，並==印出完成第i次週期的時間==，可與週期理論值做比較
- 木塊的運動須考量虎克定律($F=-kx$)，在進行加速度、速度、位置的計算，可由以下程式第5~7行，在計算虎克定律所產生的彈力大小，==spring.axis為彈簧的總長向量，扣除彈簧原長即是彈簧伸長量==，特別注意當彈簧比原長長(spring.axis - vec(L0, 0, 0)>0)時，所受的力為負，即是回復力，反之彈簧長度比原長短時，彈力向右。
- 最後更新箭頭的位置、大小資訊，==箭頭位置隨著木塊的位置移動(以下程式第15~16行)== ，==箭頭大小分別為速度、加速度大小(以下程式第17~18行)==
- 更新畫圖資訊，將木塊位置、速度、加速度資訊畫出，==特別注意我們習慣畫彈簧伸長量的圖來表示簡諧運動，故在木塊位置顯示中，我們需要扣除木塊的一半寬==，才是真正彈簧伸長量對時間的圖。

while i < 3:
    rate(1000)

    # 計算彈簧長度、伸長量、回復力
    spring.axis = block.pos - spring.pos - vec(0.5*size, 0, 0)
    # 虎克定律 F = -kx
    F = -k * (spring.axis - vec(L0, 0, 0))

    # 計算木塊加速度, 更新速度、位置
    block.a = F/m
    block.v += block.a*dt
    block.pos += block.v*dt

    # 更新代表速度、加速度的箭頭位置、方向、長度
    arrow_v.pos = block.pos + vec(0, size, 0)
    arrow_a.pos = block.pos + vec(0, 2*size, 0)
    arrow_v.axis = block.v
    arrow_a.axis = block.a

    # 畫出 x-t, v-t, a-t 圖
    xt.plot(pos=(t, block.pos.x - size/2))
    vt.plot(pos=(t, block.v.x))
    at.plot(pos=(t, block.a.x))
    xt_tot.plot(pos=(t, block.pos.x - size/2))
    vt_tot.plot(pos=(t, block.v.x))
    at_tot.plot(pos=(t, block.a.x))


    # 判斷木塊是否回到出發點, 計算回到出發點的次數
    if block.pos.x > R + size/2 and block.v.x > 0:
        i += 1
        print("第 {:d} 次回到右端點，經過時間為 {:6f} s".format(i, t))

    # 更新時間
    t += dt

:::spoiler **簡諧運動(一組彈簧組) 完整程式碼**

"""
Auther   :吳姿瑩
Date     :2022/12/29
version  :python3.10.6 vpython7
chapter  :spring彈簧
"""

from vpython import *

"""
 1. 參數設定, 設定變數及初始值
"""
thin = 0.1          # 地板厚度
m = 4               # 木塊質量 4 kg
size = 0.8          # 木塊邊長 0.8 m
R = 5               # 振幅 5 m
k = 1               # 彈性常數 1 N/m
L0 = R + size       # 彈簧原長
i = 0               # 木塊回到初位置的次數
t = 0               # 時間
dt = 0.001          # 時間間隔
print("週期理論值為 {:.6f} s".format(2*pi*sqrt(m/k)))



"""
 2. 畫面設定
"""
# 畫面設定，視窗、地板、木塊
scene = canvas(title="Simple Harmonic Motion", caption='spring motion', width=640, height=400, x=0, y=0,center=vec(0, 0.3, 0), background=vec(204/255,204/255,255/255))
floor = box(pos=vec(0, -(size+thin)/2, 0), size=vec(2*L0, thin, R), texture=textures.wood)
wall = box(pos=vec(-L0, 0, 0), size=vec(0.1, size, R), texture=textures.metal)
block = box(pos=vec(R+size/2, 0, 0), size=vec(size, size, size), texture=textures.rough, v=vec(0, 0, 0))

#彈簧相關資訊
spring = helix(pos=vec(-L0, 0, 0), radius=0.2*size, thickness=0.05*size, color=color.green)
spring.axis = block.pos - vec(size/2, 0, 0) - spring.pos

# 產生表示速度、加速度的箭頭
arrow_v = arrow(pos=block.pos + vec(0, size, 0), axis=vec(0, 0, 0), shaftwidth=0.3*size, color=color.black)
arrow_a = arrow(pos=block.pos + vec(0, 2*size, 0), axis=vec(0, 0, 0), shaftwidth=0.3*size, color=color.blue)

# 繪圖部分
gd = graph(title="x-t plot", width=600, height=450, x=0, y=850, xtitle="<i>t</i>(s)", ytitle="<i>x</i>(m)")
gd_v = graph(title="v-t plot", width=600, height=450, x=0, y=1300, xtitle="<i>t</i>(s)", ytitle="<i>v</i>(m/s)")
gd_a = graph(title="a-t plot", width=600, height=450, x=0, y=1750, xtitle="<i>t</i>(s)", ytitle="<i>a</i>(m/s<sup>2</sup>)")
gd_tot = graph(title="plot", width=600, height=450, x=0, y=2200, xtitle="<i>t</i>(s)",  ytitle="red: <i>x</i>(m), black: <i>v</i>(m/s), blue: <i>a</i>(m/s<sup>2</sup>)")

xt = gcurve(graph=gd, color=color.red)
vt = gcurve(graph=gd_v, color=color.black)
at = gcurve(graph=gd_a, color=color.blue)

xt_tot = gcurve(graph=gd_tot, color=color.red, label="red x-t plot")
vt_tot = gcurve(graph=gd_tot, color=color.black, label="black v-t plot")
at_tot = gcurve(graph=gd_tot, color=color.blue, label="blue a-t plot")


"""
 3. 物體運動部分, 重複3個週期
"""
#等待click才會開始執行
T = text(text='Please click to run! ',pos=vec(-3, 5, 0), align='left', height = 0.5, color=color.black)
scene.waitfor('click')

while i < 3:
    rate(1000)

    # 計算彈簧長度、伸長量、回復力
    spring.axis = block.pos - spring.pos - vec(0.5*size, 0, 0)
    # 虎克定律 F = -kx
    F = -k * (spring.axis - vec(L0, 0, 0))

    # 計算木塊加速度, 更新速度、位置
    block.a = F/m
    block.v += block.a*dt
    block.pos += block.v*dt

    # 更新代表速度、加速度的箭頭位置、方向、長度
    arrow_v.pos = block.pos + vec(0, size, 0)
    arrow_a.pos = block.pos + vec(0, 2*size, 0)
    arrow_v.axis = block.v
    arrow_a.axis = block.a

    # 畫出 x-t, v-t, a-t 圖
    xt.plot(pos=(t, block.pos.x - size/2))
    vt.plot(pos=(t, block.v.x))
    at.plot(pos=(t, block.a.x))
    xt_tot.plot(pos=(t, block.pos.x - size/2))
    vt_tot.plot(pos=(t, block.v.x))
    at_tot.plot(pos=(t, block.a.x))


    # 判斷木塊是否回到出發點, 計算回到出發點的次數
    if block.pos.x > R + size/2 and block.v.x > 0:
        i += 1
        print("第 {:d} 次回到右端點，經過時間為 {:f} s".format(i, t))

    # 更新時間
    t += dt

:::
:::spoiler **學生活動**
1. 試求m=1 kg,k=1 N/m的週期理論值？
2. 試求程式中完成一週期的時間為何？
3. 試問，當木塊在左右端點時，速度、加速度何者最大？
4. 試問，當木塊在平衡點時(彈簧原長處)，速度、加速度何者最大？
:::
:::spoiler **學生作業**
1. 1. Homework:設置==兩個不同顏色的彈簧分別繫在4kg 木塊上，彈簧力常數分別為k=1, 0.25 N/m==，並畫出 ==$v$和$a$，速度、加速度箭頭視覺化==
Requirements要求:
(1)需畫出木塊不重疊的移動
(2)需畫出兩個彈簧搭配顏色的(1)x-t圖 (2)v-t圖 (3)a-t圖
(3)需畫出各自木塊的速度、加速度視覺化箭頭
2. 請將==程式code以及瀏覽器結果截圖== 至google classroom繳交作業。
3. 完成今日協作平台更新
4. 將code上傳至雲端硬碟中(Lesson15 vpython彈簧運動（簡諧運動）)
:::

## :memo:參考網址
1. python 輸入和輸出: https://docs.python.org/zh-tw/3/tutorial/inputoutput.html
2. vpython 的textures材質顏色：https://www.glowscript.org/docs/VPythonDocs/textures.html
3. vpython textures 列表：https://www.glowscript.org/#/user/GlowScriptDemos/folder/Examples/program/Textures-VPython
4. vpython 彈簧螺線管使用：https://www.glowscript.org/docs/VPythonDocs/helix.html
5. vpython 箭頭arrow：https://www.glowscript.org/docs/VPythonDocs/arrow.html