Create masses.h using online isotope data from wolfram.com

Author : Martin Torheyden (mtorhey at vt period edu)

The sources of data are described at http://reference.wolfram.com/mathematica/note/IsotopeDataSourceInformation.html

The sources of data as listed on March 3, 2008:

Atomic Mass Data Center."NUBASE." 2003. »
Firestone, R.B."The Berkeley Laboratory Isotopes Project's: Exploring the Table of Isotopes." 2000. »
Raghavan, P."Table of Nuclear Moments." Atomic Data and Nuclear Data Tables 42, no.2 (1989) : 189 - 291.
Sansonetti, J.E.and W.C.Martin."NIST Handbook of Basic Atomic Spectroscopic Data." 2005. »
United States National Institute of Standards and Technology."Atomic Weights and Isotopic Compositions Elements." 2005. »
United States National Nuclear Data Center, Brookhaven National Laboratory."Nuclear Wallet Cards" 2007. »
United States National Nuclear Data Center, Brookhaven National Laboratory."NuDat 2.3." 2007. »

extractPsiSymbol[isotopesymbol_]:=Module[
{chemsymb,massnumber,psisymbol},
chemsymb=ToUpperCase[isotopesymbol[[1,2]]];
massnumber=isotopesymbol[[1,1,2]];
psisymbol=StringJoin[chemsymb,ToString[massnumber]];
Return[psisymbol];
];

getIsotopeData[atomname_]:=Module[
{atomsymbol,isotopename,i,isotopesymbol,psisymbol,isotopemass,isotopeabundance,isotopeHalfLife,isotopedata},
atomsymbol=ToUpperCase[ElementData[atomname,"Symbol"]];
(*Print["atomsymbol=",atomsymbol];*)
isotopename=IsotopeData[atomname];
isotopedata={atomsymbol,{}};
For[i=1,i≤Length[isotopename],i++,
isotopesymbol=IsotopeData[isotopename[[i]],"Symbol"];
(*Print["isotopesymbol=",isotopesymbol];*)
psisymbol=extractPsiSymbol[isotopesymbol];
(*Print["psisymbol=",psisymbol];*)
isotopemass=IsotopeData[isotopename[[i]],"AtomicMass"];
(*Print["isotopemass=",isotopemass];*)
isotopeabundance=IsotopeData[isotopename[[i]],"IsotopeAbundance"];
(*Print["isotopeabundance=",isotopeabundance];*)
isotopeHalfLife=IsotopeData[isotopename[[i]],"HalfLife"];
(*Print["isotopeHalfLife=",isotopeHalfLife];*)
If[NumberQ[isotopemass]&&NumberQ[isotopeabundance],
isotopedata[[2]]=Append[isotopedata[[2]],{psisymbol,isotopemass,isotopeabundance,isotopeHalfLife}];
];
];
(*Print["isotopedata=",isotopedata];*)
Return[isotopedata];
];

getIsotopes[]:=Module[
{atomname,atomsymbol,i,isotopes},
atomname=ElementData[];
(*Print["atomname=",atomname];*)
isotopes={};
For[i=1,i≤Length[atomname],i++,
(*getIsotopeData[atomname[[i]]];*)
isotopes=Append[isotopes,getIsotopeData[atomname[[i]]]];
];
Return[isotopes];
];

InsertListintoListofLists[OldlistofLists_List,ToInsert_List,pos_]:=Module[
{NewlistofLists,i},
NewlistofLists={};
For[i=1,i≤pos-1,i++,
NewlistofLists=Append[NewlistofLists,OldlistofLists[[i]]];
];
NewlistofLists=Append[NewlistofLists,ToInsert];
For[i=pos,i≤Length[OldlistofLists],i++,
NewlistofLists=Append[NewlistofLists,OldlistofLists[[i]]];
];
Return[NewlistofLists];
];

isstableIsotope[isotope_List]:=Module[
{isstable,i},
isstable=False;
For[i=1,i≤Length[isotope[[2]]],i++,
If[isotope[[2,i,4]]===Infinity,
isstable=True;
];
];
Return[isstable];
];

isnaturallyoccurringIsotope[isotope_List]:=Module[
{sum,isnaturallyoccurring,i},
sum=0.0;
For[i=1,i≤Length[isotope[[2]]],i++,
sum=sum+isotope[[2,i,3]];
];
(*Print["sum=",sum];*)
If[sum>0.1,
isnaturallyoccurring=True,
isnaturallyoccurring=False;
];
Return[isnaturallyoccurring];
];

getMostStableIsotope[isotope_List]:=Module[
{moststableindex,i},
moststableindex=1;
For[i=2,i≤Length[isotope[[2]]],i++,
If[isotope[[2,i,4]]>isotope[[2,moststableindex,4]],
moststableindex=i;
];
];
Return[isotope[[2,moststableindex]]];
];

getMostAbundantIsotope[isotope_List]:=Module[
{mostabundantindex,i},
mostabundantindex=1;
For[i=2,i≤Length[isotope[[2]]],i++,
If[isotope[[2,i,3]]>isotope[[2,mostabundantindex,3]],
mostabundantindex=i;
];
];
Return[isotope[[2,mostabundantindex]]];
];

getMostCommonIsotope[isotope_List]:=Module[
{mostcommonisotope},
If[isnaturallyoccurringIsotope[isotope],
mostcommonisotope=getMostAbundantIsotope[isotope],
mostcommonisotope=getMostStableIsotope[isotope];
];
Return[mostcommonisotope];
];

getTotalNumberofIsotopes[isotopes_List]:=Module[
{i,numofisotopes},
numofisotopes=0;
For[i=1,i≤Length[isotopes],i++,
numofisotopes=numofisotopes+Length[isotopes[[i,2]]];
(*Print[isotopes[[i,1]]," Number of isotopes =",Length[isotopes[[i,2]]]];*)
];
Return[numofisotopes];
];

getMostCommonIsotopeMassList[isotopes_List]:=Module[
{masslist,i},
masslist={};
For[i=1, i≤Length[isotopes],i++,
masslist=Append[masslist,{isotopes[[i,1]],getMostCommonIsotope[isotopes[[i]]][[2]]}];
];
(*Print[masslist];*)
Return[masslist];
];

getIsotopeMassList[isotopes_List]:=Module[
{masslist,i,j,mostcommonisotope},
masslist={};
For[i=1, i≤Length[isotopes],i++,
mostcommonisotope=getMostCommonIsotope[isotopes[[i]]];
masslist=Append[masslist,{isotopes[[i,1]],mostcommonisotope[[2]]}];
For[j=1,j≤Length[isotopes[[i,2]]],j++,
masslist=Append[masslist,{isotopes[[i,2,j,1]],isotopes[[i,2,j,2]]}];
];
];
(*Print[masslist];*)
(*Print["Length[masslist]=",Length[masslist]];*)
Return[masslist];
];

createIsotopeAliases[isotopes_List,aliases_List]:=Module[
{i,j,k,subst,isot,isotopesnew},
isotopesnew=isotopes;
For[i=1,i≤Length[aliases],i++,  (* loop over different aliases *)
For[j=1,j≤Length[isotopesnew], j++,  (* loop over different atomtypes *)
For[k=1,k≤Length[isotopesnew[[j,2]]],k++, (* loop over different isotopesnew *)
(*Print[isotopesnew[[j,2,k]]];*)
If[aliases[[i,1]]===isotopesnew[[j,2,k,1]],
isot=isotopesnew[[j]];
(*Print["isot=",isot];*)
subst={aliases[[i,2]],isotopesnew[[j,2,k,2]],isotopesnew[[j,2,k,3]],isotopesnew[[j,2,k,4]]};
(*Print["subst=",subst];*)
  isot[[2]]=InsertListintoListofLists[isot[[2]],subst,k+1];
(*Print["isot=",isot];*)

isotopesnew=Drop[isotopesnew,{j}];
isotopesnew=InsertListintoListofLists[isotopesnew,isot,j];
(*Print["isotopesnew[[j]]=",isotopesnew[[j]]];*)
j=j+1;
];

];
];
];
Return[isotopesnew];
];

printPsiElementData[]:=Module[
{isotopes,mostcommonisotopemasslist,isotopemasslist,mostcommonisotopemasslistLabels,mostcommonisotopemasslistValues,isotopemasslistLabels,isotopemasslistValues,i,j,labelstoberemoved},
isotopes=getIsotopes[];
isotopes=createIsotopeAliases[isotopes,{{"H2","D"},{"H3","T"}}];
isotopes=InsertListintoListofLists[isotopes,{"X",{{"X0",0.0000,1.0,Infinity}}},1];
mostcommonisotopemasslist=getMostCommonIsotopeMassList[isotopes];
(*Print["mostcommonisotopemasslist=",mostcommonisotopemasslist];*)
isotopemasslist=getIsotopeMassList[isotopes];
(* Isotope X and X0 should not be in the isotopemasslist *)
labelstoberemoved={"X","X0"};
Print["Length before X removal = ",Length[isotopemasslist]];
For[j=1,j≤Length[labelstoberemoved], j++,
For[i=1,i≤Length[isotopemasslist],i++,
If[isotopemasslist[[i,1]]===labelstoberemoved[[j]],
isotopemasslist=Drop[isotopemasslist,{i}];
];
];
];
Print["Length after X removal = ",Length[isotopemasslist]];
(*Print["isotopemasslist=",isotopemasslist];*)

(* Now all the necessary data are in mostcommonisotopemasslist and isotopemasslist. Printing out. *)
mostcommonisotopemasslistLabels={};
mostcommonisotopemasslistValues={};
For[i=1,i≤Length[mostcommonisotopemasslist],i++,
mostcommonisotopemasslistLabels=Append[mostcommonisotopemasslistLabels,mostcommonisotopemasslist[[i,1]]];
mostcommonisotopemasslistLabels[[i]]=StringJoin["\"",mostcommonisotopemasslistLabels[[i]],"\""];
(*StringJoin[*)
mostcommonisotopemasslistValues=Append[mostcommonisotopemasslistValues,mostcommonisotopemasslist[[i,2]]];
];
Print["#define LAST_ATOMIC_INDEX (",Length[mostcommonisotopemasslist]-1,")"];
Print["static const char *atomic_labels[] =",mostcommonisotopemasslistLabels];
Print["static double an2masses[] =",mostcommonisotopemasslistValues];

isotopemasslistLabels={};
isotopemasslistValues={};
For[i=1,i≤Length[isotopemasslist],i++,
isotopemasslistLabels=Append[isotopemasslistLabels,isotopemasslist[[i,1]]];
isotopemasslistLabels[[i]]=StringJoin["\"",isotopemasslistLabels[[i]],"\""];
isotopemasslistValues=Append[isotopemasslistValues,isotopemasslist[[i,2]]];
];
Print["#define LAST_MASS_INDEX (",Length[isotopemasslist]-1,")"];
Print["static const char *mass_labels[] =",isotopemasslistLabels];
Print["static double atomic_masses[] =",isotopemasslistValues];
Print[""];
Print["Length[mostcommonisotopemasslist]=",Length[mostcommonisotopemasslist]];
Print["Length[mostcommonisotopemasslistLabels]=",Length[mostcommonisotopemasslistLabels]];
Print["Length[mostcommonisotopemasslistValues]=",Length[mostcommonisotopemasslistValues]];
Print["Length[isotopemasslist]=",Length[isotopemasslist]];
Print["Length[isotopemasslistLabels]=",Length[isotopemasslistLabels]];
Print["Length[isotopemasslistValues]=",Length[isotopemasslistValues]];
];

printPsiElementData[];