In [1]:
%matplotlib inline
In [2]:
from __future__ import division, print_function
from ivisual import *
from math import *
import matplotlib.pyplot as plt
In [3]:
scene=canvas(title="Sticky Collision Between Ball and Barbell")
#set these values
R=1 #m radius of barbell
vi=10 #m/s initial x-vel of incoming particle
yi=1 #m initial y of incoming particle
vibarbell=0 #m/s initial x-vel of cm of barbell
mass1=1 #kg particle
mass2=1 #kg barbell
mass3=1 #kg barbell
w=0 #initial omega of barbell
angle=90*pi/180 #initial angle of barbell in rad
time=2 #s animation time
###################################
#radius of ball
r=0.05*R
#scene
scene.autoscale=False
scene.range=4*R
#create spheres
m1=sphere(radius=r, pos=vector(-4*R,yi,0), color=color.red, make_trail=True, retain=400)
m2=sphere(radius=r, pos=vector(0,R,0), color=color.blue, make_trail=True, retain=400)
m3=sphere(radius=r, pos=vector(0,-R,0), color=color.blue, make_trail=True, retain=400)
cmbarbell=sphere(radius=r, pos=vector(-3*R,0,0), color=color.white, make_trail=True, retain=400)
cm=sphere(radius=r, pos=vector(0,0,0), color=color.cyan, make_trail=True, retain=400)
cmbarbell.visible=1
cm.visible=1
#set masses
m1.m=mass1
m2.m=mass2
m3.m=mass3
cmbarbell.m=m2.m+m3.m
cm.m=m1.m+m2.m+m3.m
#set barball variables
omega=w
theta=angle
dtheta=0
Icm=m1.m*(mag(m1.pos-cm.pos))**2+m2.m*(mag(m2.pos-cm.pos))**2+m3.m*(mag(m3.pos-cm.pos))**2
Icmbarbell=m2.m*(mag(m2.pos-cmbarbell.pos))**2+m3.m*(mag(m3.pos-cmbarbell.pos))**2
#set velocities and momenta
m1.v=vector(vi,0,0)
m1.p=m1.m*m1.v
cmbarbell.v=vector(vibarbell,0,0)
cmbarbell.p=cmbarbell.m*cmbarbell.v
cm.p=m1.p+cmbarbell.p
cm.v=cm.p/(cm.m)
#set positions of barbell spheres based on the angle of the barbell
m2.pos=cmbarbell.pos+vector(R*cos(theta),R*sin(theta),0)
m3.pos=cmbarbell.pos+vector(R*cos(theta+pi),R*sin(theta+pi),0)
cmbarbell.pos=(m2.m*m2.pos+m3.m*m3.pos)/(m2.m+m3.m)
cm.pos=(m1.m*m1.pos + m2.m*m2.pos+m3.m*m3.pos)/(m1.m+m2.m+m3.m)
#draw the barbell rod
rod=cylinder(radius=r/5, pos=m3.pos, axis=m2.pos-m3.pos, color=color.yellow)
#calculate initial angular momenta
m1.Ltrans=cross(m1.pos,m1.p)
cmbarbell.Ltrans=cross(cmbarbell.pos,cmbarbell.p)
cmbarbell.Lrot=Icmbarbell*vector(0,0,omega)
cm.Ltrans=cross(cm.pos,cm.p)
cm.Lrot=Icm*vector(0,0,omega)
Lsys=m1.Ltrans+cmbarbell.Ltrans+cmbarbell.Lrot
Ltrans=m1.Ltrans+cmbarbell.Ltrans
Lrot=cmbarbell.Lrot
Ktrans=1/2*mag(cmbarbell.p)*mag(cmbarbell.p)/cmbarbell.m+1/2*m1.m*mag(m1.v)*mag(m1.v)
Krot=1/2*mag(Lrot)*mag(Lrot)/Icm
Ktotal=Ktrans+Krot
#print initial angular momenta
print()
print("Before Collision---------------")
print("m1.Ltrans =",m1.Ltrans)
print("cmbarbell.Ltrans =",cmbarbell.Ltrans)
print("cmbarbell.Lrot =",cmbarbell.Lrot)
print("Lsys =",Lsys)
print("K=",1/2*cm.p*cm.p/cm.m)
print("Ktrans=",Ktrans)
print("Krot=",Krot)
print("Ktotal=",Ktotal)
#set clock variables
t=0
dt=0.005
#boolean for collision
beforeCollision=True
while t<time :
rate(20)
#update CM position and angle of barbell
cm.pos=cm.pos+cm.p/cm.m*dt
theta=theta+omega*dt
#update balls
if beforeCollision:
cmbarbell.pos=cmbarbell.pos+cmbarbell.p/cmbarbell.m*dt
m1.pos=m1.pos+m1.p/m1.m*dt
m2.pos=cmbarbell.pos+R*vector(cos(theta),sin(theta),0)
m3.pos=cmbarbell.pos+R*vector(cos(theta+pi),sin(theta+pi),0)
else:
m1.pos=cm.pos+d1relcm*vector(cos(theta+phi0),sin(theta+phi0),0)
m2.pos=cm.pos+d2relcm*vector(cos(theta),sin(theta),0)
m3.pos=cm.pos+d3relcm*vector(cos(theta+pi),sin(theta+pi),0)
cmbarbell.pos=(m2.m*m2.pos+m3.m*m3.pos)/(m2.m+m3.m)
#update rod
rod.pos=m3.pos
rod.axis=m2.pos-m3.pos
#update clock
t=t+dt
#see if particle collides with the barbell
if(beforeCollision):
a=m1.pos-cm.pos
b=m2.pos-cm.pos
phi=acos(dot(a,b)/mag(a)/mag(b))
phiIntersect=atan(r*mag(b)/dot(a,b))
#collision
if phi<phiIntersect and beforeCollision:
phi0 = phiIntersect #angle between vector from CM to m1 and the barbell
d1relcm=mag(cm.pos-m1.pos) #distance between cm and m1
d2relcm=mag(cm.pos-m2.pos) #distance between cm and m2
d3relcm=mag(cm.pos-m3.pos) #distance between cm and m3
Ltrans=cross(cm.pos, cm.p)
Lrot=Lsys-Ltrans
Icm=m1.m*(mag(m1.pos-cm.pos))**2+m2.m*(mag(m2.pos-cm.pos))**2+m3.m*(mag(m3.pos-cm.pos))**2
omega=Lrot.z/Icm
Ktrans=1/2*mag(cm.p)*mag(cm.p)/cm.m
Krot=1/2*mag(Lrot)*mag(Lrot)/Icm
Ktotal=Ktrans+Krot
print()
print("After Collision--------------")
print("Ltrans =",Ltrans)
print("Lrot =",Lrot)
print("Lsys =",Lsys)
print("Icm = ",Icm)
print("omega about cm =",omega)
print("Ktrans=",Ktrans)
print("Krot=",Krot)
print("Ktotal=",Ktotal)
beforeCollision=False
In [ ]: