# Physics 3730/6720 Assignment 6

Reading and references:

For each of the following exercises create the specified files with your
answer(s). Submit your homework using the course **submit**
utility as usual.

#### Exercise 1. Monte Carlo simulations of radioactive decay

For this exercise you will write a code
called **decay.py** that simulates
radioactive decay of a bunch of atoms. You will do this using random
numbers, so this "experiment" can be called a Monte Carlo
simulation. The inputs to the code are (in this order, one per line):
- Number of atoms
**N**.
- Lifetime
**tau**.
- Simulation time interval
**dt**.
- Duration of the experiment
**T**.
- Random number seed
**seed**.

The code should start by seeding the random number generator. Then,
during a time interval **dt**, it should visit every atom once.
The probability that the atom decays during that time interval
is **dt/tau** (as long as this ratio is much less than one). So
for each atom, it should get a random number **r** on the
interval **[0,1)**. If that number is less than **dt/tau**,
then the atom decays. Otherwise, it does not. Do this for all of the
atoms during that time interval. As you go along, count the number of
decays. Reduce the number of atoms by that number. Then go on to the
next time interval and repeat the whole process.
Your code should print two numbers per line: the time and the number
of atoms that decayed in the time interval at that time. To be
precise, print the time at the end of the interval.

#### Exercise 2. Exponential decay

Run your program with the input
file **~p6720/examples/radioactive/indecay**. Plot the output
number of decays vs time (use points), and then, in the same plot, plot
the theoretical prediction for the number of decays (use lines),

N*dt/tau exp(-t/tau) .

Convert the plot to a pdf file
called **decay.pdf** for submitting. You
may use gnuplot or (if you are brave) pyplot for making the plot from
your calculation. If you use gnuplot, you can combine two plots by
separating the plot specifications with a comma. That is, if **plot
A** and **plot B** make two separate plots, then
**plot A, B** combines them into one plot. To plot the function
above, write it as a function of **x**, not **t**. If you need
to convert from Postscript to pdf, use the Unix
utility **epstopdf**. If you use pyplot, save the figure to a pdf
file using **plt.savefig("figure.pdf")**.

#### Exercise 3. Poisson distribution

Run your program with the input
file **~p6720/examples/radioactive/inpoisson**. Check that the
number of atoms has not changed very much during this experiment.
Create a histogram that shows the number of time intervals in your
experiment that had 0 decays, the number that had 1 decay, etc. You
may use the course **hist** utility or try your hand
with **matplotlib.pyplot.hist**. You should arrange so that the
histogram bins are centered on the integers.

The theoretical expectation is that the histogram is proportional to a
Poisson distribution, based on the expected average number of decays
in the interval **dt**, namely

d = N*dt/tau

The theoretical expression is
N(k) = M d^{k} e^{-d}/k!

where **M = T/dt** is the number of time intervals in your
experiment and **k!** is the factorial **k(k-1)(k-2)...1**. So
on the same plot as your histogram, plot this prediction (as points).
The easiest way to do this is to write a small Python code that writes
out **k** and **N(k)** for integer **k** from 0 to 20 or so.
Note that **math.factorial(k)** evaluates the factorial. Put the
output numbers in a file, so you can use it as input
to **gnuplot**.
Convert the plot to a pdf file
called **poisson.pdf** for submitting.

*This page is maintained by:
Carleton DeTar detar@physics.utah.edu
*** Last modified 19 October 2017 **