"""
VPython教學: 13-1.相疊的木塊完全非彈性碰撞
Ver. 1: 2017/6/16
Ver. 2: 2018/1/4 修改為 Python 3.X 版程式碼
Ver. 3: 2018/2/26 增加 v-t 圖
Ver. 4: 2019/9/8
作者: 王一哲
"""
from vpython import *
"""
1. 參數設定, 設定變數及初始值
(1) m1 = 0.2, v1 = 0.0, m2 = 0.1, v2 = 2.0
(2) m1 = 0.2, v1 = 1.0, m2 = 0.1, v2 = 0.0
(3) m1 = 0.2, v1 = 0.0, m2 = 0.1, v2 = 2.5 => b2 會飛出去
(4) m1 = 0.2, v1 = 2.0, m2 = 0.1, v2 = 0.0 => b1 到畫面右側時還沒有達到等速度
"""
d1 , h1, w1 = 1.8, 0.2, 0.2
d2 , h2, w2 = 0.2, 0.2, 0.2
m1, v1, c1 = 0.2, 0.0, color.red
m2, v2, c2 = 0.1, 2.0, color.green
xmax, xmin = 2.0, -2.0
g = 9.8
mu = 0.1
dt = 0.0005
t = 0
bx = 0
i, te = 0, -1
"""
2. 畫面設定
"""
scene = canvas(title="Two Blocks", width=800, height=300, center=vec(0, 0.4, 0), background=vec(0, 0.6, 0.6))
floor = box(pos=vec(0, -0.5*h2, 0), size=vec(xmax - xmin, 0.05, 0.8), color=color.blue)
b1 = box(pos=vec(xmin + d1/2, 0, 0), size=vec(d1, h1, w1), color = c1, v=vec(v1, 0, 0))
if(v2 >= v1): bx = xmin + 0.5*d2
else: bx = xmin + d1 - 0.5*d2
b2 = box(pos=vec(bx, h1, 0), size=vec(d2, h2, w2), color=c2, v=vec(v2, 0, 0))
gd = graph(title="<i>E</i> - <i>t</i> plot", x=0, y=300, width=600, height=450, xtitle="<i>t</i> (s)",
ytitle="red: <i>K</i><sub>1</sub>, green: <i>K</i><sub>2</sub>, blue: <i>E</i> (J)")
kt1 = gcurve(graph=gd, color=c1)
kt2 = gcurve(graph=gd, color=c2)
et = gcurve(graph=gd, color=color.blue)
gd2 = graph(title="<i>v</i> - <i>t</i> plot", x=0, y=750, width=600, height=450, xtitle="<i>t</i> (s)",
ytitle="red: <i>v</i><sub>1</sub>, green: <i>v</i><sub>2</sub> (m/s)")
vt1 = gcurve(graph=gd2, color=c1)
vt2 = gcurve(graph=gd2, color=c2)
"""
3. 物體運動部分, 重複執行直到: (1) b1 抵達邊緣時, (2) b2 抵達 b1 邊緣時
"""
while((b1.pos.x <= xmax - d1/2) and (b2.pos.x + d2/2 <= b1.pos.x + d1/2 + 0.001)):
rate(500)
if(b2.v.x > b1.v.x):
force = mu * m2 * g
b1.a = vec(force / m1, 0, 0)
b2.a = vec(-force / m2, 0, 0)
elif(b2.v.x < b1.v.x):
force = mu * m2 * g
b1.a = vec(-force / m1, 0, 0)
b2.a = vec(force / m2, 0, 0)
else:
force = 0
b1.a = vec(0, 0, 0)
b2.a = vec(0, 0, 0)
if(abs(b2.v.x - b1.v.x) < 0.0005 and i == 0):
te = t
i = 1
b1.v += b1.a * dt
b2.v += b2.a * dt
b1.pos += b1.v * dt
b2.pos += b2.v * dt
k1 = 0.5 * m1 * mag2(b1.v)
k2 = 0.5 * m2 * mag2(b2.v)
e = k1 + k2
kt1.plot(pos = (t, k1))
kt2.plot(pos = (t, k2))
et.plot(pos = (t, e))
vt1.plot(pos = (t, b1.v.x))
vt2.plot(pos = (t, b2.v.x))
t += dt
print("v1 = ", b1.v.x)
print("v2 = ", b2.v.x)
print("te = ", te)
print("end")