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

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

from visual import *

scene.autoscale=0
scene.width=scene.height=800

print "Rutherford scattering"
print "----------------------------"
print ""

qe=1.6e-19	#charge of electron
NA=6.02e23 #Avogadro's number
const=9e9 #Elec force constant

KeV=1.0e7 #eV
K=KeV*qe #J
qAlpha=2.0*qe
qGold=79.0*qe
mGold=197.0/NA/1000.0
mAlpha=4.0/NA/1000.0

pAlpha=sqrt(2.0*mAlpha*K)

rf=2.0*const*qAlpha*qGold*mAlpha/pAlpha**2
ri=const*qAlpha*qGold/(0.01*K)

scene.range=1.2*ri

t=0
dt=0.001*(mAlpha*ri/pAlpha)

#create alpha particles
alphas=[]
Nalphas=50
bmax=ri/4.0
for i in arange(0,Nalphas,1):
	alphas.append(sphere(pos=vector(-ri,-bmax+i*2*bmax/(Nalphas*1.0),0), radius=0.01*ri, color=color.red))
	alphas[i].p=vector(pAlpha,0,0)

gold=sphere(pos=vector(0,0,0), radius=0.01*ri, color=color.yellow)

gold.p=vector(0,0,0)

while 1:
	rate(100)
	t=t+dt
	
	for i in range(len(alphas)):
		part=alphas[i]
		r=part.pos-gold.pos
		rmag=mag(r)
		runit=r/rmag

		Falpha=const*qAlpha*qGold/rmag**2*runit
		part.p=part.p+Falpha*dt
		part.pos=part.pos+part.p/mAlpha*dt