YouTube
Simulation [gnuplot]
Model of a Bouncing Ball
Equations of motion \begin{eqnarray} \label{eq:eom2} \left\{ \begin{array}{l} \frac{dX}{dt}&=&V_X\\ \frac{dY}{dt}&=&V_Y\\ \frac{d{V}_{X}}{dt}&=&-g\sin\alpha\\ \frac{d{V}_{Y}}{dt}&=&-g\cos\alpha \end{array} \right. \end{eqnarray} Inelastic collision \begin{eqnarray} V_{X}^{\rm{out}}&=-\color{red}{e}V_{X}^{\rm{in}}\\ V_{Y}^{\rm{out}}&=-\color{red}{e}V_{Y}^{\rm{in}} \end{eqnarray} Rotation matrix \begin{equation} R(\alpha) = \left[ \begin{array}{cc} \cos\alpha & -\sin\alpha \\ \sin\alpha & \cos\alpha \end{array} \right] \end{equation} Rotation using the matrix \begin{eqnarray} \left[\begin{array}{c} x \\ y\end{array}\right] &=& R(\alpha) \left[\begin{array}{c}X \\ Y \end{array}\right] &=& \left[ \begin{array}{c} X\cos\alpha - Y\sin\alpha \\ X\sin\alpha + Y\cos\alpha \end{array} \right]\\ \left[\begin{array}{c} v_x \\ v_y\end{array}\right] &=& R(\alpha) \left[\begin{array}{c}V_X \\ V_Y \end{array}\right] &=& \left[ \begin{array}{c} V_X\cos\alpha - V_Y\sin\alpha \\ V_X\sin\alpha + V_Y\cos\alpha \end{array} \right]\\ \end{eqnarray}
Case 1: $\alpha=20^\circ$
PLT file ①
# Setting ----------------------------------------
reset
set term gif animate delay 5 size 1280, 720
set output "BallBox.gif"
set nokey
unset grid
A = -250 # A<=x<=B
B = 250
C = -60 # C<=y<=D
D = 180
set xr[A:B]
set yr[C:D]
set xl "{/TimesNewRoman:Italic=24 x}"
set yl "{/TimesNewRoman:Italic=24 y}"
set tics front font 'Times New Roman,20'
set xtics 50
set ytics 50
set size ratio -1
# Parameter ----------------------------------------
m = 1.5 # mass [kg]
g = 9.8 # gravitational acceleration [m/s2]
r = 4.0 # Radius of the ball
V = 40 # velocity [m/s]
a = 20*pi/180 # Slope angle [rad] (-90<a<90)
b = 70*pi/180 # Projection angle [rad]
e = 0.8 # The coefficient of restitution [-]
ep = 2 # Epsilon
dt = 0.001 # Step [s]
dh = dt/6.0
cut = 200 # Decimation
dis = 800 # Start to disappear
N = 5 # The number of balls
cnt = 0 # The number of balls being framed out
array X[N] # Position of balls
array Y[N]
array VX[N] # Velocity of balls
array VY[N]
array x[N] # Position of balls after a rotation
array y[N]
array vx[N] # Velocity of balls after a rotation
array vy[N]
# Color of balls
array color[N] = ["royalblue", "red", "orange", "green", "black"]
# Border line
Lu = 130
Ld = 0
Lr = 150
Ll = -150
# Goal
if(a>0){
GOALx = Ll + r + ep
} else { if(a<0) {
GOALx = Lr - r - ep
} else { if(a==0) {
# Not define GOALx
}}}
GOALy = Ld + r + ep
#Rotation matrix R(a)
# R[1]=R11, R[2]=R12, R[3]=R21, R[4]=R22
array R[4] = [cos(a), -sin(a), sin(a), cos(a)]
# Functions ----------------------------------------
# Equations of Motion
f1(x, y, vx, vy) = vx # dx/dt
f2(x, y, vx, vy) = vy # dy/dt
f3(x, y, vx, vy) = -g*sin(a) # dvx/dt
f4(x, y, vx, vy) = -g*cos(a) # dvy/dt
# 4th order Runge-Kutta (Define rk_i(x, y, vx, vy))
do for[i=1:4]{
rki = "rk"
fi = "f".sprintf("%d", i)
rki = rki.sprintf("%d(x, y, vx, vy) = (\
k1 = %s(x, y, vx, vy),\
k2 = %s(x + dt*k1/2., y + dt*k1/2., vx + dt*k1/2., vy + dt*k1/2.),\
k3 = %s(x + dt*k2/2., y + dt*k2/2., vx + dt*k2/2., vy + dt*k2/2.),\
k4 = %s(x + dt*k3, y + dt*k3, vx + dt*k3, vy + dt*k3),\
dh * (k1 + 2*k2 + 2*k3 + k4))", i, fi, fi, fi, fi)
eval rki
}
# X->x, Y->y (x=R(a)X)
x(a, b) = R[1]*a + R[2]*b
y(a, b) = R[3]*a + R[4]*b
# x->X, y->Y (X=R(-a)x)
X(a, b) = R[4]*a - R[2]*b
Y(a, b) = -R[3]*a + R[1]*b
# Time
Time(t) = sprintf("{/TimesNewRoman:Italic t} = %.1f", t)
# Plot ----------------------------------------
# Initial Value
t = 0
b_inc = 36 # Increment for b
do for[j=1:N]{
X[j] = X(0, 60)
Y[j] = Y(0, 60)
VX[j] = V*cos(b + (b_inc*(j-int(N/2))) * pi/180)
VY[j] = V*sin(b + (b_inc*(j-int(N/2))) * pi/180)
# Vector rotation (Rotation matrix R(a))
x[j] = x(X[j] , Y[j] )
y[j] = y(X[j] , Y[j] )
vx[j] = x(VX[j], VY[j])
vy[j] = y(VX[j], VY[j])
}
# Draw initiate state for 70 steps
do for [i=1:70] {
# Time
set label 1 Time(t) left at graph 0.05, 0.93 font 'TimesNewRoman, 25' front
# Balls
do for[j=1:N]{
set obj j circ at x[j], y[j] size r fc rgb color[j] fs solid noborder
}
# Draw ground and balls
if(a != 0){
plot tan(a)*(x-x(0, Ld)) + y(0, Ld) with filledcurve x1 lc rgb 'gray50', \
tan(a)*(x-x(0, Lu)) + y(0, Lu) with filledcurve x2 lc rgb 'gray50', \
-1/tan(a)*(x-x(Ll, 0)) + y(Ll, 0) with filledcurve y1 lc rgb 'gray50', \
-1/tan(a)*(x-x(Lr, 0)) + y(Lr, 0) with filledcurve y2 lc rgb 'gray50'
} else {
plot Ld with filledcurve x1 lc rgb 'gray50', \
Lu with filledcurve x2 lc rgb 'gray50', \
((x<Ll)||(x>Lr) ? Lu : 1/0) with filledcurve y1=Ld lc rgb 'gray50'
}
}
# Update until all of balls are framed out
i = 0 # The number of loops
while(1){
i = i +1
t = t + dt
set label 1 Time(t)
do for[j=1:N]{
# 4th order Runge-Kutta
temp_X = X[j] + rk1(X[j], Y[j], VX[j], VY[j])
temp_Y = Y[j] + rk2(X[j], Y[j], VX[j], VY[j])
temp_VX = VX[j] + rk3(X[j], Y[j], VX[j], VY[j])
temp_VY = VY[j] + rk4(X[j], Y[j], VX[j], VY[j])
X[j]=temp_X; Y[j]=temp_Y; VX[j]=temp_VX; VY[j]=temp_VY
# Judge whether balls bounce or not
if(Y[j] > Lu-r){
Y[j] = Lu-r
VY[j] = -e*VY[j]
}
if(Y[j] < Ld+r){
Y[j] = Ld+r
VY[j] = -e*VY[j]
}
if(X[j] > Lr-r){
X[j] = Lr-r
VX[j] = -e*VX[j]
}
if(X[j] < Ll+r){
X[j] = Ll+r
VX[j] = -e*VX[j]
}
# Vector rotation (Rotation matrix R(a))
x[j] = x(X[j] , Y[j] )
y[j] = y(X[j] , Y[j] )
vx[j] = x(VX[j], VY[j])
vy[j] = y(VX[j], VY[j])
# Update objects
set obj N*i+j circ at x[j], y[j] size r fc rgb color[j] fs solid noborder
# Make old objects trajectory of the ball
set obj N*(i-1)+j at x[j], y[j] size 0.1
# Start to disappear
if(i>=dis){
unset object 1+N*(i-dis)+j
}
}
# Decimate and plot
if(i%cut==0){
replot
}
# Count balls meeting a condition
do for[j=1:N]{
if(a > 0){
if(X[j]<GOALx && Y[j]<GOALy){
cnt = cnt + 1 # whether a ball stands still or not
}
} else { if(a < 0){
if(X[j]>GOALx && Y[j]<GOALy){
cnt = cnt + 1 # whether a ball stands still or not
}
} else {
if(Y[j]<GOALy){
cnt = cnt + 1 # whether a ball is grounded or not
}
}}
}
# Exit from the loop when all of balls meet the condition
if(cnt == N){
break
} else {
cnt =0
}
}
set out
GIF file ①
$V=40\ \mathrm{m/s},\ a=20^\circ,\ b=70^\circ,\ e=0.8,\ \varepsilon=0.1,\\L_u=1320,\ L_d=0,\ L_r=150,\ L_l=-150$
Case 2: $\alpha=45^\circ$
PLT file ②
A = -200 B = 200 C = -200 D = 200
V = 40 a = 45*pi/180 b = 50*pi/180 e = 0.8 ep = 2
Lu = 120 Ld = -120 Lr = 120 Ll = -120
GIF file ②
$V=40\ \mathrm{m/s},\ a=45^\circ,\ b=50^\circ,\ e=0.8,\ \varepsilon=2,\\L_u=120,\ L_d=-120,\ L_r=120,\ L_l=-120$
Case 3: $\alpha=0^\circ$
GIF file ④
$V=40\ \mathrm{m/s},\ a=0^\circ,\ b=85^\circ,\ e=0.8,\ \varepsilon=2,\\L_u=130,\ L_d=0,\ L_r=150,\ L_l=-150$
Extra
Case 4: $\alpha=5^\circ$
GIF file ③
$V=60\ \mathrm{m/s},\ a=5^\circ,\ b=85^\circ,\ e=0.5,\ \varepsilon=2,\\L_u=130,\ L_d=0,\ L_r=150,\ L_l=-150$
