anonymous / rutherford_thomson.py
Created

Embed URL

HTTPS clone URL

SSH clone URL

You can clone with HTTPS or SSH.

Download Gist

VPython program to simulate scattering of alpha particles off gold atoms.

View rutherford_thomson.py
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91
from __future__ import division
from visual import *
from visual.graph import *
from math import *
 
size=1e-15
size2=1e-10
b=0.005*size2
start=100*size2
k=9e9
e=1.6e-19
 
alpha=sphere(pos=(-start,b,0),radius=0.1*size,color=color.red)
 
alpha.mass=4*1.67e-27
alpha.charge=+2*e
alpha.size=1e-12
alpha.energy=4e6*1.6e-19
alpha.speed=sqrt(2*alpha.energy/alpha.mass)
 
gold=sphere(pos=(0,0,0), radius=size, color=color.yellow)
 
gold.mass=197*1.67e-27
gold.charge=+79*e
 
shell=sphere(pos=gold.pos, radius=size2, color=color.green)
shell.mass=gold.mass/1837
shell.charge=-gold.charge
 
alpha.momentum=alpha.mass*alpha.speed*vector(1,0,0)
gold.momentum=gold.mass*vector(0,0,0)
 
alpha.trail=curve(color=alpha.color)
alpha.trail.append(pos=alpha.pos)
gold.trail=curve(color=gold.color)
gold.trail.append(pos=gold.pos)
 
r=alpha.pos-gold.pos
rhat=r/r.mag
r_min=r.mag
 
init_r=alpha.pos.mag
 
print(r,rhat,r.mag2)
 
dt=0.01*r.mag/alpha.speed
t=0
 
while alpha.pos.mag<2*init_r:
if r.mag>size:
F_alpha=k*alpha.charge*gold.charge*rhat/r.mag2
F_gold=-F_alpha
else:
F_alpha=k*alpha.charge*gold.charge*rhat*(r.mag/size)/(size*size)
F_gold=-F_alpha
 
if r.mag>size2:
F_alpha=F_alpha+k*alpha.charge*shell.charge*rhat/r.mag2
F_gold=-F_alpha
check=F_alpha+F_gold
#print(F_alpha.y)
 
alpha.momentum=alpha.momentum+F_alpha*dt
gold.momentum=gold.momentum+F_gold*dt
 
alpha.pos=alpha.pos+(alpha.momentum/alpha.mass)*dt
alpha.trail.append(pos=alpha.pos)
gold.pos=gold.pos+(gold.momentum/gold.mass)*dt
shell.pos=gold.pos
gold.trail.append(pos=gold.pos)
 
r=alpha.pos-gold.pos
rhat=r/r.mag
 
alpha.speed=alpha.momentum.mag/alpha.mass
dt=0.01*r.mag/alpha.speed
if r.mag<r_min:
r_min=r.mag
t=t+dt
 
if alpha.pos.x>0:
angle=degrees(atan(abs((alpha.pos.y-b)/alpha.pos.x)))
else:
angle=180-degrees(atan(abs((alpha.pos.y-b)/alpha.pos.x)))
 
print(b,alpha.pos.y/b,angle)
print(r.mag/size,r_min/size,r_min/size2)
 
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Something went wrong with that request. Please try again.