Catching a Nuclear Smuggler

One of the terrorist threats our nation faces is an attack with a nuclear weapon smuggled onto US soil. In this laboratory we will use a Monte Carlo method to investigate some of the parameters that determine our ability to find such a device. The primary method is detecting the bomb's radioactivity;all nuclear-weapon materials like uranium and plutonium emit high-energy photons called gamma rays. The more gamma rays emitted,the easier it is to uncover the material. A terrorist can reduce this 'signature' by surrounding the radioactive material with lead or some other dense metal to block the passage of the gamma rays. This effect is called shielding. However,there is another shielding effect that we must consider if want want to detect the bomb's gamma-ray signature. Consider a sphere of 235-U (a uranium nucleus with a total of 235 neutrons and protons) metal. A gamma ray emitted by one of the uranium nuclei at the center of the sphere has to pass through a considerable amount of material before it gets out. The uranium of the sphere shields the nuclei in its interior. This 'self-attenuation' reduces the number of gamma rays that can leave the sphere and shrinks the radioactive signature. Below you will be asked to write a program to calculate the fraction of gamma rays that actually get out of a sphere of uranium. This is the self-attenuation of the nuclear material.

1.0 Some Necessary Commands

In this section several Mathematica commands that will be useful in your simulation are demonstrated.  You should execute the commands in the cells below to investigate the new functions.

1.1 Generating random numbers.

Random numbers can be generated with the following command.

xmin  =  -3;
xmax  =  +3;
Random[Real, {xmin, xmax}]

1.2 Loops.

There are two useful methods for doing iterative calculations in Mathematica. The first is the Do loop below. It works well when the iterator is simple and when you need to use nested loops.

\(\(\( (*\ 
      one - dimensional\ \(\(loop\)\(.\)\)\ *) \)\(
\)\(Do[
y\  = \ x + 3; 
Print["y=", 
        y, "."], 
{x, \(-2\), 4, 
        2.0}
]\)\)\)

\(\(\( (*\ 
      two - dimensional\ \(\(loop\)\(.\)\)\ *) \)\(
\)\(Do[
z\  = \ x*y; 
Print["x*y=", 
        z, " x=", x, " y=", 
        y, "."], 
{x, \(-2\), 4, 2.0}, {y, 0, 3, 
        1.5}]\)\)\)

The second loop method is the While function shown below. It is useful when the decision to continue the loop is more complex. It is a bit messier for nested loops.

\(\(\( (*\ 
      define\ the\ initial\ \(\(values\)\(.\)\)\ *) \)\(
\)\(\(x\  = \ \(-0.1\);\)

    \(test\  = \ Sin[x\^2];\)

    While[
      test\  < \ 0.8, 
x\  = \ 
        x + 0.2; 
test\  = \ 
        Sin[x\^2]; 
Print["x=", 
        x, "."]
]\)\)\)

1.3 An acceptance-rejection method in one dimension.

An acceptance-rejection method can be used to generate random numbers that have an exponential distribution. The While loop shown below can be used to test for acceptance or rejection of the point obtained with the Random function. The test in the While loop puts a limit on the number of iterations to prevent 'runaway' loops.

\(\(\( (*\ 
      Set\ things\ \(\(up\)\(.\)\)\ *) \)\(
\)\(\(argument\  = \ 0.75;\)

    \(maxloops\  = \ 5;\)

    \(ProbMax\  = \ 1.0;\)

    \(ProbMin = 0.0;\)

 (*\ 
      Apply\ acceptance - 
        rejection\ to\ select\ \(\(events\)\(.\)\)\ *) 

    \(Do[
TestDecay\  = \ 
          Random[Real, {ProbMin, ProbMax}]; 
\ \ \ \ \ probability\  = \ 
          Exp[\(-argument\)]; 
If[
          TestDecay\  < \ probability, 
          Print["Take this event."]\ ]; 
If[
          TestDecay\  > \ probability, 
          Print["Skip this one."]\ ], 
{nloops, 1, 
          maxloops}
];\)\)\)\)

1.4 Histogramming a distribution.

When we investigate statistical processes like nuclear decay we will often want to observe the frequency distribution of a series of Monte Carlo runs or simulations. The following code demonstrates how one can generate the frequency distribution in Mathematica.

\(\(\( (*\ 
      The\ next\ two\ commands\ can\ be\ executed\ only\ once
	\ \ or\ you\ will\ have\ to\ restart\ Mathematica . \ 
          They
	\ \ load\ new\ commands\ in\ from\ files\ on\ \(\(disk\)\(.\)\)\ *) \)\(
\)\(<< Graphics`Graphics`

    << Statistics`DescriptiveStatistics`\)\)\)

Next, enter some data. The command below shows how to enter an array of data in Mathematica.

data  = {0.0789609, -0.471797, -0.698464, -0.158695, 
          0.0848231, 2.39845, 1.36031, 1.07166, 1.06801, 
          0.954635, -1.47636, -1.03402, -0.172193, -0.0546502, -1.40035, 1.89033, -0.164824, 1.31812, 0.477657, 
          1.37246, 0.713923, 
          1.05531, -0.961621, -0.541414, -0.789451, 0.751251, 
          0.692603, 0.206209, 0.785431, 0.886077, 1.02105, -0.405319, 
          1.97813, 0.230765, -0.796516, 1.21065, -1.89966, 0.367685, 
          1.57743, -0.0878903, 0.240861, 0.737706, 
          0.0976328, -0.340478, 1.5002, 0.430084, 0.315298, 
          1.32552, -0.389734, -0.66366, 0.48413, 
          0.839891, -1.63951, -0.32257, 0.326307, 
          2.23439, -2.09491, -0.164018, 0.595124, -0.529316, 
          0.288771, 0.0797093, 1.0018, -0.795769, 0.966422, -1.45939, 
          0.95543, 0.682462, 0.122714, -0.388816, 1.59474, -1.39029, 
          1.24754, 
          1.38261, -0.233177, -1.18613, -0.307277, -0.229948, -1.34569, 0.0661231, 0.533871, 0.337921, -0.690251, -0.55572, 
          0.519533, 0.157858, -0.374646, -1.6623, 0.266872, 
          0.894534, -0.588928, -0.525838, -0.661636, 
          0.448432, -0.120686, -0.961073, -0.150906, 1.08303, 
          1.96417, 0.136119};

Make a histogram of the data array.

Histogram[data]

2.0 Selecting Points Randomly Inside a Sphere.

When an atomic nucleus inside our uranium sphere emits a gamma ray, it can come from any point inside the sphere. 

a. Write a short code using a Do loop to generate points randomly in a cube of side r centered at the origin and then select the ones that are within the volume of sphere of radius r centered at the origin .

b. What is the ratio of the volume of sphere of radius r to a cube of side r?

c. Calculate the number of points inside your sphere divided by the total number of 'throws'. Is the value what you expected?

3.0 Simulating Isotropic Emission.

When an atomic nucleus emits a gamma ray, it can come out in any direction. To simulate this process consider the following questions. 

a. What is the range covered by cosθ where θ is the angle between the direction of the gamma ray and z axis?

b. Generate a table of random values of cos θ. Make a histogram of your results. Does the histogram look right? Why or why not?

c. What is the range covered by the angle ϕ?

d. Generate a table of random values of ϕ. Make a histogram of your results. Does the histogram look right? Why or why not?

4.0 Simulating the Effect of Shielding.

From observation we know that the chance that a gamma ray will penetrate some material decreases exponentially with the thickness of the material.

I(r)  = (Subscript[I, 0])(e^(-μr))

where r is the thickness of materials, I is the intensity, and μ is called the range parameter and is the product of the stopping power and the density of the material.

a. For μ=0.88 (cm^(-1))and r=1.4 cm, write a code using the acceptance-rejection method that will select gamma rays that penetrate shielding of thickness r. Test it out for 100 events.

b. Modify your code in Part a so you can calculate the number of total gamma rays sent through the shielding and the number of gamma rays that penetrated the shielding. Use these results to calculate the fraction of gamma rays that made it through.

c. Calculate the fraction of gamma rays that penetrate the shielding for several values of r and plot these results as a function of r. Does your plot agree with I(r)?

5.0 Calculating Self-Attenuation.

The attenuation of a beam of gamma rays in some material is an exponentially decreasing function described by

I(r) = (Subscript[I, 0]) (e^(-μr))

where I(r) is the intensity of gamma rays still left in the beam, μ is called the range parameter, r is how far the gamma rays have penetrated the material, and (Subscript[I, 0])is the intensity of gamma rays in the beam at r=0.

We will now develop the algorithm for  determining the self-attenuation of a sphere of nuclear material. Use the following parameters to define the problem.

nthrows =20		- number of Monte Carlo rolls of the dice at a single point in the sphere.
mU = 25000.0		- mass of 235-U pit in grams.
rhoU = 19.05		- density of 235-U grams per cubic centimeter.
range = 0.876		- range parameter in U in (cm^(-1)).

The questions below are intended to guide you in developing the algorithm. It is important to understand the solution to the problem before you start writing code.

a. What is the radius of the uranium sphere?

b. Consider a gamma ray that is emitted isotropically from some random point in the uranium sphere. Make a sketch showing the uranium sphere, the emission point of the gamma ray, and the trajectory of the gamma ray. Assume the gamma ray leaves the sphere.

c. Now consider an algorithm where you step along the trajectory of the gamma ray. At each step, how would you calculate the distance from the origin to the current position of the gamma ray at that particular step along the trajectory?

d. What is the condition for deciding if the gamma ray has reached the surface of the sphere?

e. Once the gamma ray has reached the surface, what is the expression for the probability the gamma ray has not been absorbed by the uranium sphere?

f. Using the ideas discussed above, carefully describe in words and symbols the algorithm for emitting a gamma ray from a random position inside the uranium sphere, tracking it's passage to the surface of the sphere, and deciding if that gamma ray leaves the sphere or is absorbed.

g. Implement the algorithm you describe in Part f and calculate the self-attenuation of the uranium sphere.