# Physics 3730/6720 Assignment 6

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.
• 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. (Do this with a for loop with one iteration per surviving atom.) 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 dk 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.