# -*- coding: iso-8859-15 -*- # This program is public domain # Author: Paul Kienzle # Based on spreadsheet by Les Slaback (1998). r""" Calculate expected neutron activation from time spent in beam line. Accounts for burnup and 2n, g production. This is only a gross estimate. Many effects are not taken into account, such as self-shielding in the sample and secondary activation from the decay products. **Introduction to neutron activation terminology** See a text!! These are just a few notes to refresh the memory of those who already know the topic. Reactions: The most common reaction is the simple addition of a neutron, thereby increasing its atomic number by one but leaving it as the same element. Not all activation products are radioactive. Those that are not are not included in this database. Beware that there are production chains that depend upon these intermediate products. This database does not address those more complicated processes. There are exceptions to the simple addition process such as the n,alpha reaction of Li-6. These then result in a different element as the activation product. These are identified in the reaction column. Radioactive products also have the potential of undergoing a neutron reaction. This is accounted for in this database for selected products (generally those that result in significant half-life products). These reactions and products should only be important at very high fluences, i.e., >1E16 n/cm2. Neutron energy: The majority of the reactions are initiated by thermal neutrons. As a practical matter the number of thermal neutrons is usually measured with a cadmium filter. This excludes the neutrons above the cadmium absorption threshold, called the cadmium cutoff. For most materials the resulting measured 'thermal' neutron fluence is adequate to determine the thermal neutron activation. A few materials have large cross-sections for the neutrons above the cadmium cutoff. Hence a 'cadmium ratio' is needed to predict the number of neutrons present above this cutoff, as seen by the specific reaction of interest. For the same neutron spectrum two different elements will have different cadmium ratios. Similarly, for the same reaction but different neutron environments the ratio will also vary. Copper is a good material on which to predict the thermal-to-epithermal fluence ratio since its two cross-sections are about equal. If you use the ratio of another reaction, correct the activity ratio of the nuclide by the cross-section ratio in order to derive the extimated fluence ratio. Take care. By definition, the cadmium ratio is an activity ratio, not a fluence ratio. [Note: This ratio is divided into the specfied thermal neutron fluence rate to get the epithermal rate.] Fast neutron reactions: There are also those reactions that depend on higher energy neutrons. These can be of a great variety: n,gamma; n,p; n,d; n,alpha; n,triton; etc. Any they are highly variable in how high the neutron energy must be to initiate the reaction. In general the cross- sections are relatively small, e.g., 10's of millibarns and less as compared to 1000's of barns in some cases for thermal reactions. Plus the fast neutron fluences are much smaller so that in general the amount of the activation product is much less than for thermal production processes. Barn: As they say- just how big is the broad side of that barn? For neutrons it is 1 E-24 cm2, e.g. 0.000000000000000000000001 square cm. Millibarns are of course for the _______ (leprechauns?). Calulation: (number of atoms)*(cross-section)*(no. of neutrons)*(decay correction)*C.F. Number of atoms: mass*avagadro's number*isotopic abundance/mole.weight where the mass is that of the target element, not the whole sample. Decay correction: Since the radioactive product decays while it is being made one must correct for this using (1-exp(-0.693*t/hlflf)) where t is the exposure time and hlflf is the product halflife C.F: includes unit conversion factors, e.g., convert to microcuries Example:: >>> from periodictable import activation >>> env = activation.ActivationEnvironment(fluence=1e5, Cd_ratio=70, fast_ratio=50, location="BT-2") >>> sample = activation.Sample("Co30Fe70", 10) >>> sample.calculate_activation(env, exposure=10, rest_times=[0, 1, 24, 360]) >>> sample.show_table() ----------------- activity (uCi) ------------------ isotope product reaction half-life 0 hrs 1 hrs 24 hrs 360 hrs -------- --------- -------- ---------- ------------ ------------ ------------ ------------ Co-59 Co-60 act 5.272 y 0.0004959 0.0004959 0.0004958 0.0004933 Co-59 Co-60m+ act 10.5 m 1.664 0.03169 --- --- -------- --------- -------- ---------- ------------ ------------ ------------ ------------ total 1.664 0.0322 0.0005084 0.0005049 -------- --------- -------- ---------- ------------ ------------ ------------ ------------ >>> print("%.3f"%sample.decay_time(0.001)) # number of hours to reach 1 nCi 2.053 The default rest times used above show the sample activity at the end of neutron activation and after 1 hour, 1 day, and 15 days. The neutron activation table, *activation.dat*,\ [#Shleien1998]_ contains details about the individual isotopes, with interaction cross sections taken from from IAEA-273\ [#IAEA1987]_. Activation can be run from the command line using:: $ python -m periodictable.activation FORMULA where FORMULA is the chemical formula for the material. .. [#Shleien1998] Shleien, B., Slaback, L.A., Birky, B.K., 1998. Handbook of health physics and radiological health. Williams & Wilkins, Baltimore. .. [#IAEA1987] IAEA (1987) Handbook on Nuclear Activation Data. TR 273 (International Atomic Energy Agency, Vienna, Austria, 1987). http://cds.cern.ch/record/111089/files/IAEA-TR-273.pdf Code original developed for spreadsheet by Les Slaback of NIST. """ # Comments on ACT_2N_X.xls from ACT_CALC.TXT # # Fast neutron cross-sections added 2/1/94 # # One database error detected and corrected at this time. # # Corrections made 3/4/94 # # 1. Halflife for entry 21 entered. (previous entry blank) # 2. Cross-section for entry 311 entered (previous entry blank) # 3. 2 beta mode entries - added parent lambda's # 4. Changed range of look up function in 2 columns to include last row. # 5. With 9E15 y halflife the calculation always returned 0 microcuries! # The 9E15 exceeds a math precision resulting in zero being assigned to the # the function 1-e(-x). For any value of x<1E-10 the approximation x+x^2/2 # is more accurate, and for values of X<1E-16 this approximation returns a # non-zero value while 1-e(-x)=0. # Appropriate changes have been made, but in fact only two entries are # affected by this, one with a 12 Ty halflife and one with a 9000Ty t1/2. # The equations for "b" mode production were not changed because they are # more complex, and none of the halflives currently in the database related # to this mode are a problem. Take care if you add more with very long halflives. # [PAK: These have since been updated to use expm1()] # 6. The cross-section for the reaction Sr-88(n,gamma)Sr-89 was reduced from # .058 to .0058 b. The entry of .058 in IAEA273 is in error, based on a # number of other references (including IAEA156). # 7. The unit for the halflife of Pm-151 coorected from 'm' to 'h' # # Changes made in April 1994 # # 1. Burnup cross-sections added to the database and nuclides produced by # two neutron additions (2n,gamma) have been added. The activation equations # have been changed to account for burnup. This does not become significant # until exposures exceed 1e17 n/cm2 (e.g., 1e10 n/cm2/sec*1000 hrs), even for # those with very large cross-sections. Note that 'burnup' can be viewed as # loss of the intended n,gamma product or as a 2n,gamma production mode. Both # effects are included in the database and equations. # 2.ACT_CALC.WQ1 does not have the burnup equations or related cross-section data. # ACT_2N.WQ1 has the added equations and data, but requires manual entry of the # cross-section database index numbers. # ACT_2N_X.WQ1 allows direct entry of the chemical element to retrieve ALL # entries from the database for that element # Results from both have been compared to assure that the new spreadsheet # is correct. Also the Au-197(2n,g)Au-199 reaction has been checked against # a textbook example (Friedlander,Kennedy,Miller:Nuclear Chemistry). # 3.An educational note related to this addition: # Computing the equation [exp(-x)-exp(-y)] in that syntax is better than using # [exp(-x)*(1-exp(x-y))]. The latter format blows up when large # values of X and Y are encountered due to computational limitations for these # functions in a PC. But this problem was only encountered in excess of # 1E22 n/cm2. # 4. 61 (2n,g) reactions have been added. You will note many more # radionuclides with secondary capture cross-sections. Many of these go to # stable nuclides (particularly the larger x-section ones - logically enough). # The isomer-ground state pairs that have 2n modes to the same resulting # nuclide are treated one of several ways. # - if the obvious dominate pathway is only through the ground state then # only production via that mode is included. # - if the combination of parent halflives and production cross-sections # make it unclear as to which is the dominant production mode then both # are calculated. Beware, these are not necessarily additive values. # Sometimes the ground state also reflects production via the isomer. # See the specific notes for any particular reaction. # 5.None of the entries for production via 'b' (beta decay ingrowth from a # neutron induced parent) account for burnup. Those with significant burnup # cross-sections have a specific note reminding of this. This is an issue # only at high fluence rates. # 6.Note that the database entries have different meanings for 'b' and '2n' # production modes. For instance, the first set of cross-sections (proceeding # horizontally for a particular database entry) is not the cross-section of # the 2n reaction, but that to produce the parent. The second set of cross- # sections which are the burnup cross-sections for n,g; n,p; n,alpha; etc. # reactions are the production cross-sections for the 2n reactions. # 7.If you want to determine how much burnup is occurring do one of the # following: # - do the calculation separately in ACT_CALC and ACT_2N. # - do the calculation at a low fluence, e.g., 1E7, prorate to the fluence # of interest, and compare to the result calculated directly. Make # certain the same exposure time is used. You cannot prorate this # parameter. # # Additions/changes made July 1997 # # The following n,p and n,alpha reactions have both a thermal cross-section (as # per IAEA273) and a fast cross-section. In all cases the database entry was for # the fast cross-section but indicated it was a thermal reaction. That has been # corrected (and verified against IAEA156) and a second entry made for the thermal # reaction. Despite the entry in IAEA 273 there is some question as to whether # the thermal induced reaction is energetically possible. Dick Lindstrom's # calculation shows that the coulomb barrier should prevent the thermal reaction # from being possible. # # The entries changed, and added, are for the following: # # 35Cl (n,p) # 35Cl (n,alpha) # 33S (n,p) # 39K (n,alpha) # 40Ca (n,p) # 58Ni (n,alpha) from math import exp, log, expm1 import os from .formulas import formula as build_formula from . import core LN2 = log(2) def table_abundance(iso): """ Isotopic abundance in % from the periodic table package. """ return iso.abundance def IAEA1987_isotopic_abundance(iso): """ Isotopic abundance in % from the IAEA, as provided in the activation.dat table. Note: this will return an abundance of 0 if there is no neutron activation for the isotope even though for isotopes such as H[1], the natural abundance may in fact be rather large. IAEA 273: Handbook on Nuclear Activation Data, 1987. """ try: return iso.neutron_activation[0].abundance except AttributeError: return 0 class Sample: """ Sample properties. *formula* : chemical formula Chemical formula. Any format accepted by :func:`.formulas.formula` can be used, including formula string. *mass* : float | g Sample mass. *name* : string Name of the sample (defaults to formula). """ def __init__(self, formula, mass, name=None): self.formula = build_formula(formula) self.mass = mass # cell F19 self.name = name if name else str(self.formula) # cell F20 self.activity = {} # The following are set in calculation_activation self.environment = None # type: "ActivationEnvironment" self.exposure = 0. self.rest_times = () def calculate_activation(self, environment, exposure=1, rest_times=(0, 1, 24, 360), abundance=table_abundance): """ Calculate sample activation (uCi) after exposure to a neutron flux. *environment* is the exposure environment. *exposure* is the exposure time in hours (default is 1 h). *rest_times* are deactivation times in hours (default [0, 1, 24, 360]). *abundance* is a function that returns the relative abundance of an isotope. By default it uses :func:`table_abundance` with natural abundance defined in :mod:`periodictable.mass`, but there is the alternative :func:`IAEA1987_isotopic_abundance` in the activation data table. """ self.activity = {} self.environment = environment self.exposure = exposure self.rest_times = rest_times for el, frac in self.formula.mass_fraction.items(): if core.isisotope(el): A = activity(el, self.mass*frac, environment, exposure, rest_times) self._accumulate(A) else: for iso in el.isotopes: iso_mass = self.mass*frac*abundance(el[iso])*0.01 if iso_mass: A = activity(el[iso], iso_mass, environment, exposure, rest_times) self._accumulate(A) def decay_time(self, target, tol=1e-10): """ After determining the activation, compute the number of hours required to achieve a total activation level after decay. """ if not self.rest_times or not self.activity: return 0 # Find the smallest rest time (probably 0 hr) k, t_k = min(enumerate(self.rest_times), key=lambda x: x[1]) # Find the activity at that time, and the decay rate data = [ (Ia[k], LN2/a.Thalf_hrs) for a, Ia in self.activity.items() # TODO: not sure why Ia is zero, but it messes up the initial value guess if it is there if Ia[k] > 0.0 ] # Need an initial guess near the answer otherwise find_root gets confused. # Small but significant activation with an extremely long half-life will # dominate at long times, but at short times they will not affect the # derivative. Choosing a time that satisfies the longest half-life seems # to work well enough. guess = max(-log(target/Ia)/La + t_k for Ia, La in data) # With times far from zero the time resolution in the exponential is # poor. Adjust the start time to the initial guess, rescaling intensities # to the activity at that time. adj = [(Ia*exp(-La*(guess-t_k)), La) for Ia, La in data] #print(adj) # Build f(t) = total activity at time T minus target activity and its # derivative df/dt. f(t) will be zero when activity is at target f = lambda t: sum(Ia*exp(-La*t) for Ia, La in adj) - target df = lambda t: sum(-La*Ia*exp(-La*t) for Ia, La in adj) #print("data", data, []) t, ft = find_root(0, f, df, tol=tol) percent_error = 100*abs(ft)/target if percent_error > 0.1: #return 1e100*365*24 # Return 1e100 rather than raising an error msg = ( "Failed to compute decay time correctly (%.1g error). Please" " report material, mass, flux and exposure.") % percent_error raise RuntimeError(msg) # Return time at least zero hours after removal from the beam. Correct # for time adjustment we used to stablize the fit. return max(t+guess, 0.0) def _accumulate(self, activity): for el, activity_el in activity.items(): el_total = self.activity.get(el, [0]*len(self.rest_times)) self.activity[el] = [T+v for T, v in zip(el_total, activity_el)] def show_table(self, cutoff=0.0001, format="%.4g"): """ Tabulate the daughter products. *cutoff=1* : float | uCi The minimum activation value to show. *format="%.1f"* : string The number format to use for the activation. """ # TODO: need format="auto" which picks an appropriate precision based on # cutoff and/or activation level. # Track individual rows with more than 1 uCi of activation, and total activation # Replace any activation below the cutoff with '---' rows = [] total = [0]*len(self.rest_times) for el, activity_el in sorted_activity(self.activity.items()): total = [t+a for t, a in zip(total, activity_el)] if all(a < cutoff for a in activity_el): continue activity_str = [format%a if a >= cutoff else "---" for a in activity_el] rows.append([el.isotope, el.daughter, el.reaction, el.Thalf_str]+activity_str) footer = ["", "", "", "total"] + [format%t if t >= cutoff else "---" for t in total] # If no significant total activation then don't print the table if all(t < cutoff for t in total): print("No significant activation") return # Print the table header, with an overbar covering the various rest times # Print a dashed separator above and below each column header = ["isotope", "product", "reaction", "half-life"] \ + ["%g hrs"%vi for vi in self.rest_times] separator = ["-"*8, "-"*9, "-"*8, "-"*10] + ["-"*12]*len(self.rest_times) cformat = "%-8s %-9s %8s %10s " + " ".join(["%12s"]*len(self.rest_times)) width = sum(len(c)+1 for c in separator[4:]) - 1 if width < 16: width = 16 overbar = "-"*(width//2-8) + " activity (uCi) " + "-"*((width+1)//2-8) offset = sum(len(c)+1 for c in separator[:4]) - 1 print(" "*(offset+1)+overbar) print(cformat%tuple(header)) print(cformat%tuple(separator)) # Print the significant table rows, or indicate that there were no # significant rows if the total is significant but none of the # individual isotopes if rows: for r in rows: print(cformat%tuple(r)) else: print("No significant isotope activation") print(cformat%tuple(separator)) # If there is more than one row, or if there is enough marginally # significant activation that the total is greater then the one row # print the total in the footer if len(rows) != 1 or any(c != t for c, t in zip(rows[0][4:], footer[4:])): print(cformat%tuple(footer)) print(cformat%tuple(separator)) def find_root(x, f, df, max=20, tol=1e-10): r""" Find zero of a function. Returns when $|f(x)| < tol$ or when max iterations have been reached, so check that $|f(x)|$ is small enough for your purposes. Returns x, f(x). """ fx = f(x) for _ in range(max): #print(f"step {_}: {x=} {fx=} df/dx={df(x)} dx={fx/df(x)}") if abs(fx) < tol: break x -= fx / df(x) fx = f(x) return x, fx def sorted_activity(activity_pair): """Interator over activity pairs sorted by isotope then daughter product.""" return sorted(activity_pair, key=lambda x: (x[0].isotope, x[0].daughter)) class ActivationEnvironment: """ Neutron activation environment. The activation environment provides details of the neutron flux at the sample position. *fluence* : float | n/cm^2/s Thermal neutron fluence on sample. For COLD neutrons enter equivalent thermal neutron fluence. **Warning**: For very high fluences, e.g., >E16 to E17 n/cm2, the equations give erroneous results because of the precision limitations. If there is doubt simply do the calculation at a lower flux and proportion the result. This will not work for the cascade reactions, i.e., two neutron additions. *Cd_ratio* : float Neutron cadmium ratio. Use 0 to suppress epithermal contribution. This is to account for those nuclides that have a significant contribution to the activation due to epithermal neutrons. This is tabulated in the cross-section database as the 'resonance cross- section'. Values can range from 4 to more than 100. ....... Use 0 for your initial calculation .......... If you do a specific measurement for the nuclide and spectrum of interest you simply apply a correction to the thermal based calculation, i.e., reduce the fluence by the appropriate factor. This computation can be based on a Cd ratio of a material that has no significant resonance cross-section, or that has been corrected so that it reflects just the thermal to epithermal fluence ratio. The computation simply adds a portion of this resonance cross-section to the thermal cross-section based on the presumption that the Cd ratio reflects the fluence ratio. Copper is a good material on which to equate the Cd ratio to the thermal-epithermal fluence ratio. For other materials, correct for the thermal:epithermal cross-section ratio. At the NBSR this ranges from 12 in mid-core, to 200 at RT-4, to >2000 at a filtered cold neutron guide position. *fast_ratio* : float Thermal/fast ratio needed for fast reactions. Use 0 to suppress fast contribution. This is very reaction dependent. You in essence need to know the answer before you do the calculation! That is, this ratio depends upon the shape of the cross-section curve as well as the spectrum shape above the energy threshold of the reaction. But at least you can do 'what-if' calculations with worse case assumptions (in the absence of specific ratios). The fast cross-sections in this database are weighted for a fast maxwellian spectrum so the fast/thermal ratio should be just a fluence correction (i.e., a fluence ratio), not an energy correction. Fast neutron cross-sections from IAEA273 are manually spectrum weighted. Those from IAEA156 are fission spectrum averaged as tabulated. Use 50 (for NBSR calculations) as a starting point if you simply exploring for possible products. """ def __init__(self, fluence=1e5, Cd_ratio=0., fast_ratio=0., location=""): self.fluence = fluence # cell F13 self.Cd_ratio = Cd_ratio # cell F15 self.fast_ratio = fast_ratio # cell F17 self.location = location # cell F21 # Cell Q1 @property def epithermal_reduction_factor(self): """ Used as a multiplier times the resonance cross section to add to the thermal cross section for all thermal induced reactions. """ return 1./self.Cd_ratio if self.Cd_ratio >= 1 else 0 COLUMN_NAMES = [ "_symbol", # 0 AF "_index", # 1 AG "Z", # 2 AH "symbol", # 3 AI "A", # 4 AJ "isotope", # 5 AK "abundance", # 6 AL "daughter", # 7 AM "_Thalf", # 8 AN "_Thalf_unit", # 9 AO "isomer", # 10 AP "percentIT", # 11 AQ "reaction", # 12 AR "fast", # 13 AS "thermalXS", # 14 AT "gT", # 15 AU "resonance", # 16 AV "Thalf_hrs", # 17 AW "Thalf_str", # 18 AX "Thalf_parent", # 19 AY "thermalXS_parent", # 20 AZ "resonance_parent", # 21 BA "comments", # 22 BB ] INT_COLUMNS = [1, 2, 4] BOOL_COLUMNS = [13] FLOAT_COLUMNS = [6, 11, 14, 15, 16, 17, 19, 20, 21] UNITS_TO_HOURS = {'y': 8760, 'd': 24, 'h': 1, 'm': 1/60, 's': 1/3600} def activity(isotope, mass, env, exposure, rest_times): """ Compute isotope specific daughter products after the given exposure time and rest period. Activations are listed in *isotope.neutron_activation*. Most of the activations (n,g n,p n,a n,2n) define a single step process, where a neutron is absorbed yielding the daughter and some prompt radiation. The daughter itself will decay during exposure, yielding a balance between production and emission. Any exposure beyound about eight halflives will not increase activity for that product. Activity for daughter products may undergo further neutron capture, reducing the activity of the daughter product but introducing a grand daughter with its own activity. The data tables for activation are only precise to about three significant figures. Any changes to the calculations below this threshold, e.g., due to slightly different mass or abundance, are therefore of little concern. Differences in formulas compare to the NCNR activation spreadsheet: * Column M: Use ln(2) rather than 0.693 * Column N: Use ln(2) rather than 0.693 * Column O: Rewrite to use expm1(x) = exp(x) - 1:: = L (1 - exp(-M Y43)/(1-(M/N) ) + exp(-N Y43)/((N/M)-1) ) = L (1 - N exp(-M Y43)/(N-M) + M exp(-N Y43)/(N-M) ) = L/(N-M) ( (N-M) - N exp(-M Y43) + M exp(-N Y43) ) = L/(N-M) ( N(1 - exp(-M Y43)) + M (exp(-N Y43) - 1) ) = L/(N-M) ( -N expm1(-M Y43) + M expm1(-N) ) = L/(N-M) (M expm1(-N Y43) - N expm1(-M Y43)) * Column X: Rewrite to use expm1(x) = exp(x) - 1:: = W ((abs(U)<1e-10 and abs(V)<1e-10) ? (V-U + (V-U)(V+U)/2) : (exp(-U)-exp(-V))) = W (exp(-U) - exp(-V)) = W exp(-V) (exp(V-U) - 1) = W exp(-U) (1 - exp(U-V)) = W exp(-V) expm1(V-U) = -W exp(-U) expm1(U-V) = (U > V) ? (W exp(-V) expm1(V-U)) : (-W exp(-U) expm1(U-V)) Differences in the data tables: * AW1462 (W-186 => W-188 2n) t1/2 in hrs is not converting days to hours * AK1495 (Au-198 => Au-199 2n) target should be Au-197 * AN1428 (Tm-169 => Tm-171 2n) t1/2 updated to Tm-171 rather than Tm-172 * AN1420 (Er-162 => Ho-163 b) t1/2 updated to 4570 y from 10 y * AT1508 (Pb-208 => Pb-209 act) Thermal (b) x 1000 to convert from mbarns to barns """ # TODO: is the table missing 1-H => 3-H ? # 0nly activations which produce radioactive daughter products are # included. Because 1-H => 2-H (act) is not in the table, is this why # there is no entry for 1-H => 3-H (2n). result = {} if not hasattr(isotope, 'neutron_activation'): return result for ai in isotope.neutron_activation: # Ignore fast neutron interactions if not using fast ratio if ai.fast and env.fast_ratio == 0: continue # Column D: elemental % mass content of sample # mass fraction and abundance already included in mass calculation, so not needed # Column E: target nuclide and comment # str(isotope), ai.comment # Column F: Nuclide Produced # ai.daughter # Column G: Half-life # ai.Thalf_str # Column H: initial effective cross-section [barn] # env.epithermal_reduction_factor:$Q$1 = 1/env.Cd_ratio:$F$15 initialXS = ai.thermalXS + env.epithermal_reduction_factor*ai.resonance # Column I: reaction # ai.reaction # Column J: fast? # ai.fast # Column K: effective reaction flux [n/cm^2/s] # env.fluence:$F$13 env.fast_ratio:$F$17 flux = env.fluence/env.fast_ratio if ai.fast else env.fluence # Column L: root part of activation calculation [uCi] # mass:$F$19 # Decay correction portion done in column M # The given mass is sample mass * sample fraction * isotope abundance # The constant converts from Bq to uCi via avogadro's number with # N_A[atoms] / 3.7e10[Bq/Ci] * 1e6 [uCi/Ci] ~ 1.627605611e19 Bq_to_uCi = 1.6276e19 #Bq_to_uCi = constants.avogadro_number / 3.7e4 # 1/(cm^2 s) (Bq s) (barn cm^2/barn) g / (g/mol) ((1/mol) / (Bq/uCi)) # TODO: use isotope.mass rather than isotope.isotope # Using isotope number rather than isotope mass to estimate the number # atoms introduces some error, overestimating Li activation by 0.25% # and underestimating Fe activation by 0.11% root = flux * (initialXS * 1e-24) * mass / isotope.isotope * Bq_to_uCi # Column M: 0.69/t1/2 [1/h] lambda of produced nuclide lam = LN2/ai.Thalf_hrs #print(ai.thermalXS, ai.resonance, env.epithermal_reduction_factor) #print(isotope, "D", mass, "F", ai.daughter, "G", ai.Thalf_str, # "H", initialXS, "I", ai.reaction, "J", ai.fast, "K", flux, # "L", root, "M", lam) # Column Y: activity at the end of irradiation [uCi] if ai.reaction == 'b': # Column N: 0.69/t1/2 [1/h] lambda of parent nuclide parent_lam = LN2 / ai.Thalf_parent # Column O: Activation if "b" mode production # 2022-05-18 PAK: addressed the following # Note: problems resulting from precision limitation not addressed # in "b" mode production # by rewriting equation to use expm1: # activity = root*(1 - exp(-lam*exposure)/(1 - (lam/parent_lam)) # + exp(-parent_lam*exposure)/((parent_lam/lam)-1)) # Let x1=-lam*exposure x2=-parent_lam*exposure a=x1/x2=lam/parent_lam # activity = root * (1 - exp(x1)/(1-a) + exp(x2)/(1/a - 1)) # = root * (1 - x2*exp(x1)/(x2-x1)) + x1*exp(x2)/(x2-x1)) # = root * ((x2-x1) - x2*exp(x1) + x1*exp(x2)) / (x2-x1) # = root * (x2*(1-exp(x1)) + x1*(exp(x2)-1)) / (x2-x1) # = root * (x1*expm1(x2) - x2*expm1(x1)) / (x2-x1) # = root * (lam*expm1(x2) - parent_lam*expm1(x1)) / (parent_lam - lam) # Checked for each b-mode production that small halflife results are # unchanged to four digits and Eu[151] => Gd[152] no longer fails. activity = root/(parent_lam - lam) * ( lam*expm1(-parent_lam*exposure) - parent_lam*expm1(-lam*exposure)) #print("N", parent_lam, "O", activity) elif ai.reaction == '2n': # Column N: 0.69/t1/2 [1/h] lambda of parent nuclide parent_lam = LN2 / ai.Thalf_parent # Column P: effective cross-section 2n product and n, g burnup [barn] # Note: This cross-section always uses the total thermal flux effectiveXS = ai.thermalXS_parent + env.epithermal_reduction_factor*ai.resonance_parent # Column Q: 2n mode effective lambda of stable target [1/h] lam_2n = (flux*3600)*(initialXS*1e-24) # Column R: radioactive parent [1/h] parent_activity = (env.fluence*3600)*(effectiveXS*1e-24) + parent_lam # Column S: resulting product [1/h] product_2n = lam if ai.reaction == '2n' else 0 # Column T: activity if 2n mode activity = root*lam*(parent_activity-parent_lam)*( (exp(-lam_2n*exposure) / ((parent_activity-lam_2n)*(product_2n-lam_2n))) + (exp(-parent_activity*exposure) / ((lam_2n-parent_activity)*(product_2n-parent_activity))) + (exp(-product_2n*exposure) / ((lam_2n-product_2n)*(parent_activity-product_2n))) ) #print("N", parent_lam, "P", effectiveXS, "Q", lam_2n, # "R", parent_activity, "S", product_2n, "T", activity) else: # Provide the fix for the limitied precision (15 digits) in the # floating point calculation. For neutron fluence rates above # 1e16 the precision in certain cells needs to be improved to # avoid erroneous results. Also, burnup for single capture # reactions (excluding 'b') is included here. # See README file for details. # Column P: effective cross-section 2n product and n, g burnup [barn] # Note: This cross-section always uses the total thermal flux effectiveXS = ai.thermalXS_parent + env.epithermal_reduction_factor*ai.resonance_parent # Column U: nv1s1t [neutrons] U = (flux*3600)*(initialXS*1e-24)*exposure # Column V: nv2s2t+L2*t [neutrons] V = ((env.fluence*3600)*(effectiveXS*1e-24)+lam)*exposure # Column W: L/(L-nvs1+nvs2) [unitless] W = lam/(lam-(flux*3600)*(initialXS*1e-24)+(env.fluence*3600)*(effectiveXS*1e-24)) # Column X: W*(exp(-U)-exp(V)) if U,V > 1e-10 else W*(V-U+(V-U)*(V+U)/2) # [PAK 2024-02-28] Rewrite the exponential difference using expm1() #Xp = exp(-V)*expm1(V-U) if U>V else -exp(-U)*expm1(U-V) X = W*exp(-V)*expm1(V-U) if U>V else -W*exp(-U)*expm1(U-V) # Column Y: O if "b" else T if "2n" else L*X activity = root*X #print(f"{ai.isotope}=>{ai.daughter} {U=} {V=} {W=} {Xp=} {X=} {activity=} {lam=} {flux*initialXS*3600*1e-24+env.fluence*effectiveXS*3600*1e-24}") #print(f"{ai.isotope}=>{ai.daughter} {U=} {V=} {W=} {X=} {activity=}") if activity < 0: msg = "activity %g less than zero for %g"%(activity, isotope) raise RuntimeError(msg) #print(ai.thermalXS_parent, ai.resonance_parent, exposure) #print("P", effectiveXS, "U", U, "V", V, "W", W, "X", # precision_correction, "Y", activity) # columns: F32 H K L U V W X #data = env.fluence, initialXS, flux, root, U, V, W, precision_correction #print(" ".join("%.5e"%v for v in data)) # TODO: chained activity (e.g., ) result[ai] = [activity*exp(-lam*Ti) for Ti in rest_times] #print([(Ti, Ai) for Ti, Ai in zip(rest_times, result[ai])]) return result class ActivationResult: r""" *isotope* : Activation target for this result ({symbol}-{A}) *abundance* : float | % IAEA 1987 isotopic abundance of activation target *symbol* : Element symbol for isotope *A* : int Number of protons plus neutrons in isotope *Z* : int Number of protons in isotope *reaction* : Activation type - "act" for (n,gamma) - "n,p" for (n,proton) - "n,a" for (n,alpha) - "2n" for activation of a daughter (e.g., 95Mo + n => 95Nb + n => 96Nb) - "n,2n" for neutron catalyzed release of a neutron - "b" for beta decay of a daughter (e.g., 98Mo + n => 99Mo => Tc-99) *comments* : Notes relating to simplifications or assumptions in the database. For most situations these do not affect your results. *daughter* : Daughter product from activation ({symbol}-{A}{isomer}) *isomer* : Daughter product isotope annotation m, m1, m2: indicate metastable states. Decay may be to the ground state or to another nuclide. \+: indicates radioactive daughter production already included in daughter listing several parent t1/2's required to acheive calculated daughter activity. All activity assigned at end of irradiation. In most cases the added activity to the daughter is small. \*: indicates radioactive daughter production NOT calculated, approximately secular equilibrium s: indicates radioactive daughter of this nuclide in secular equilibrium after several daughter t1/2's t: indicates transient equilibrium via beta decay. Accumulation of that nuclide during irradiation is separately calculated. *Thalf_hrs* : float | hours Half-life of daughter in hours *Thalf_str* : Half-life of daughter as string, such as "29.0 y" *Thalf_parent* : float | hours Half-life of parent in chained "2n" or "b" reaction *fast* : bool Indicates whether the reaction is fast or thermal. If fast then the fluence is reduced by the specified fast/thermal ratio. When *fast_ratio* is zero in the environment this activation will not appear. *thermalXS*, *resonance*, *thermalXS_parent*, *resonance_parent* : float | barns Activation database values for computing the reaction cross section. *gT*, *percentIT* : Unused. Database notes: Most of the cross-section data is from IAEA 273. This is an excellent compilation. I highly recommend its purchase. The fast neutron cross-section data entered in the spreadsheet is weighted by an U-235 maxwellian distributed fission spectrum. Allthermal reactions producing nuclides with a half-life in excess of 1 second, and some less than 1 second, are included. Also included are radioactive daughters with halflives substantially longer than the parent produced by the neutron induced reaction, i.e., those daughters not in secular equilibrium. See the database in for more detailed notes relating to the database entries. [Note: 150 fission spectrum averaged fast neutron reactions added from IAEA156 on 2/1/94. All those that are tabulated as measured have been entered. Those estimated by calculation have not been entered. In practice this means that the convenient, observable products are in this database. Fast reactions included are n,p; n,alpha; n,2n; and n,n'. Reaction = b : This is the beta produced daughter of an activated parent. This is calculated only for the cases where the daughter is long lived relative to the parent. In the reverse case the daughter activity is reasonably self evident and the parent nuclide is tagged to indicate a radioactive daughter. The calculated activity of the beta produced daughter is through the end of irradiation. Contributions from the added decay of the parent after the end of irradiation are left for the user to determine, but are usually negligible for irradiations that are long relative to the parent halflife. During irradiation parent activity (A1) is: A1 = K [1 - exp(-L1*t)] where K is the parent saturation activity and L1 is the decay constant. For the beta produced daughter the activity (A2) is: A2 = K [1 - exp(-L1*t) * L2/(L2-L1) + exp(-L2*t) * L1/(L2-L1)] """ def __init__(self, **kw): self.__dict__ = kw def __repr__(self): return f"ActivationResult({self.isotope},{self.reaction},{self.daughter})" def __str__(self): return f"{self.isotope}={self.reaction}=>{self.daughter}" def init(table, reload=False): """ Add neutron activation levels to each isotope. """ from math import floor # TODO: importlib.reload does not work for iso.neutron_activation attribute if 'neutron_activation' not in table.properties: table.properties.append('neutron_activation') elif not reload: return else: # Reloading activation table so clear the existing data. for el in table: for iso in el.isotopes: if hasattr(el[iso], 'neutron_activation'): del el[iso].neutron_activation # We are keeping the table as a simple export of the activation data # from the ncnr health physics excel spreadsheet so that it is easier # to validate that the table contains the same data. Unfortunately some # of the cells involved formulas, which need to be reproduced when loading # in order to match full double precision values. activations = {} path = os.path.join(core.get_data_path('.'), 'activation.dat') with open(path, 'r') as fd: for row in fd: #print(row, end='') columns = row.split('\t') if columns[0].strip() in ('', 'xx'): continue columns = [c[1:-1] if c.startswith('"') else c for c in columns] #print columns for c in INT_COLUMNS: columns[c] = int(columns[c]) for c in BOOL_COLUMNS: columns[c] = (columns[c] == 'y') for c in FLOAT_COLUMNS: s = columns[c].strip() columns[c] = float(s) if s else 0. # clean up comment column columns[-1] = columns[-1].replace('"', '').strip() kw = dict(zip(COLUMN_NAMES, columns)) iso = (kw['Z'], kw['A']) act = activations.setdefault(iso, []) # TODO: use NuBase2020 for halflife # NuBase2020 uses different isomer labels. # Strip the (+, *, s, t) tags # Convert m1/m2 to m/n for # In-114 In-116 Sb-124 Eu-152 Hf-179 Ir-192 # Convert Eu-152m2 to Eu-152r # Convert m to n for # Sc-46m Ge-73m Kr-83m Pd-107m Pd-109m Ag-110m Ta-182m Ir-194m Pb-204m # Convert m to p for # Sb-122m Lu-177m # Recreate Thalf_hrs column using double precision. # Note: spreadsheet is not converting half-life to hours in cell AW1462 (186-W => 188-W) kw['Thalf_hrs'] = float(kw['_Thalf']) * UNITS_TO_HOURS[kw['_Thalf_unit']] #print(f"T1/2 {kw['Thalf_hrs']} +/- {kw['Thalf_hrs_unc']}") # Recreate Thalf_parent by fetching from the new Thalf_hrs # e.g., =IF(OR(AR1408="2n",AR1408="b"),IF(AR1407="b",AW1406,AW1407),"") # This requires that the parent is directly before the 'b' or 'nb' # with its activation list already entered into the isotope. # Note: 150-Nd has 'act' followed by two consecutive 'b' entries. if kw['reaction'] in ('b', '2n'): parent = act[-2] if act[-1].reaction == 'b' else act[-1] kw['Thalf_parent'] = parent.Thalf_hrs else: #assert kw['Thalf_parent'] == 0 kw['Thalf_parent'] = None # Half-lives use My, Gy, Ty, Py value, units = float(kw['_Thalf']), kw['_Thalf_unit'] if units == 'y': if value >= 1e12: value, units = value/1e12, 'Ty' elif value >= 1e9: value, units = value/1e9, 'Gy' elif value >= 200e3: # above 200,000 years use My value, units = value/1e6, 'My' elif value >= 20e3: # between 20,000 and 200,000 use ky value, units = value/1e3, 'ky' formatted = f"{value:g} {units}" #if formatted.replace(' ', '') != kw['Thalf_str'].replace(' ', ''): # print(f"{kw['_index']}: old {kw['Thalf_str']} != {formatted} new") kw['Thalf_str'] = formatted # Strip columns whose names start with underscore kw = dict((k, v) for k, v in kw.items() if not k.startswith('_')) # Create an Activation record and add it to the isotope act.append(ActivationResult(**kw)) # Check abundance values #if abs(iso.abundance - kw['abundance']) > 0.001*kw['abundance']: # percent = 100*abs(iso.abundance - kw['abundance'])/kw['abundance'] # print "Abundance of", iso, "is", iso.abundance, \ # "but activation.dat has", kw['abundance'], "(%.1f%%)"%percent # Plug the activation products into the table for (Z, A), daughters in activations.items(): table[Z][A].neutron_activation = tuple(daughters) def demo(): # pragma: nocover import sys import argparse import periodictable as pt parser = argparse.ArgumentParser(description='Process some data with mass, flux, exposure, and decay options.') parser.add_argument('-m', '--mass', type=float, default=1, help='Specify the mass value (default: 1)') parser.add_argument('-f', '--flux', type=float, default=1e8, help='Specify the flux value (default: 1e8)') parser.add_argument('-e', '--exposure', type=float, default=10, help='Specify the exposure value (default: 10)') parser.add_argument('-d', '--decay', type=float, default=5e-4, help='Specify the decay value (default: 5e-4)') parser.add_argument('--cd-ratio', type=float, default=70, help='Specify the Cd ratio value (default: 70)') parser.add_argument('--fast-ratio', type=float, default=50, help='Specify the fast ratio value (default: 50)') parser.add_argument('formula', nargs='?', type=str, default=None, help='Specify the formula as a positional argument') args = parser.parse_args() formula = args.formula if formula is None: # Make sure all elements compute #formula = "".join(str(el) for el in pt.elements)[1:] formula = build_formula([(1/el.mass, el) for el in pt.elements][1:]) # Use an enormous mass to force significant activation of rare isotopes mass, fluence = 1e15, 1e8 env = ActivationEnvironment(fluence=args.flux, Cd_ratio=args.cd_ratio, fast_ratio=args.fast_ratio, location="BT-2") sample = Sample(formula, mass=args.mass) abundance = IAEA1987_isotopic_abundance #abundance=NIST2001_isotopic_abundance, sample.calculate_activation( env, exposure=args.exposure, rest_times=(0, 1, 24, 360), abundance=abundance, ) decay_time = sample.decay_time(args.decay) print(f"{args.mass} g {formula} for {args.exposure} hours at {args.flux} n/cm^2/s") print(f"Time to decay to {args.decay} uCi is {decay_time} hours.") sample.calculate_activation( env, exposure=args.exposure, rest_times=(0, 1, 24, 360, decay_time), abundance=abundance, ) sample.show_table(cutoff=0.0) ## Print a table of flux vs. activity so we can debug the ## precision_correction value in the activity() function. ## Note that you also need to uncomment the print statement ## at the end of activity() that shows the column values. #import numpy as np #sample = Sample('Co', mass=10) #for fluence in np.logspace(3, 20, 20-3+1): # env = ActivationEnvironment(fluence=fluence) # sample.calculate_activation( # env, exposure=exposure, rest_times=[0], # abundance=IAEA1987_isotopic_abundance) if __name__ == "__main__": demo() # pragma: nocover