In this article I will briefly explain radioactive decay and the concept of half-lives, and then go on to write Python code to illustrate the underlying mathematics.
Radioactive Decay and Half-Lives
All atoms of a particular chemical element have a fixed number of protons in their nuclei, and this is what defines an atom as being of a particular element. However, nuclei can also contain varying numbers of neutrons, and each combination of proton and neutron numbers is called a nuclide.
Some of these nuclides are unstable and therefore radioactive, and are known as radionuclides. There are 339 naturally occurring nuclides on earth (many more have been created in particle accelerators or nuclear reactors) of which 252 are believed to be stable as no decay has ever been observed. The other 87 are unstable of which 34 are primordial radionuclides, the rest being created by the decay of primordial radionuclides as part of a so-called decay chain.
The decay of an individual atom of a radioactive nuclide is random and unpredictable. However, the time for which the probability of an atom decaying is 0.5 can be established by measuring a sufficiently large sample of the nuclide; this time is known as the nuclide's half-life.
Given a reasonable sample size, the nuclide will decay by approximately 50% during its half-life. The remaining amount of the nuclide will decay by approx. 50% during the next half-life, leaving 25% of the original amount and so on.
This graph shows the amount remaining across 5 half-lives of a radionuclide. The actual units being measured are not shown as they could be the amount of radiation, or the mass or number of atoms of the nuclide remaining.
That was a very brief introduction to the concept of the half-life as applied to radioactive decay. The main purpose of this article is to examine the mathematical aspects and to demonstrate them with some Python code.
The Project
This project consists of the following files:
- radioactivenuclides.py
- halflifeutilities.py
- halflife.py
The files can be downloaded as a zip from the Downloads page, or you can clone/download the Github repository if you prefer.
Source Code Links
The radioactivenuclides.py file contains a dictionary of the 34 primordial radionuclides, including their half-lives.
halflifeutilities.py includes functions for the following tasks:
- Print the nuclides in radioactivenuclides.py
- Calculate the amounts of a radionuclide left across a number of half-lives
- Print the above list
- Calculate the amount of a radionuclide left after a specified time
- Calculate the amount of time taken for a radionuclide to decay to a certain amount
halflife.py simply calls the functions in the previous file.
radioactivenuclides.py
Firstly let's take a look at the data in radioactivenuclides.py.
radioactivenuclides.py (truncated - showing 3 of 34)
radioactive_nuclides = { "tellurium128": { "element": "tellurium", "elementsymbol": "Te", "atomicmass": 128, "significand": 2.2, "exponent": 24 }, "xenon124": { "element": "xenon", "elementsymbol": "Xe", "atomicmass": 124, "significand": 1.8, "exponent": 22 }, "krypton78": { "element": "krypton", "elementsymbol": "Kr", "atomicmass": 78, "significand": 9.2, "exponent": 21 },
This file is simply data coded as dictionaries inside a dictionary. As I mentioned it contains the 34 primordial radionuclides but you can add to them if you wish.
The keys to the outer dictionary are the names of the elements concatenated with the nuclide's atomic mass, ie. the total number of protons and neutrons in the nucleus.
The inner dictionaries contain the following data:
- element - the element's name
- elementsymbol - the 1 or 2 letter element symbol
- atomicmass - the the number of protons and neutrons in the nucleus
- significand and exponent - the components of the scientific form of the nuclide's half-life in years, for example tellurium 128's half-life is 2.2 x 1024
halflifeutilities.py
This is the halflifeutilities.py file.
halflifeutilities.py
import math import radioactivenuclides as ran def _sci_not_to_num(significand, exponent): return significand * 10 ** exponent def printnuclides(): """ Print details of the the contents of radioactivenuclides as a table """ print("-" * 71) print("| Key ", end="") print("| Nuclide symbol ", end="") print("| Nuclide name ", end="") print("| Half-life (years) |") print("-" * 71) for key, value in ran.radioactive_nuclides.items(): print(f"| {key.ljust(14, ' ')}|", end="") print(f" {(str(value['atomicmass']) + value['elementsymbol']).ljust(15, ' ')}|", end="") print(f" {(value['element'] + ' ' + str(value['atomicmass'])).ljust(15, ' ')}|", end="") print(f" {(str(value['significand']) + 'undefined10^' +str(value['exponent']) ).ljust(18, ' ')}|") print("-" * 71) def decay_list(nuclidekey, halflives, startingvalue = 1048576): """ For the specified nuclide creates a list of amounts remaining after each half-life has elapsed """ dl = [] nuclide = ran.radioactive_nuclides[nuclidekey] halflife = _sci_not_to_num(nuclide['significand'], nuclide['exponent']) for i in range(0, halflives + 1): remaining = 0.5**i dl.append({ "years_elapsed": halflife * i, "remaining_decimal": remaining, "remaining_amount": startingvalue * remaining }) return dl def decay_table(nuclidekey, halflives, startingvalue = 1048576): """ For specified nuclide prints table of half lives and years, and amounts remaining as decimals and amounts. """ dl = decay_list(nuclidekey, halflives, startingvalue) tablewidth = 75 nuclide = ran.radioactive_nuclides[nuclidekey] halflife = _sci_not_to_num(nuclide['significand'], nuclide['exponent']) print(f" {nuclide['element']} {nuclide['atomicmass']}\n") print(f" Half-Life {nuclide['significand']}undefined10^{nuclide['exponent']} years\n") print("-" * tablewidth) print("| Elapsed | Remaining |") print("-" * tablewidth) print("| Half-Lives ", end="") print("| Years ", end="") print("| Decimal | Amount |") print("-" * tablewidth) for index, item in enumerate(dl): print(f"| {index:<11d}|", end="") print(f" {item['years_elapsed']:<17e}|", end="") print(f" {item['remaining_decimal']:<19.12f} |", end="") print(f" {item['remaining_amount']:<17.8g} |") print("-" * tablewidth) def at_time(nuclidekey, years, startingvalue = 1048576): """ Calculates amount of nuclide remaining after specified number of years. Result is a dictionary with keys remaining_decimal and remaining_amount """ nuclide = ran.radioactive_nuclides[nuclidekey] halflife = _sci_not_to_num(nuclide['significand'], nuclide['exponent']) exponent = -(years / halflife) remaining_decimal = 2**exponent remaining_amount = startingvalue * remaining_decimal return { "remaining_decimal": remaining_decimal, "remaining_amount": remaining_amount } def time_to(nuclidekey, remaining_decimal): """ Calculates the amount of time for a nuclide to decay to the specified proportion """ nuclide = ran.radioactive_nuclides[nuclidekey] halflife = _sci_not_to_num(nuclide['significand'], nuclide['exponent']) factor = -(math.log2(remaining_decimal)) return halflife * factor
_sci_not_to_num
As I just mentioned the half-lives of nuclides in radioactivenuclides.py are stored as significands and exponents of the scientific notation. This function calculates the actual number from these values for use in calculations.
printnuclides
This function simply prints the contents of radioactivenuclides.py in a table format. The "Nuclide symbol" consists of the atomicmass and elementsymbol values concatenated. "Nuclide name" is the concatenation of element and 128 with a space between them.
The half-lives are shown in the format 2.2×10^24. Python uses the format 2.2e+24 which I dislike so I stored half lives as significands and exponents to make output simpler. I could have done things the other way round, string the actual value and converting to significands and exponents just for output. However, this can generate slightly inaccurate results with decimals such as .000000000001 or .999999999999.
decay_list
The decay_list function takes a nuclide key (as used in the radioactive_nuclides dictionary), a number of half lives and a starting value. It then calculates the fraction remaining after each half-life, adding a dictionary of the number of half-lives elapsed, and the fraction and actual amount left to a list. The startingvalue and remaining_amount can represent radiation, mass or number of atoms.
The fraction remaining is calculated by raising 0.5 to the power of the number of half-lives elapsed, which is represented by the loop index. The following table shows how this works.
| Half-Lives | 0.5 ^ half-lives |
|---|---|
| 0 | 1 |
| 1 | 0.5 |
| 2 | 0.25 |
| 3 | 0.125 |
| 4 | 0.0625 |
| 5 | 0.03125 |
| 6 | 0.015625 |
| 7 | 0.0078125 |
| 8 | 0.00390625 |
decay_table
Here we print the data generated by decay_list. After printing the name of the nuclide and its half-life we print some headings and then the data itself within a loop.
at_time
This function takes a nuclide, number of years and a starting value and calculates the remainder (both fraction and amount) after the specified number of years. This table displays an example of the calculation for the arguments shown.
| nuclidekey | tellurium128 |
|---|---|
| years | 1,100,000,000,000,000,000,000,000 |
| startingvalue | 1,048,576 |
| halflife | 2,200,000,000,000,000,000,000,000 |
| exponent = -(years / halflife) | 0.5 |
| remaining_decimal = 2^exponent | 0.7071067811865476 |
| remaining_mount | 741455.2001894653 |
The number of years is 0.5 of the half-life, and if you look at the graph above you can see the curve is at about 0.7 at this point on the x-axis.
time_to
The time_to function takes a nuclide key and a decimal fraction, and calculates how long it takes the radionuclide to decay to that fraction. This is an example of arguments and calculated values.
| nuclidekey | tellurium128 |
|---|---|
| remaining_decimal | 0.5 |
| halflife | 2,200,000,000,000,000,000,000,000 |
| factor = -(math.log2(remaining_decimal)) | 1 |
| halflife * factor | 2,200,000,000,000,000,000,000,000 |
I chose 0.5 for the value of remaining_decimal to give an obvious result which is the half-life, and as you can see we have reached the correct result.
halflife.py
The halflife.py code is short and simple, and just runs the functions from halflifeutilities.py.
halflife.py
import halflifeutilities as hlu import radioactivenuclides as ran def main(): print("-------------------------------------") print("| codedrome.com |") print("| Half Life of Radioactive Nuclides |") print("-------------------------------------\n") hlu.printnuclides() # hlu.decay_table("uranium235", 8, 1048576) # nuclide = ran.radioactive_nuclides["tellurium128"] # halflife = nuclide['significand'] * 10 ** nuclide['exponent'] # at = hlu.at_time("tellurium128", halflife * 0.5, 1048576) # print(at) # tt = hlu.time_to("tellurium128", 0.5) # print(f"{tt} Years") # print(f"{tt:,.0f} Years") # print(f"{tt / 1000000:,.0f} Million Years") # print(f"{tt / 1000000000:,.0f} Billion Years") # print(f"{tt / 1000000000000:,.0f} Trillion Years") if __name__ == "__main__": main()
The first two sections of code simply call functions, the second, decay_table, with your choice of arguments.
For at_time I have first picked up the half life of tellurium 128 and then multiplied it by 0.5 for the years argument, before printing the result.
After calling time_to I have printed the result in no less than 5 formats. The first just prints the result in whatever way Python sees fit, whereas the rest print it in ever-larger multiples of years.
Now we can run the code like this:
Running the Code
python3 halflife.py
This is the output of the printnuclides function. We can see here the scientific notation format I have used.
Program Output: printnuclides
-----------------------------------------------------------------------
| Key | Nuclide symbol | Nuclide name | Half-life (years) |
-----------------------------------------------------------------------
| tellurium128 | 128Te | tellurium 128 | 2.2×10^24 |
| xenon124 | 124Xe | xenon 124 | 1.8×10^22 |
| krypton78 | 78Kr | krypton 78 | 9.2×10^21 |
| xenon136 | 136Xe | xenon 136 | 2.38×10^21 |
| germanium76 | 76Ge | germanium 76 | 1.8×10^21 |
| barium130 | 130Ba | barium 130 | 1.2×10^21 |
| selenium82 | 82Se | selenium 82 | 1.1×10^20 |
| cadmium116 | 116Cd | cadmium 116 | 3.1×10^19 |
| calcium48 | 48Ca | calcium 48 | 2.3×10^19 |
| bismuth209 | 209Bi | bismuth 209 | 2.01×10^19 |
| zirconium96 | 96Zr | zirconium 96 | 2×10^19 |
| tellurium130 | 130Te | tellurium 130 | 8.8×10^18 |
| neodymium150 | 150Nd | neodymium 150 | 7.9×10^18 |
| molybdenum100 | 100Mo | molybdenum 100 | 7.8×10^18 |
| europium151 | 151Eu | europium 151 | 5×10^18 |
| tungsten180 | 180W | tungsten 180 | 1.8×10^18 |
| vanadium50 | 50V | vanadium 50 | 1.4×10^17 |
| cadmium113 | 113Cd | cadmium 113 | 7.7×10^15 |
| samarium148 | 148Sm | samarium 148 | 7×10^15 |
| neodymium144 | 144Nd | neodymium 144 | 2.29×10^15 |
| osmium186 | 186Os | osmium 186 | 2×10^15 |
| hafnium174 | 174Hf | hafnium 174 | 2×10^15 |
| indium115 | 115In | indium 115 | 4.4×10^14 |
| gadolinium152 | 152Gd | gadolinium 152 | 1.1×10^14 |
| platinum190 | 190Pt | platinum 190 | 6.5×10^11 |
| samarium147 | 147Sm | samarium 147 | 1.06×10^11 |
| lanthanum138 | 138La | lanthanum 138 | 1.02×10^11 |
| rubidium87 | 87Rb | rubidium 87 | 4.97×10^10 |
| rhenium187 | 187Re | rhenium 187 | 4.12×10^10 |
| lutetium176 | 176Lu | lutetium 176 | 3.764×10^10 |
| thorium232 | 232Th | thorium 232 | 1.406×10^10 |
| uranium238 | 238U | uranium 238 | 4.471×10^9 |
| potassium40 | 40K | potassium 40 | 1.25×10^9 |
| uranium235 | 235U | uranium 235 | 7.04×10^8 |
-----------------------------------------------------------------------
If you comment out printnuclides in main and uncomment decay_table this is the output.
Program Output: decay_table
uranium 235
Half-Life 7.04×10^8 years
---------------------------------------------------------------------------
| Elapsed | Remaining |
---------------------------------------------------------------------------
| Half-Lives | Years | Decimal | Amount |
---------------------------------------------------------------------------
| 0 | 0.000000e+00 | 1.000000000000 | 1048576 |
| 1 | 7.040000e+08 | 0.500000000000 | 524288 |
| 2 | 1.408000e+09 | 0.250000000000 | 262144 |
| 3 | 2.112000e+09 | 0.125000000000 | 131072 |
| 4 | 2.816000e+09 | 0.062500000000 | 65536 |
| 5 | 3.520000e+09 | 0.031250000000 | 32768 |
| 6 | 4.224000e+09 | 0.015625000000 | 16384 |
| 7 | 4.928000e+09 | 0.007812500000 | 8192 |
| 8 | 5.632000e+09 | 0.003906250000 | 4096 |
---------------------------------------------------------------------------
Uncomment the next 4 lines of code in main which gives us this. Note that the result is returned as a dictionary.
Program Output: at_time
{'remaining_decimal': 0.7071067811865476, 'remaining_amount': 741455.2001894653}
Finally run the last 6 lines of code.
Program Output: time_to
2.2e+24 Years
2,200,000,000,000,000,016,777,216 Years
2,200,000,000,000,000,000 Million Years
2,200,000,000,000,000 Billion Years
2,200,000,000,000 Trillion Years
I hope you found this article interesting and useful. I have used radioactive decay to illustrate the concept of exponential decay but the principle has many other applications so it is useful to both understand it and be able to implement it in code.
For updates on the latest posts please follow CodeDrome on Twitter