from vpython import *

#GlowScript 3.1 VPython

# Rutherford Scattering Exercise

# Define Constants

q_e = 1.6e-19 # Charge on a proton [Coulombs]

m_p = 1.7e-27 # Mass of proton [kg]

oofpez = 9e9 # oofpez = One Over Four Pi Epsilon Zero [N m**2 / C**2]

K0 = 5.5e6 * 1.6e-19 # Energy equivalent of 5.5MeV in J

b = 5e-15 # Impact Parameter [m]

# Model the Au nucleus as a yellow sphere of radius 4 fm, located at the origin

Au = sphere(pos=vector(0,0,0), radius=4e-15, color=color.yellow, make_trail=True)

Au.opacity = 0.7 # Make Au a little opaque so its trail can be seen eaiser.

Au.m = (79+118) * m_p # Mass of Gold nucleus, 79 protons and 118 neutrons [kg]

Au.q = 79 * q_e # Charge of Gold nucleus [Coulombs]

# Define momentum of Au. Since initially at rest, velocity = <0,0,0>

Au.p = Au.m*vector(0,0,0)

# Model the Alpha particle as a magenta sphere of radius 1 fm

Alpha = sphere(radius=1e-15, color=color.magenta, make_trail=True)

# Set the initial position of the Alpha particle to be a large distance to the left

# of the Au nucleus and above the x axis by an amount equal to the impact parameter.

Alpha.pos = vector(-100e-15,b,0)

Alpha.m = (2+2) * m_p # Mass of Alpha particle, 2 protons and 2 neutrons [kg]

Alpha.q = 2 * q_e # Charge of Alpha particle [Coulombs]

# Define initial momentum of Alpha particle using its initial kinetic energy.

Alpha.p = vector(sqrt(2*Alpha.m*K0),0,0) # Initially moving in +x direction.

Alpha.p1 = Alpha.p # Save the initial momentum for later

# Compute the initial separation of the Alpha and Au and store it in variable r0

r = Alpha.pos - Au.pos

r0 = mag(r) # Save the initial separation

# Set up the graph displays

gdx = gdisplay(width=600, height=200, title='x Momentum vs Time', background=vector(0.8, 0.8, 0.8))

p_Au_x_graph = gcurve(color=color.yellow, label="Au")

p_Alpha_x_graph = gcurve(color=color.magenta, label="Alpha")

p_total_x_graph= gcurve(color=color.cyan, label="total")

# Other curves as needed

gdy = gdisplay(width=600, height=200, title='y Momentum vs Time', background=vector(0.8, 0.8, 0.8))

p_Au_y_graph = gcurve(color=color.yellow, label="Au")

p_Alpha_y_graph = gcurve(color=color.magenta, label="Alpha")

p_total_y_graph= gcurve(color=color.cyan, label="total")

g_Energy = gdisplay(width=600, height=200, title='Energy vs Time', background=vector(0.8, 0.8, 0.8))

p_Au_KE = gcurve(color=color.yellow, label="Au KE")

p_Alpha_KE = gcurve(color=color.magenta, label="Alpha KE")

p_PE = gcurve(color=color.black, label="Potential Energy")

p_Total_E= gcurve(color=color.cyan, label="Total Energy")

# Create a variable d to hold the minimum separation between the Au and Alpha

# Set it to the initial separation.

d = r0

# Set up the time counter and time step

t = 0

dt = 5e-25

# Continue until separation distance (r) is greater than the initial separation (r0),

# that is until the Alpha has finished interacting with the Au and is at least

# as far away as it was in the beginning.

while mag(r) <= r0:

rate(10000)

# Compute position of the alpha particle relative to the gold nucleus

r = Alpha.pos - Au.pos

# A - Compute the force on the alpha particle due to the gold nucleus

F =

# B - Update the momentum of the Alpha particle

Alpha.p =

# C - Update the position of the Alpha particle

Alpha.pos =

# D - Update the momentum of the Gold nucleus

Au.p=

# E - Update the position of the Gold nucleus

Au.pos =

t = t + dt

# Update the distance of closest approach

d = min(d,mag(r))

# Compute the kinetic and potential energies

K_Au = mag(Au.p)**2/(2*Au.m)

K_Alpha = mag(Alpha.p)**2/(2*Alpha.m)

U = oofpez*Au.q*Alpha.q/mag(r)

# Update the energy graph

p_Au_KE.plot(pos = (t, K_Au))

p_Alpha_KE.plot(pos = (t, K_Alpha))

p_PE.plot(pos = (t, U))

p_Total_E.plot(pos = (t, K_Au + K_Alpha + U))

# Update the x-component of momentum graph

p_Au_x_graph.plot(pos = (t, Au.p.x))

p_Alpha_x_graph.plot(pos = (t, Alpha.p.x))

p_total_x_graph.plot(pos = (t , Alpha.p.x + Au.p.x))

# Update the y-component of momentum graph

p_Au_y_graph.plot(pos = (t, Au.p.y))

p_Alpha_y_graph.plot(pos = (t, Alpha.p.y))

p_total_y_graph.plot(pos = (t, Alpha.p.y + Au.p.y))

# Use the dot product to find the scattering angle

D = acos(dot(Alpha.p , Alpha.p1) / (mag(Alpha.p) * mag(Alpha.p1))) # radians

# Print summary of results

print("Initial kinetic energy of alpha: ", K0/1.6e-13, " MeV")

print("Impact Parameter: ", b/1e-15, " fm")

print("Distance of closest approach: " , d/1e-15, " fm")

print("Scattering angle: ", degrees(D), " degrees")