Charged_particle_collision.ipynb

Kernel: Python 2 (SageMath)

from vpython import *
from math import sqrt

# Written by Ruth Chabay, licensed under Creative Commons 4.0.
# All uses permitted, but you must not claim that you wrote it, and
# you must include this license information in any copies you make.
# For details see http://creativecommons.org/licenses/by/4.0

## impact parameter
##b = 0
b = 15e-15  
##b = 80e-15

# alpha particle charge/e, mass/mp
projectileproperties = (2,4) 
##projectileproperties = (-1,(9e-31/1.7e-27))    ## electron

# oxygen nucleus charge/3, mass/mp
targetproperties = (8,16) 
##targetproperties = (2,4) # alpha particle
##targetproperties = (1,(9e-31/1.7e-27))    ## positron

## make radii bigger for visual impact
rpscale = 2 
parroVisible = 1

scene.width = 1000
scene.height = 650
scene.x = scene.y = 0
scene.background = color.white
scene.fov = 0.01
range0 = 200e-15
scene.range = range0
xstart = 0.95*range0

kcoul = 9e9
qe = 1.6e-19
mproton = 1.7e-27
rproton = 1.3e-15*rpscale
alpha = sphere(pos=vector(-xstart,b,0), color=color.red, make_trail=True, interval=40)
target = sphere(pos=vector(0,0,0), color=color.blue, make_trail=True, interval=40)
alpha.mass = projectileproperties[1]*mproton
alpha.radius = (alpha.mass/mproton)**(1./3.)*rproton
alpha.q = projectileproperties[0]*qe
ke = 1e6*qe
alpha.p = vector(sqrt(2.*alpha.mass*ke),0,0)
target.mass = targetproperties[1]*mproton
target.radius = (target.mass/mproton)**(1./3.)*rproton
target.q = targetproperties[0]*qe
target.p = vector(0,0,0)
dt = (5.*xstart/(mag(alpha.p)/alpha.mass)/1e5)*10
ptot = alpha.p+target.p
vcm = ptot/(alpha.mass+target.mass)

pscale = 40e-15/4e-20
paarro = arrow(pos=alpha.pos, axis=alpha.p*pscale, color=color.cyan, shaftwidth = 0.5*alpha.radius, fixedwidth=1, visible=parroVisible)
ptarro = arrow(pos=target.pos, axis=target.p*pscale, color=color.magenta, shaftwidth = 0.5*alpha.radius, fixedwidth=1, visible=parroVisible)

##alpha.p = alpha.p-alpha.mass*vcm
##target.p = target.p-target.mass*vcm


while True:
    rate(500)
    r12 = alpha.pos-target.pos
    F = (kcoul*alpha.q*target.q/mag(r12)**2)*norm(r12)
    alpha.p = alpha.p + F*dt
    target.p = target.p - F*dt
    alpha.pos = alpha.pos + (alpha.p/alpha.mass)*dt
    target.pos = target.pos + (target.p/target.mass)*dt
    paarro.pos = alpha.pos
    paarro.axis = alpha.p*pscale
    ptarro.pos = target.pos
    ptarro.axis = target.p*pscale