#copyright (c) Aaron Titus, 2004
#licensed under the gpl license

#http://linus.highpoint.edu/~atitus/mandi/vpython

from visual import *

print "Perfectly Inelastic collision with free barbell."

print "\nClick to start the animation."

R=1.0
r=0.05*R

scene.autoscale=false
scene.range=3*R
scene.width=800
scene.height=400

m1=sphere(radius=r, pos=vector(-2*R,R,0), color=color.red)
m2=sphere(radius=r, pos=vector(0,R,0), color=color.blue)
m3=sphere(radius=r, pos=vector(0,-R,0), color=color.blue)
cm=sphere(radius=r, pos=vector(0,0,0), color=color.white)
cm.visible=0
rod=cylinder(radius=r/5., pos=m3.pos, axis=m2.pos-m3.pos, color=color.yellow)

m1.m=1.
m2.m=1.
m3.m=1.
cm.m=m2.m+m3.m

cm.pos=(m2.m*m2.pos+m3.m*m3.pos)/cm.m

m1.p=vector(1.,0,0)
cm.p=vector(0,0,0)
m2.p=vector(0,0,0)
m3.p=vector(0,0,0)

omega=0
I=m2.m*R**2+m3.m*R**2

m1.L=cross(m1.pos,m1.p)
cm.L=cross(cm.pos,cm.p)
Lrot=vector(0,0,I*omega)

t=0
dt=0.01
dtheta=0
theta2=atan2(m2.pos.y,m2.pos.x)
theta3=atan2(m3.pos.y,m3.pos.x)

scale=1
Lscale=1

showVectors=true

if showVectors:
	pVector1=arrow(pos=m1.pos, axis=m1.p*scale, color=m1.color)
	pVector2=arrow(pos=m2.pos, axis=m2.p*scale, color=m2.color)
	pVector3=arrow(pos=m3.pos, axis=m3.p*scale, color=m3.color)
	pVectorCM=arrow(pos=cm.pos, axis=cm.p*scale, color=cm.color)
	LVector1=arrow(pos=m1.pos, axis=m1.L*Lscale, color=color.green)
	LVectorCM=arrow(pos=cm.pos, axis=cm.L*Lscale, color=color.green)
	LVectorRot=arrow(pos=cm.pos, axis=Lrot*scale, color=color.green)

beforeCollision=true

scene.mouse.getclick()

while 1:
	rate(30)

	if showVectors:
		pVector1.pos=m1.pos
		pVector1.axis=m1.p*scale
		pVectorCM.pos=cm.pos
		pVectorCM.axis=cm.p*scale
		LVector1.pos=m1.pos
		LVector1.axis=m1.L*Lscale
		LVectorCM.pos=cm.pos
		LVectorCM.axis=cm.L*Lscale
		LVectorRot.pos=cm.pos
		LVectorRot.axis=Lrot*Lscale
	
	t=t+dt
	
	dtheta=omega*dt
	theta2=theta2+dtheta
	theta3=theta3+dtheta
	cm.pos=cm.pos+cm.p/cm.m*dt
	if beforeCollision:
		m1.pos=m1.pos+m1.p/m1.m*dt
	else:
		theta1=theta2+r/R
		m1.pos=cm.pos+vector(R*cos(theta1),R*sin(theta1),0)
	m2.pos=cm.pos+vector(R*cos(theta2),R*sin(theta2),0)
	m3.pos=cm.pos+vector(R*cos(theta3),R*sin(theta3),0)
	rod.pos=m3.pos
	rod.axis=m2.pos-m3.pos

	m1.L=cross(m1.pos,m1.p)
	cm.L=cross(cm.pos,cm.p)
	Lrot=vector(0,0,omega*I)
			
	if (mag(m1.pos-m2.pos)<r or mag(m1.pos-m3.pos)<r) and beforeCollision:
		#initial angular momentum
		m1.Li=m1.L
		
		print m1.L

		#final momentum after collision
		cm.p=m1.p

		#final angular momentum
		cm.L=cross(cm.pos,cm.p)

		#angular velocity
		Lrot=(m1.Li-cm.L)
		omega=1/I*Lrot.z
		
		pVector1.visible=0
		LVector1.visible=0
		
		beforeCollision=false