From Aravir |
# Import the library(s)
#########################################
thk = 0.1 #ketebalan
side = 7.0 #lebar
s2 = 2*side – thk
s3 = 2*side + thk
#wallR=tembok kanan; wallL=tembok kiri; …dst
wallL = box (pos=vector(-side, 0, 0), length=thk, height=s2, width=s3, color = color.red)
wallB = box (pos=vector(0, -side, 0), length=s3, height=thk, width=s3, color = color.blue)
wallT = box (pos=vector(0, side, 0), length=s3, height=thk, width=s3, color = color.blue)
wallBK = box(pos=vector(0, 0, -side), length=s2, height=s2, width=thk, color = (0.7,0.7,0.7))
##########################################
# Buat Bola
no_particles=37
ball_radius=1.
maxpos=side-.5*thk-ball_radius
maxv=7.0
#ball_list=[] adalah perintah membuat list (semacam array di delphi) bernama #ball_list, berfungsi untuk menampung objek bernama ball yang terbuat dari object #bernama sphere
ball=sphere(color=color.green,radius=jari)
ball.pos=maxpos*vector(uniform(-1,1),uniform(-1,1),uniform(-1,1))
ball.velocity=maxv*vector(uniform(-1,1),uniform(-1,1),uniform(-1,1))
#ball.velocity=vector(0,0,0)
#ball.mass=jari
ball_list.append(ball)
#ball_list[0].velocity=vector(14,0,0)
###########################################
# Setup Histogram Plot/buat grafik
###########################################
graphwindow=gdisplay(xtitle=’v’,ytitle=’N’,ymax=no_particles/2)
velocity_dist=ghistogram(bins=arange(0,2*maxv,maxv/5),accumulate=1,average=1)
expected_distribution=gcurve(color=color.green)
for vx in arange(0,2*maxv,maxv/20):
expected_distribution.plot(pos = (vx,.27*no_particles*exp(-vx**2/maxv**2*3/2)))
##########################################
# Time loop for moving Ball(s)
###########################################
timestep = 0.05
while (1==1):
rate(100)
######################################
# Loop over list of particles
# Move and check wall collisions
######################################
for ball in ball_list:
suhu=mag(ball.velocity)/maxv
ball.color=(suhu/2,1-suhu/2,0)
#move ball
ball.pos = ball.pos + ball.velocity*timestep
#check right wall collisions
if ball.x > maxpos:
#reflect velocity
ball.velocity.x = -ball.velocity.x
#reflect position
ball.x=2*maxpos-ball.x
#left wall
if ball.x < -maxpos:
ball.velocity.x = -ball.velocity.x
ball.x=-2*maxpos-ball.x
# roof
if ball.y > maxpos:
ball.velocity.y = -ball.velocity.y
ball.y=2*maxpos-ball.y
#floor
if ball.y < -maxpos:
ball.velocity.y = -ball.velocity.y
ball.y=-2*maxpos-ball.y
#back wall
if ball.z > maxpos:
ball.velocity.z = -ball.velocity.z
ball.z=2*maxpos-ball.z
#front wall
if ball.z < -maxpos:
ball.velocity.z = -ball.velocity.z
ball.z=-2*maxpos-ball.z
######################################
# Ball Collision Detection
######################################
#loop through all pairs
for i in range(no_particles):
for j in range(i+1,no_particles):
distance=mag(ball_list[i].pos-ball_list[j].pos)
#check collision
if distance<(ball_list[i].radius+ball_list[j].radius):
#unit vector in collision direction
direction=norm(ball_list[j].pos-ball_list[i].pos)
vi=dot(ball_list[i].velocity,direction)
vj=dot(ball_list[j].velocity,direction)
#mi=ball_list[i].radius
#mj=ball_list[j].radius
#impact velocity
exchange=vj-vi
##v2’=1/(m1+m2)*(2m1v1+(m2-m1)v2)
##v1’=1/(m1+m2)*(2m2v2+(m2-m1)v1)
#vii=((mi-mj)*vi+2.*mj*vj)/(mi+mj)
#vjj=((mi-mj)*vj+2.*mi*vi)/(mi+mj)
##exchange momentum
ball_list[i].velocity=ball_list[i].velocity + exchange*direction
ball_list[j].velocity=ball_list[j].velocity – exchange*direction
#ball_list[i].velocity=vii*direction
#ball_list[j].velocity=vjj*direction
##adjust position
#overlap=2*ball_radius-distance
overlap=ball_list[i].radius+ball_list[j].radius-distance
ball_list[i].pos=ball_list[i].pos – overlap*direction
ball_list[j].pos=ball_list[j].pos + overlap*direction
#######################################
# Plot the x velocity histogram
#######################################
v_list=[]
for ball in ball_list:
v_list.append(mag(ball.velocity))
velocity_dist.plot(data=v_list)
###########################################
From Aravir |
From Aravir |
From Aravir |