# Homework Help: A Monte Carlo script to simulate the activation and decay of cobalt

1. Aug 15, 2017

### cnet

1. The problem statement, all variables and given/known data

A sample of cobalt was previously measured in a neutron scattering experiment. As a result, the sample has become activated. Write a Monte Carlo script to simulate the activation and decay of a Co sample during and after a neutron scattering experiment. It is recommended to do this in Matlab.

The Matlab script should calculate and graph the number N1(t), N2(t) and N3(t) as a function of time t, where N1(t) is the number of atoms of 60Co, N2(t) is the number of atoms of 61Co, and N3(t) is the number of atoms of 59Co respectively.

2. Relevant equations

N(t) = N(0)e-λt
3. The attempt at a solution

To start with, I used the given information neutron ﬂux is 1×1016 n/s/cm2. The neutron absorption cross section is 98.223×10−24 cm2. The half life of 60Co is 5.27 y. The half life of 61Co is 1.5 h. Gamma decay is the only decay channel available for activated Co and the sample area is exactly 1 cm2.

However, the simulation is not just about calculation and is supposed to be a process the condition of which is subject to all variables affecting the number of each isotope. Does anyone know how to construct this simulation for the activation and decay of activated Co?

Regards,

2. Aug 15, 2017

### FactChecker

The program should step through time modeling the decay of each type at each time step. Because it is a Monte Carlo simulation, you should have a random decay of each type with the correct mean and variance for one time step.

3. Aug 16, 2017

### haruspex

It might be more efficient to generate times to the next decay, rather than step through equal time intervals.
If the RNG produces a number uniformly in [0,1), subtract from 1 to get (0,1] then apply ln() suitably in conjunction with the rate parameter.

4. Aug 17, 2017

### FactChecker

An interesting alternative. I hadn't thought of an event driven simulation. That presents some complications of scheduling events and figuring out what happens when an event adds a particle of another type which had not been included in that type's next decay event scheduling.

5. Aug 17, 2017

### haruspex

The way I was thinking of it, if there are three different types of event then you can generate a random time to the next event of any type by summing the rates. You can then decide which type it was by a second random number, comparing it to a set of thresholds spaced according to the ratios of the rates. Would that induce a bias?

Last edited: Aug 17, 2017
6. Aug 17, 2017

### FactChecker

I think you are right. I stand corrected. They are all assumed to be Poisson processes, so at each event time you can figure the next event time for all types without regard for the past. So there is no need to keep a schedule of upcoming decays events -- you only need to know the time of the next event.

Except that I don't think you can sum the rates of Poissons with different rates. You would have to figure the random time of the next decay for each type separately and determine the earliest one.

7. Aug 23, 2017

### cnet

This is difficult.

The system has 3 isotopes 60Co, 61Co and 59Co. The former two are radioactive, but 59Co is stable.

To begin with, I take that there were 106 Co atoms in the sample. So:

N(0) = 1e6

My questions are:
(1) Do I have to use the relative isotope abundances for natural occurring Co to get the initial numbers of N1(0), N2(0) and N3(0)?

(2) When the neutron beam was on, N1(t) and N2(t) for 60Co & 61Co increase as some of the atoms become activated, but at the same time they decay to stable N3(t) for 59Co via gamma decay. How to take this into account in the simulation script?

(3) Can I assume, after the neutron beam was turned off, N1(t) and N2(t) for 60Co & 61Co would just decrease via gamma decay, and N3(t) for 59Co increases?

I have a script for a single decay process, but in this case there are three isotopes and their numbers are inter-related I think.

8. Aug 23, 2017

### FactChecker

That seems like a good starting state unless you are given some other information.
If some types decay into other types, you should keep track of that. Then the calculations for the next time step will be based on the new amounts.
If that is what happens, then that is what you should model in the simulation. At a time step when the neutron beam is on you would have one set of rates and when it is off you would have a different set of rates.
Yes, they are inter-related. I don't know what simulation language you are using, but there are some that specifically model continuous process like this. Others would allow you to step through in fixed time steps. Others can adjust time steps to be shorter when things are changing rapidly and longer when not much is happening. From the sound of it, this seems like you could just use a general purpose language like C, Fortran, etc. to step through time and simulate the process.

PS. Clearly there are too many atoms to use an event-driven simulation as suggested in post #3, so a time step approach is best.

Last edited: Aug 23, 2017
9. Aug 24, 2017

### cnet

I have been told that all atoms initially in the Co sample are stable 59Co. I have written a Matlab script for the neutron activation of 59Co to 60Co/61Co using the time step approach. (attached: https://files.fm/u/kutcymnn)

But I don't know in each time step how many 59Co atoms decayed to 60Co and how many to 61Co. How can I inter-relate these numbers?

10. Aug 24, 2017

### FactChecker

Oh! I didn't realize that you didn't have that information. So this is more a nuclear physics question than a simulation question. I can not help and will withdraw.

11. Aug 24, 2017

### cnet

No, I just realised it has to be done in two steps. 59Co-60Co and then 60Co-61Co, for the neutron activation time. I have the former in the script. I don't think it's allowed for 59Co to go directly to 61Co. But adding the 2nd step into the script is a problem for me as I am not a professional programmer.

12. Aug 24, 2017

### FactChecker

If there is a sequence of decay steps, you can take care of that in the sequence of time steps. At each time step, only worry about the direct conversions of one type to another. On every step, some percentage of one type converts directly to other types. On the same step, some percentage of those other types convert directly to still other types, including the atoms it gained the time step before. Make the time steps short enough that you don't lose too much accuracy in your percentages due to breaking the decay steps into separate time steps.