%% PHYS 225 - Rutherford Lab (sample report)
% FirstName LastName
%%
disp(date)
%% Purpose
% The purpose of this lab is to create a simulation of Rutherford's
% experiment using Matlab; this will allow us to estimate the fraction of
% alpha particles that are scattered from a gold atom at steep angles.
%% Answers to handout questions
%%
% 1. Given Rutherford's findings about the structure of the atom, if you
% bombarded a gold atom with a spray of alpha particles, what fraction
% would you estimate to be scattered at steep angles? (We are looking for
% an estimate, not a complicated calculation, but you should think about
% the geometry of the situation. Note: look up appropriate quantities you
% need for your estimate). Briefly explain what quantities went into your
% estimate and how you chose them.
%%
% We can estimate the fraction of alpha particles that are scattered by
% finding the ratio between the area of the cross-section of the gold
% nucleus to the area of the entire atom's cross-section:
r = 7.3e-15; % nuclear radius in meters
R = 0.1441e-9; % atomic radius in meters
fs = pi*r^2/(pi*R^2); % fraction scattered
disp('Fraction scattered:')
disp(fs)
%%
% We find that
% the number of scattered alpha particles is $2.566x10^{-9}$ times the total
% number of particles incident on the gold atom.
%%
% 2. It can be shown that the projectiles follow hyperbolic paths and that
% the relationship between a projectile's impact parameter $b$ and its
% scattering angle $\theta$ is given by:
% $b = \frac{zZ}{2K} \frac{e^2}{4\pi \epsilon_0} cot(\theta /2)$
% where $z$ and $Z$ are the charge of the alpha
% particle and the gold nucleus, respectively, $e$ is the electron's charge,
% $\epsilon_0$ is the permittivity of vacuum, and $K$ is the kinetic energy of the alpha
% particle. Estimate $b$ for an alpha particle with kinetic energy 8.0 MeV
% that is scattered at an angle of 90°. Use this value of $b$ to refine your
% estimate for question 1.
%%
clear all
z = 2; % charge of alpha particle
Z = 79; % charge of gold nucleus
q = 1.602e-19; % electron charge in Coulombs
e0 = 8.85e-12; % permittivity of vacuum in F/m
K = 8e6*q; % kinetic energy in eV, transform to joules by multiplying by q
b = z*Z*q^2*cot(pi/4)/(2*K*4*pi*e0);
disp('Impact parameter:')
disp(b)
R = 0.1441e-9; % atomic radius in meters
disp('Fraction scattered:')
disp(b^2/R^2)
%%
% Using the formula, we find that $b = 1.423x10^{-14} m$ and the fraction of
% alpha particles that are scattered is $9.745x10^{-9}$.
%% Computation: Matlab code for simulating Rutherford's experiment
R = 14410; % atomic radius
b = 1.423; % impact parameter
n = 1e2; % number of random points (x10)
inside = 0; % counter for points inside b
outside = 0; % counter for points outside b
% Generate the random points
% Express the rectangular coordinates of x and y in terms of the
% point's radius and its angle relative to the positive x-axis. Take
% the square root of r1 so that the points will be evenly distributed
% throughout the circle.
for i=1:10
for i=1:n
r1 = rand;
r2 = rand;
x = R*sqrt(r1)*cos(2*pi*r2);
y = R*sqrt(r1)*sin(2*pi*r2);
if (x^2+y^2) < b^2 % if the point is within radius b add 1 to
inside = inside+1; % the number of points inside
else
outside = outside+1; % otherwise, add 1 to the number of points
end % outside
end
end % the second loop allows for more iterations
disp('Points inside:')
disp(inside)
disp('Points outside:')
disp(outside)
%% Results and analysis
% Code that generates table with data from several runs:
expected=9.745e-9; % from in-class calculations
nshot = [100; 1000; 5000; 10000; 62200]; % number of points 'shot'
pointsin = [0.5; 9.5; 45.3; 93; 597]; % points that fell inside b
fracin = pointsin./(nshot*1e6); % fraction of points that fell inside b
fracerror = 1-fracin/expected; % fractional error
data = horzcat(nshot, pointsin, fracin, fracerror); % all data collected for table
hFig1 = figure('Position',[0 0 600 250]);
hTable = uitable(hFig1, 'Data',data,... % make table with data
'ColumnName',{'n (in million)' 'Avg. points inside' 'Avg. fraction inside' 'Fractional error'}, ...
'Position',[10 10 420 150]);
%%
% Code for plotting average fraction inside vs. number of particles:
close(hFig1);
hFig2 = figure('Position',[1 1 500 300]);
plot(nshot,fracin,'.')
axis([0 70000 4e-9 10e-9])
hold on
plot(nshot,fracin)
title('Fig. 2 Average fraction inside vs. number of particles');
xlabel('Number of particles (x10^6)')
ylabel('Average fraction inside')
hold off
%%
% When n, the number of particles 'shot,' is small the results of the
% simulation tend to fluctuate because not enough data is being
% sampled. As n increases, however, the range of results across
% multiple trials becomes much smaller. Figure 2 shows that the average
% fraction scattered
% appears to be converging toward a fraction of approximately $9.0x10^{-9}$,
% which is close to the expected value of $9.752x10^{-9}$. Additionally,
% the cumulative total (out of 62.2 billion incident particles) yields
% a fraction of $9.598x10^{-9}$, which is also very close to the expected
% value.
%%
% Code for plotting average fractional error vs. number of particles:
close(hFig2);
hFig3 = figure('Position',[1 1 500 300]);
plot(nshot,fracerror,'.')
axis([0 70000 0 0.5])
hold on
nsqrt = [0:10:62200];
afe = 5./sqrt(nsqrt);
plot(nsqrt,afe)
title('Fig. 3 Average fractional error vs. number of particles');
xlabel('Number of particles (x10^6)')
ylabel('Average fractional error')
clear all
%%
% As displayed in Figure 3, the average fractional error decreases
% rapidly with increasing number of particles. Notably, this is very
% close to an inverse square
% root relationship. Therefore, I will postulate that
% AFE ~ $1/\sqrt(x)$ .
% According to this model, increasing the number of random points by a
% factor of $x$ causes the AFE to decrease by a factor of $1/\sqrt(x)$. With a
% significant amount of additional data, it would be possible to test
% and investigate this relationship.
%% Conclusion
% In this lab, I successfully used the Monte Carlo method in the Matlab
% computing environment to develop a simulation of the Rutherford
% experiment. When a sufficient number of points were sampled, the results
% of the simulation matched expectations very closely. Although a
% simulation such as this one cannot take the place of real data collection,
% it can help us validate the inferences that Rutherford made about the
% structure of the atom based on the data that he collected. Observing the
% strong correlation between the simulation's data and the expected results,
% I conclude that the simulation supports Rutherford's deductions.
% Additionally, I was able to form a reasonable postulate regarding the rate
% of convergence of the fractional error in the collected data. Further
% pursuit of this topic could shed light on the nature and limitations of
% the Monte-Carlo method.
%% Feedback
% As someone who had very little experience with Matlab coming
% into this course, I think this lab was a great introduction to some of
% Matlab's key features and its applications in modeling. This lab
% fulfills the third course objective because it focuses on modeling
% physics problems on the computer. In addition to the two class periods
% dedicated to this lab, I spent approximately eight hours collecting
% additional data, analyzing the data, and writing the lab report. This
% lab could be improved by adding a more detailed explanation of what is
% meant by the rate of convergence. From what I've heard, that portion of
% the lab was somewhat puzzling for myself and others.