A Monte Carlo script to simulate the activation and decay of cobalt

[cnet]

Homework Statement

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 ⁶⁰Co, N2(t) is the number of atoms of ⁶¹Co, and N3(t) is the number of atoms of ⁵⁹Co respectively.

Homework Equations

N(t) = N(0)e^-λt

The Attempt at a Solution

To start with, I used the given information neutron flux is 1×10¹⁶ n/s/cm². The neutron absorption cross section is 98.223×10⁻²⁴ cm². The half life of ⁶⁰Co is 5.27 y. The half life of ⁶¹Co is 1.5 h. Gamma decay is the only decay channel available for activated Co and the sample area is exactly 1 cm².

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,

 

[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.

 

[haruspex]

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.

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.

 

[FactChecker]

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.

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.

 

[haruspex]

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.

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?

 

[FactChecker]

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?

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.

 

[cnet]

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.

This is difficult.

The system has 3 isotopes ⁶⁰Co, ⁶¹Co and ⁵⁹Co. The former two are radioactive, but ⁵⁹Co is stable.

To begin with, I take that there were 10⁶ 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 ⁶⁰Co & ⁶¹Co increase as some of the atoms become activated, but at the same time they decay to stable N3(t) for ⁵⁹Co 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 ⁶⁰Co & ⁶¹Co would just decrease via gamma decay, and N3(t) for ⁵⁹Co 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.

 

[FactChecker]

This is difficult.

The system has 3 isotopes ⁶⁰Co, ⁶¹Co and ⁵⁹Co. The former two are radioactive, but ⁵⁹Co is stable.

To begin with, I take that there were 10⁶ 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)?

That seems like a good starting state unless you are given some other information.

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

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.

(3) Can I assume, after the neutron beam was turned off, N1(t) and N2(t) for ⁶⁰Co & ⁶¹Co would just decrease via gamma decay, and N3(t) for ⁵⁹Co increases?

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.

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.

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.

 

[cnet]

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.

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

But I don't know in each time step how many ⁵⁹Co atoms decayed to ⁶⁰Co and how many to ⁶¹Co. How can I inter-relate these numbers?

 

[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.

 

[cnet]

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.

No, I just realized 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.

 

[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.