I took the day and wrote up some code to generate simulated radiation events based on some previously theorized decay chains. The Event and decay chain data were taken from the Thorium series section of the https://en.wikipedia.org/wiki/Decay_chains website. This is a proof of concept for the computational logic rather than a final product, so there are things that can and should be done such as reading the decay chain data off of input files rather than hard coding it, and many more data output options should be added.
Anyway, here is a simulation of the first hour of events from a pure sample of 100,000 RA-224 nuclei:
Does the code work? It seems to work but I have not compared it to any real data. It’s certainly got advantages over handling the real deal. The algorithm used to calculate the time for each decay is likely to result in errors in the key regions at the beginnings and end of each half-life. This is especially problematic when considering the first part of the first half life, as that will be the time in which the most heat energy would be expected to be produced. This algorithm will need to be improved. The program is divided into four parts, here is the code:
First is DataExport.java
package Radioactivity;
import java.io.BufferedWriter;
import java.io.FileWriter;
import java.io.IOException;
import java.io.PrintWriter;
public class DataExport {
public static void main(String[] args) {
RA224 test = new RA224(100000);
double[][] data = test.puGetEnergyAndTime();
for (int x = 0; x < data.length;x++) {
scribe(Double.toString(data[x][0]) + " " + Double.toString(data[x][1]));
}
}
private static void scribe(String textline) {
PrintWriter writer = null;
BufferedWriter bw = null;
FileWriter fw = null;
try {
fw = new FileWriter("/home/user/workspace/Radioactivity/output/out.txt",true);
bw = new BufferedWriter(fw);
writer = new PrintWriter(bw);
writer.println(textline);
} catch (IOException ex) {
System.out.println(ex);
} finally {
try {writer.close();} catch (Exception ex) {}
}
}
}
Next is RA224.java which contains the decay chain logic
package Radioactivity;
public class RA224 {
// a group of partNum ra224 nuclei whose radioactive Events are calculated for maxYears
/* Variable and Function Nomenclature Prescipts:
* pv = private
* pu = public
* pvs = private static
* pus = public static
* pvsf = private static final
* pusf = public static final
*/
private int pvPartNum = 0;
private Event[] pvEvents;
public RA224(int num) {
if (num > 0) {
pvPartNum = num;
pvEvents = new Event[pvPartNum*6];
eventCalc();
} else {
System.out.println("Construction of this (RA224) failed because the (num) supplied was not greater than zero");
pvPartNum = 0;
}
}
private void eventCalc() {
for(int x = 0; x<pvPartNum; x++) {
//ra224 to rn220
//type = alpha
//energy = 5.789 MeV
//halflife = 3.6319 days
String type = "alpha";
double energy = 5.789;
double halflife = 24*60*60*3.6319;
pvEvents[x] = new Event(halflife,energy,type);
}
for(int x = 0; x<pvPartNum; x++) {
//rn220 to po216
//type = alpha
//energy = 6.404 MeV
//halflife = 55.6 seconds
String type = "alpha";
double energy = 6.404;
double halflife = 55.6;
pvEvents[x+pvPartNum] = new Event(halflife,energy,type,pvEvents[x].puGetTime());
}
for(int x = 0; x<pvPartNum; x++) {
//po216 to pb212
//type = alpha
//energy = 6.906 MeV
//halflife = 0.145 seconds
String type = "alpha";
double energy = 6.906;
double halflife = 0.145;
pvEvents[x+2*pvPartNum] = new Event(halflife,energy,type,pvEvents[x+pvPartNum].puGetTime());
}
for(int x = 0; x<pvPartNum; x++) {
//pb212 to bi212
//type = beta-
//energy = 0.570
//halflife = 10.64 hours
String type = "beta-";
double energy = 0.570;
double halflife = 60*60*10.64;
pvEvents[x+3*pvPartNum] = new Event(halflife,energy,type,pvEvents[x+2*pvPartNum].puGetTime());
}
int numPO212 = Math.toIntExact(Math.round(Math.random()*pvPartNum));
int numTI208 = pvPartNum - numPO212;
//bi212 to either po212 or ti208
//po212 type = beta-
//ti208 type = alpha
//po212 energy = 2.252 MeV
//ti208 energy = 6.208 MeV
//po212 halflife = 60.55 min
//ti208 halflife = 60.55 min
for(int x=0;x<numPO212;x++){
String type = "beta-";
double energy = 2.252;
double halflife = 60*60.55;
pvEvents[x+4*pvPartNum] = new Event(halflife,energy,type,pvEvents[x+3*pvPartNum].puGetTime());
}
for(int x=0;x<numTI208;x++){
String type = "alpha";
double energy = 6.208;
double halflife = 60*60.55;
pvEvents[x+4*pvPartNum+numPO212] = new Event(halflife,energy,type,pvEvents[x+3*pvPartNum+numPO212].puGetTime());
}
for(int x = 0; x<numPO212; x++) {
//po212 to pb208
//type = alpha
//energy = 8.955 MeV
//halflife = 299 ns
String type = "alpha";
double energy = 8.955;
double halflife = 299*Math.pow(10,-9);
pvEvents[x+5*pvPartNum] = new Event(halflife,energy,type,pvEvents[x+4*pvPartNum].puGetTime());
}
for(int x = 0; x<numTI208; x++) {
//po212 to pb208
//type = beta-
//energy = 4.999 MeV
//halflife = 3.053 min
String type = "beta-";
double energy = 4.999;
double halflife = 60*3.053;
pvEvents[x+5*pvPartNum+numPO212] = new Event(halflife,energy,type,pvEvents[x+4*pvPartNum+numPO212].puGetTime());
}
}
public double[][] puGetEnergyAndTime() {
double returned[][] = new double[pvPartNum*6][2];
for (int x = 0;x<pvPartNum*6;x++) {
returned[x][0]=pvEvents[x].puGetTime();
returned[x][1]=pvEvents[x].puGetEnergy();
}
return returned;
}
}
Next is TH232.java with decay logic for that nuclei (not used)
package Radioactivity;
public class TH232 {
// a group of partNum th232 nuclei whose radioactive Events are calculated for maxYears
/* Variable and Function Nomenclature Prescipts:
* pv = private
* pu = public
* pvs = private static
* pus = public static
* pvsf = private static final
* pusf = public static final
*/
private int pvPartNum = 0;
private Event[] pvEvents;
public TH232(int num) {
if (num > 0) {
pvPartNum = num;
pvEvents = new Event[pvPartNum*10];
eventCalc();
} else {
System.out.println("Construction of this (TH232) failed because the (num) supplied was not greater than zero");
pvPartNum = 0;
}
}
private void eventCalc() {
for(int x = 0; x<pvPartNum; x++) {
//th232 to ra228
//type = alpha
//energy = 4.081 MeV
//halflife = 1.405 * 10^10 years
String type = "alpha";
double energy = 4.081;
double halflife = 365*24*60*60*1.405*Math.pow(10,10);
pvEvents[x] = new Event(halflife,energy,type);
}
for(int x = 0; x<pvPartNum; x++) {
//ra228 to ac228
//type = beta-
//energy = 0.046 MeV
//halflife = 5.75 years
String type = "beta-";
double energy = 0.046;
double halflife = 365*24*60*60*5.75;
pvEvents[x+pvPartNum] = new Event(halflife,energy,type,pvEvents[x].puGetTime());
}
for(int x = 0; x<pvPartNum; x++) {
//ac228 to th228
//type = beta-
//energy = 2.124 MeV
//halflife = 6.35 hours
String type = "beta-";
double energy = 2.124;
double halflife = 60*60*6.35;
pvEvents[x+2*pvPartNum] = new Event(halflife,energy,type,pvEvents[x+pvPartNum].puGetTime());
}
for(int x = 0; x<pvPartNum; x++) {
//th228 to ra224
//type = alpha
//energy = 5.520 MeV
//halflife = 1.9116 years
String type = "alpha";
double energy = 5.520;
double halflife = 365*24*60*60*1.9116;
pvEvents[x+3*pvPartNum] = new Event(halflife,energy,type,pvEvents[x+2*pvPartNum].puGetTime());
}
for(int x = 0; x<pvPartNum; x++) {
//ra224 to rn220
//type = alpha
//energy = 5.789 MeV
//halflife = 3.6319 days
String type = "alpha";
double energy = 5.789;
double halflife = 24*60*60*3.6319;
pvEvents[x+4*pvPartNum] = new Event(halflife,energy,type,pvEvents[x+3*pvPartNum].puGetTime());
}
for(int x = 0; x<pvPartNum; x++) {
//rn220 to po216
//type = alpha
//energy = 6.404 MeV
//halflife = 55.6 seconds
String type = "alpha";
double energy = 6.404;
double halflife = 55.6;
pvEvents[x+5*pvPartNum] = new Event(halflife,energy,type,pvEvents[x+4*pvPartNum].puGetTime());
}
for(int x = 0; x<pvPartNum; x++) {
//po216 to pb212
//type = alpha
//energy = 6.906 MeV
//halflife = 0.145 seconds
String type = "alpha";
double energy = 6.906;
double halflife = 0.145;
pvEvents[x+6*pvPartNum] = new Event(halflife,energy,type,pvEvents[x+5*pvPartNum].puGetTime());
}
for(int x = 0; x<pvPartNum; x++) {
//pb212 to bi212
//type = beta-
//energy = 0.570
//halflife = 10.64 hours
String type = "beta-";
double energy = 0.570;
double halflife = 60*60*10.64;
pvEvents[x+7*pvPartNum] = new Event(halflife,energy,type,pvEvents[x+6*pvPartNum].puGetTime());
}
int numPO212 = Math.toIntExact(Math.round(Math.random()*pvPartNum));
int numTI208 = pvPartNum - numPO212;
//bi212 to either po212 or ti208
//po212 type = beta-
//ti208 type = alpha
//po212 energy = 2.252 MeV
//ti208 energy = 6.208 MeV
//po212 halflife = 60.55 min
//ti208 halflife = 60.55 min
for(int x=0;x<numPO212;x++){
String type = "beta-";
double energy = 2.252;
double halflife = 60*60.55;
pvEvents[x+8*pvPartNum] = new Event(halflife,energy,type,pvEvents[x+7*pvPartNum].puGetTime());
}
for(int x=0;x<numTI208;x++){
String type = "alpha";
double energy = 6.208;
double halflife = 60*60.55;
pvEvents[x+8*pvPartNum+numPO212] = new Event(halflife,energy,type,pvEvents[x+7*pvPartNum+numPO212].puGetTime());
}
for(int x = 0; x<numPO212; x++) {
//po212 to pb208
//type = alpha
//energy = 8.955 MeV
//halflife = 299 ns
String type = "alpha";
double energy = 8.955;
double halflife = 299*Math.pow(10,-9);
pvEvents[x+9*pvPartNum] = new Event(halflife,energy,type,pvEvents[x+8*pvPartNum].puGetTime());
}
for(int x = 0; x<numTI208; x++) {
//po212 to pb208
//type = beta-
//energy = 4.999 MeV
//halflife = 3.053 min
String type = "beta-";
double energy = 4.999;
double halflife = 60*3.053;
pvEvents[x+9*pvPartNum+numPO212] = new Event(halflife,energy,type,pvEvents[x+8*pvPartNum+numPO212].puGetTime());
}
}
public double[][] puGetEnergyAndTime() {
double returned[][] = new double[pvPartNum*10][2];
for (int x = 0;x<pvPartNum*10;x++) {
returned[x][0]=pvEvents[x].puGetTime();
returned[x][1]=pvEvents[x].puGetEnergy();
}
return returned;
}
}
Finally is the Event.java class that actually generates when the random events occur based on half life.
package Radioactivity;
public class Event {
// An Event calculates and contains one radiative event
/* Variable and Function Nomenclature prescripts:
* pv = private
* pu = public
* pvs = private static
* pus = public static
* pvsf = private static final
* pusf = public static final
*/
private double pvTime = 0.0; //calculated time of the event referenced from arbitrary start point in seconds
private double pvEnergy = 0.0; //user supplied event energy in MeV
private String pvType; //user supplied event type (alpha, beta, gamma, neutron)
private double pvHalfLife =0.0; //user supplied event half-life in seconds
private int pvMaxHalfLives = 20;
private double pvTimeOffset = 0.0; //user supplied offset time in seconds
public Event(double halfLife, double energy, String type) {
if (energy > 0) {
pvEnergy = energy;
} else {
System.out.println("(Event) construction failed because (Event) input (energy) must be a positive number and greater than zero");
}
pvType = type;
if (halfLife > 0) {
pvHalfLife = halfLife;
} else {
System.out.println("(Event) construction failed because input (halflife) must be a positive number and greater than zero");
}
if (pvHalfLife > 0) {
pvTime = pvCalcTime();
} else {
System.out.println("(Event) construction failed because (pvTime) cannot be calculated without a defined (halflife)");
}
pvMaxHalfLives = 20;
pvTimeOffset = 0.0;
}
public Event(double halfLife, double energy, String type,double timeOffset) {
if (timeOffset > 0) {
pvTimeOffset = timeOffset;
} else {
System.out.println("(Event) construction failed because (timeOffset) input must be a positive number and greater than zero");
}
if (energy > 0) {
pvEnergy = energy;
} else {
System.out.println("(Event) construction failed because (Event) input (energy) must be a positive number and greater than zero");
}
pvType = type;
if (halfLife > 0) {
pvHalfLife = halfLife;
} else {
System.out.println("(Event) construction failed because input (halflife) must be a positive number and greater than zero");
}
if (pvHalfLife > 0) {
pvTime = pvCalcTime();
} else {
System.out.println("(Event) construction failed because (pvTime) cannot be calculated without a defined (halflife)");
}
pvMaxHalfLives = 20;
}
public Event(double halfLife, double energy, String type, int maxHalfLives) {
if (energy > 0) {
pvEnergy = energy;
} else {
System.out.println("(Event) construction failed because (Event) input (energy) must be a positive number and greater than zero");
}
pvType = type;
if (halfLife > 0) {
pvHalfLife = halfLife;
} else {
System.out.println("(Event) construction failed because input (halflife) must be a positive number and greater than zero");
}
if (maxHalfLives > 0) {
pvMaxHalfLives = maxHalfLives;
} else {
System.out.println("Using default (pvMaxHalfLives) of 20 because the (maxHalfLives) supplied to the constructor was not a positive number and greater than zero");
pvMaxHalfLives = 20;
}
if (pvHalfLife > 0) {
pvTime = pvCalcTime();
} else {
System.out.println("(Event) construction failed because (pvTime) cannot be calculated without a defined (halflife)");
}
pvTimeOffset = 0.0;
}
public Event(double halfLife, double energy, String type,double timeOffset, int maxHalfLives) {
if (timeOffset > 0) {
pvTimeOffset = timeOffset;
} else {
System.out.println("(Event) construction failed because (timeOffset) input must be a positive number and greater than zero");
}
if (energy > 0) {
pvEnergy = energy;
} else {
System.out.println("(Event) construction failed because (Event) input (energy) must be a positive number and greater than zero");
}
pvType = type;
if (halfLife > 0) {
pvHalfLife = halfLife;
} else {
System.out.println("(Event) construction failed because input (halflife) must be a positive number and greater than zero");
}
if (maxHalfLives > 0) {
pvMaxHalfLives = maxHalfLives;
} else {
System.out.println("Using default (pvMaxHalfLives) of 20 because the (maxHalfLives) supplied to the constructor was not a positive number and greater than zero");
pvMaxHalfLives = 20;
}
if (pvHalfLife > 0) {
pvTime = pvCalcTime();
} else {
System.out.println("(Event) construction failed because (pvTime) cannot be calculated without a defined (halflife)");
}
}
private double pvCalcTime() {
double seed1 = Math.random();
double seed2 = Math.random();
double weightedLife = Math.pow(2,pvMaxHalfLives);
double testWeightedLife = seed1 * weightedLife;
double t = 0;
int x = 0;
while (t == 0) {
if (x == (pvMaxHalfLives-1)) {
t = (x+seed2) * pvHalfLife;
} else if (testWeightedLife <= weightedLife*(1+(x*2))/Math.pow(2,x+1)){
t = (x+seed2) * pvHalfLife;
}
x+=1;
}
return (t+pvTimeOffset);
}
public double puGetTime() {
if (pvTime > 0) {
return pvTime;
} else {
System.out.println("(puGetTime) failed because this (Event) does not have a calculated (pvTime)");
return 0;
}
}
public double puGetHalfLife() {
if (pvHalfLife > 0) {
return pvHalfLife;
} else {
System.out.println("(puGetHalfLife) failed because this (Event) does not have a proper (pvHalflife) assigned");
return 0;
}
}
public String puGetType() {
return pvType;
}
public double puGetEnergy() {
if (pvEnergy > 0) {
return pvEnergy;
} else {
System.out.println("(puGetEnergy) failed because this (Event) does not have a proper (pvEnergy) assigned");
return 0;
}
}
public void puSetHalfLife(double halflife) {
if (halflife > 0) {
pvHalfLife = halflife;
} else {
System.out.println("(puSetHalfLife) failed because the (halflife) provided is not greater than zero");
pvHalfLife = 0;
}
}
public void puSetHalfLife(double halflife, int maxhalflives) {
if (halflife > 0) {
pvHalfLife = halflife;
if (maxhalflives > 0) {
pvMaxHalfLives = maxhalflives;
pvTime = pvCalcTime();
} else {
System.out.println("(puSetHalfLife) failed because the (maxhalflives) provided is not greater than zero");
pvMaxHalfLives = 0;
pvTime = 0.0;
}
} else {
System.out.println("(puSetHalfLife) failed because the (halflife) provided is not greater than zero");
pvHalfLife = 0;
}
}
public void puSetType(String newtype) {
pvType = newtype;
}
public void puSetEnergy(double energy) {
if (energy > 0) {
pvEnergy = energy;
} else {
System.out.println("(puSetEnergy) failed because the (energy) provided is not greater than zero");
pvEnergy = 0;
}
}
}
