Rutherford scattering

Rutherfords spredningsforsøg er et eksempel på en elastisk kollision, da atomkernerne ikke rammer hinanden, men kun er udsat for hinandens elektriske ladning.
Kraften der virker mellem kernerne er beskrevet ved Coulumbs lov: \begin{equation} \vec{F} = \frac{1}{4 \pi \epsilon_o} \frac{q_1 q_2}{r^2} \hat{r} \end{equation}


% Rutherford scattering experiment
% (c) Morten Dam Joergensen, 2011
%
% Background: [F&Y] p. 1294, [Rel] Chap. 6
% 
% Goals:
%  - Recreate the famous Rutherford scattering experiment, where alpha
%  particles are scattered on gold nuclei.
%  The scattering angles observed showed Rutherford that the mass of the
%  atom was concentrated at a small point with a large charge density.
% 
%  - Use and improve the simulation to determine:
% 
%  --- Use a Histogram for the following questions,
%     - The initial momentum of the alpha particle
%     - The the final momentum of the alpha particle a good distance away from
%     the gold nucleus.
%     - The final momentum of the gold nucleus long after the interaction
%     - What is the final kinetic energy of the alpha particle
%     - What is the final kinetic energy of the gold nucleus
%     - Histogram the closest distance between the alpha particle and the
%     gold nucleus (hint, use the max(val1, val2) function!
%     - For the radii given below, do the particles touch?
%     - Estimate a better delta-t value based on the initial velocity of
%     the alpha particle (remember that delta-t should be much smaller)

% --- Changes to the plotting:
%     - Create a plot that shows the trajectories of the alpha particle at
%     the following angles: angles = -10:2:10; (like fig 39.13 in [F&Y])
%     - Try changing the following parameters:
%       - Target material (Try Al, Cu, Pb, U)
%       - Kinetic energy of the alpha particle (1, 5, 10, 20, ... MeV keep it non relativistic)
%       - Try specific imparact parameters or use the angle variables, to
%       see what happens at specific values.
%       - Try playing with dt, what values makes sense.


clear;clc; figure(1); clf; figure(2); clf;

eV = 1.602e-19; % (J) Electronvolt in joules
MeV = 1e6 .* eV;  % (J) MeV in joules
c = 2.99792e8; % (m/s) Speed of light 

konst = 9e9; % (N*m^2/C^2) Coulomb constant: 1/(4pi*eps0) 
e = 1.602176565e-19; % (C) Elementary charge in coulombs

% Coulomb's law (electrostatic force)
f_elec = @(p1, p2) (konst .* (p1.Z*e .* p2.Z*e) ./ (norm(p1.pos-p2.pos).^2)) .* (p1.pos-p2.pos)/norm(p1.pos-p2.pos);

% Relativistic helper functions
beta = @(v) v/c; % Velocity in terms of light speed
gamma = @(v) 1/sqrt(1 - beta(v)^2); % Lorentz' gamma factor in terms of velocity
gamma_KE = @(KE, mass) (KE + (mass * c^2)) / (mass * c^2);  % Lorentz' gamma factor:    gamma=(E0 + KE)/E0 = (mc^2 + KE)/mc^2
beta_from_gamma = @(g) sqrt(1 - 1/g^2); % Beta from gamma factor:  beta = sqrt( 1 - 1/gamma^2 )
speed_from_KE = @(KE, mass) c * beta_from_gamma(gamma_KE(KE*MeV, mass)) * gamma_KE(KE*MeV, mass);

% Velocities in terms of kinetic energy, rest mass and angle
velocity_from_KE_rel = @(KE, mass, angle) [speed_from_KE(KE, mass) * cosd(angle), speed_from_KE(KE, mass) * sind(angle)];

% non-relativistic velocity (quasi-relativistic perhaps...  try comparing with the relativistic version for various KE)
velocity_from_KE = @(KE, mass, angle) [sqrt((2*KE * MeV)/mass) * cosd(angle), sqrt((2*KE * MeV)/mass) * sind(angle)];

N_tries = 200;
% final_scattering_angles = zeros(N_tries);
for tries = 1:N_tries
    
     % Glitter and gold (197Au)
    Au = struct();
    Au.N = 118; % Neutron count
    Au.Z = 79; % Proton count
    Au.q = Au.Z*e; % (C) Charge
    Au.pos = [5e-13, 0]; % Initial position
    Au.radius = 8e-15; % (m) Radius
    Au.p = [0, 0]; % (m*v) Initial momentum
    Au.mass = 3.3e-25; % (kg) rest mass (Rutherford estimate)

    % Alpha particle (Helium nucleus, 2H^2+)
    Alpha = struct();
    Alpha.N = 2; % Neutron count
    Alpha.Z = 2; % Proton count
    Alpha.q = Alpha.Z*e; % (C) Charge
    
    Alpha.pos = [0, (rand-0.5)*1e-13]; % Inital position (socastic)
    
    Alpha.radius = 2e-15; % (m) Radius
    Alpha.mass = 6.64e-27; % (kg) rest mass
    Alpha.angle = 0; % (degrees) angle from horisontal plane
    Alpha.KE = 10; % Kinetic energy in MeV
    Alpha.vel = velocity_from_KE(Alpha.KE , Alpha.mass, Alpha.angle); % (m/s)
    Alpha.vel_rel = velocity_from_KE_rel(Alpha.KE, Alpha.mass, Alpha.angle);% (m/s)
    Alpha.p = Alpha.vel .* Alpha.mass;

    dt = 1e-22; % (Very) small timestep!

    stop_dist = 0.9e-12; % When to stop simulation, is the the best way to do it?

    while norm(Au.pos-Alpha.pos) <= stop_dist

        % Columb force
        f_c = f_elec(Alpha, Au);

        % Update momentum
        Alpha.p = Alpha.p + dt.* f_c;
        Au.p = Au.p - dt.* f_c;

        % Update position of alpha particle
        Alpha.pos = Alpha.pos + dt.* (Alpha.p ./ Alpha.mass);

        % Update position of Gold nucleus
        Au.pos = Au.pos + dt.* (Au.p./Au.mass);

        % Draw
        figure(1); clf
        axis equal
        axis([-1.1*stop_dist, 1.1*stop_dist, -1.1*stop_dist, 1.1*stop_dist])

        hold on
        plot([-1.1*stop_dist, 1.1*stop_dist],[ 0, 0],'--k')
        plot([5e-13, 5e-13],[ -1.1*stop_dist, 1.1*stop_dist],'--k')
        drawNucleus(Alpha.Z, Alpha.N, Alpha.pos,'b');
        drawNucleus(Au.Z, Au.N, Au.pos,'k');
        text(-stop_dist, -stop_dist,['Scattering angle=',num2str(rad2deg( atan2(Alpha.p(2), Alpha.p(1)) ))])
  
        
%       pause(0.001)
    end
    % See http://en.wikipedia.org/wiki/Atan2
    final_scattering_angles(tries) = abs(rad2deg(atan2(Alpha.p(2), Alpha.p(1)))); %(atand(Alpha.p(2)/Alpha.p(1)));
    figure(2)
    hist(final_scattering_angles,20)
end
Download: rutherfordscattering.m
function g = drawNucleus(Z,N, pos,drawstring)
   r_nuclei = 1.3e-15; %m Radius of a single nuclei
   scale_factor = 1;
   
   vol_nuclei = 4/3 * pi * r_nuclei.^3;
   vol_total = (Z+N) * vol_nuclei;
   
   packing_fraction = 0.74; % http://en.wikipedia.org/wiki/Close-packing_of_spheres
   
   vol_packed = ((1-packing_fraction) * vol_total) + vol_total;
   
   r_total = (3*vol_packed / 4*pi).^(1/3);
   
   
    [x,y,z] = cylinder(r_total * scale_factor,200);
    g = plot(pos(1) + x(1,:), pos(2) + y(1,:), drawstring);
end
Download function drawnucleus.m