3
$\begingroup$

For teaching purposes I am trying to create some plots of radioactive decay chains, whereby each node on the graph represents the isotope with a color corresponding to the isotopes stability and with an edge weight that corresponds to the probability of decay.

All of this information is available as collated entity data for isotopes from Wolfram:

IsotopeData["Properties"]

{"AtomicMass", "AtomicNumber", "BindingEnergy", "BranchingRatios","DaughterNuclides", "DecayEnergies", "DecayModes","DecayModeSymbols", "DecayProducts", "ExcitedStateEnergies","ExcitedStateHalfLives", "ExcitedStateLifetimes", "ExcitedStateParities", "ExcitedStateSpins", "ExcitedStateWidths", "FullSymbol", "HalfLife", "IsotopeAbundance", "Lifetime", "MagneticMoment", "MassExcess", "MassNumber", "Memberships", "Name", "NeutronNumber", "Parity", "QuadrupoleMoment", "QuantumStatistics", "Spin", "Stable", "StandardName", "Symbol", "Width"}

One example in the Documentation gives a neat way of extracting all daughter nuclides and plotting them as a graph:

DaughterNuclides[s_List] := DeleteCases[Union[Apply[Join, Map[IsotopeData[#, "DaughterNuclides"] &, DeleteCases[s, _Missing]]]], _Missing];
ReachableNuclides[s_List] := FixedPoint[Union[Join[#, DaughterNuclides[#]]] &, s];

verts = ReachableNuclides[{Entity["Isotope", "Uranium232"]}];
DaughterNuclidesQ[s1_,s2_] := (s1 =!= s2 && MemberQ[DaughterNuclides[{s1}], s2]);
RelationGraph[DaughterNuclidesQ, verts, Sequence[VertexLabels -> {Entity["Isotope", "Bismuth210"] -> Row[{Superscript["", 210], "Bi"}], Entity["Isotope", "Bismuth212"] -> 
Row[{Superscript["", 212], "Bi"}], Entity["Isotope", "Lead204"] -> Row[{Superscript["", 204], "Pb"}], Entity["Isotope", "Lead206"] -> Row[{Superscript["", 206], "Pb"}], Entity["Isotope", "Lead208"] -> Row[{Superscript["", 208], "Pb"}], Entity["Isotope", "Lead210"] -> Row[{Superscript["", 210], "Pb"}], Entity["Isotope", "Lead212"] -> Row[{Superscript["", 212], "Pb"}], Entity["Isotope", "Mercury200"] -> 
Row[{Superscript["", 200], "Hg"}], Entity["Isotope", "Mercury204"] -> 
Row[{Superscript["", 204], "Hg"}], Entity["Isotope", "Mercury206"] -> 
Row[{Superscript["", 206], "Hg"}], Entity["Isotope", "Polonium210"] -> Row[{Superscript["", 210], "Po"}], Entity["Isotope", "Polonium212"] -> Row[{Superscript["", 212], "Po"}], Entity["Isotope", "Polonium216"] -> Row[{Superscript["", 216], "Po"}], Entity["Isotope", "Radium220"] -> Row[{Superscript["", 220], "Ra"}], Entity["Isotope", "Radium224"] -> Row[{Superscript["", 224], "Ra"}], Entity["Isotope", "Radon216"] -> Row[{Superscript["", 216], "Rn"}], Entity["Isotope", "Radon220"] -> Row[{Superscript["", 220], "Rn"}], Entity["Isotope", "Thallium206"] -> Row[{Superscript["", 206], "Tl"}], Entity["Isotope", "Thallium208"] -> Row[{Superscript["", 208], "Tl"}], Entity["Isotope", "Thorium228"] -> Row[{Superscript["", 228], "Th"}], Entity["Isotope", "Uranium232"] -> Row[{Superscript["", 232], "U"}]}, PlotRangePadding -> 0.65, ImageSize -> 300, PlotTheme -> "Scientific"]]

This gives the following graph:

Graph of Decays

How can I modify this code so that the color of each node corresponds to half-life and the thickness of the lines corresponds to decay probability?

$\endgroup$

1 Answer 1

6
$\begingroup$

I struggled with this for a while, but have come up with the following solution which although slow is suitable for creating images for teaching. I use IsotopeData to get hold of half-lives and decay probabilities and then use this to style edges and to make graphics for the vertices. Perhaps somebody may find it useful (the reason I am posting). If anyone has suggestions or improvements please feel free to modify! The code is somewhat slow I think because of calling all of the data from Wolfram, but seems to work ok for my purposes.

I have used the solution from this link to get the plot colours of the scientific theme.

decayPlot[isotopename_] := Module[{DaughterNuclides, ReachableNuclides,   DaughterNuclidesQ, isoSymbol, verts, halflives, minmaxlog, shapecalc, vertshapes, branchratio, edgerules, edgestyle,colorvals},
colors = (("DefaultPlotStyle" /. (Method /. Charting`ResolvePlotTheme["Scientific", ListLinePlot])) /. Directive[x_, __] :> x); 
DaughterNuclides[s_List] := DeleteCases[Union[Apply[Join, Map[IsotopeData[#, "DaughterNuclides"] &, DeleteCases[s, _Missing]]]], _Missing];
ReachableNuclides[s_List] := FixedPoint[Union[Join[#, DaughterNuclides[#]]] &, s];
DaughterNuclidesQ[s1_,s2_] := (s1 =!= s2 && MemberQ[DaughterNuclides[{s1}], s2]);
isoSymbol[isotope_] := Row[{Superscript["", IsotopeData[isotope, "MassNumber"]],IsotopeData[isotope, "AtomicSymbol"]}];
verts = ReachableNuclides[{Entity["Isotope", isotopename]}];
halflives = (QuantityMagnitude[UnitConvert[#["HalfLife"], "s"]]) & /@verts;
minmaxlog = Log[MinMax[Select[halflives, # != Infinity &]]];
colorvals={0, 0.25, 0.5, 0.75, 1};
shapecalc[isotope_] := Graphics[{EdgeForm[Thickness[0.06]], If[QuantityMagnitude[UnitConvert[isotope["HalfLife"], "s"]] == Infinity, White, 
  ColorData["Pastel"][1 - (Log[QuantityMagnitude[UnitConvert[isotope["HalfLife"], "s"]]] - minmaxlog[[1]])/Total[minmaxlog]]],, Disk[], Black,Text[Style[isoSymbol[isotope], FontSize -> 6], {0, 0}]}, ImageSize -> 20];
vertshapes = Map[# -> shapecalc[#] &, verts];
branchratio[edge_] := QuantityMagnitude[edge[[1]]["BranchingRatios"][[Position[edge[[1]]["DaughterNuclides"], edge[[2]]][[1, 1]]]]]/100;
edgerules = EdgeRules[RelationGraph[DaughterNuclidesQ, verts]];
edgestyle = Map[(#[[1]] \[DirectedEdge] #[[2]] ->If[branchratio[#] < 0.05, {Dashed,Thickness[0.002]}, {Thickness[0.01]}]) &, edgerules];
Legended[RelationGraph[DaughterNuclidesQ, verts, VertexLabels -> None, VertexShape -> vertshapes, VertexSize -> 0.6, PlotRangePadding -> 0.65`, ImageSize -> 300, PlotTheme -> "Scientific", EdgeStyle -> edgestyle, GraphLayout -> "LayeredDigraphEmbedding"], {SwatchLegend[Join[{White}, Map[ColorData["Pastel"][#] &, colorvals]],Join[{"Stable"}, Map[Quantity[ScientificForm[Exp[minmaxlog[[2]] - #*(minmaxlog[[2]] - minmaxlog[[1]])], 3], "s"] &, colorvals]]],LineLegend[{Directive[{colors[[1]], Thickness[0.01]}],Directive[{colors[[1]], Dashed, Thickness[0.005]}]}, {"Probable Decay", "Improbable Decay"}]}]]

decayPlot["Uranium235"]

This then gives:

Decay Uranium 235

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.