"""Physical scale of salting: fission-product abundance pattern (ENDF/B-VIII.0 thermal yields), mass of the surface convection zone versus Teff (Amard et al. 2019 STAREVOL tracks), and the mass of each element needed to raise a photosphere by a given factor.""" import os, re, csv, zipfile, urllib.request import numpy as np from common import * NFY_URL = 'https://www.nndc.bnl.gov/endf-b8.0/zips/ENDF-B-VIII.0_nfy.zip' TRACK_URL = 'https://cdsarc.cds.unistra.fr/ftp/J/A+A/631/A77/tracks/M00%sZ13V00.dat' MASSES = ['0p60', '0p70', '0p80', '0p90', '1p00', '1p10', '1p20', '1p30', '1p40', '1p50'] # first beta-stable isobar reached from the neutron-rich side, for each mass chain CHAIN_END = {80: 'se', 81: 'br', 82: 'se', 83: 'kr', 84: 'kr', 85: 'rb', 86: 'kr', 87: 'rb', 88: 'sr', 89: 'y', 90: 'zr', 91: 'zr', 92: 'zr', 93: 'nb', 94: 'zr', 95: 'mo', 96: 'zr', 97: 'mo', 98: 'mo', 99: 'ru', 100: 'mo', 101: 'ru', 102: 'ru', 103: 'rh', 104: 'ru', 105: 'pd', 106: 'pd', 107: 'ag', 108: 'pd', 109: 'ag', 110: 'pd', 111: 'cd', 112: 'cd', 113: 'cd', 114: 'cd', 115: 'in', 116: 'cd', 117: 'sn', 118: 'sn', 119: 'sn', 120: 'sn', 121: 'sb', 122: 'sn', 123: 'sb', 124: 'sn', 125: 'te', 126: 'te', 127: 'i', 128: 'te', 129: 'xe', 130: 'te', 131: 'xe', 132: 'xe', 133: 'cs', 134: 'xe', 135: 'ba', 136: 'xe', 137: 'ba', 138: 'ba', 139: 'la', 140: 'ce', 141: 'pr', 142: 'ce', 143: 'nd', 144: 'nd', 145: 'nd', 146: 'nd', 147: 'sm', 148: 'nd', 149: 'sm', 150: 'nd', 151: 'eu', 152: 'sm', 153: 'eu', 154: 'sm', 155: 'gd', 156: 'gd', 157: 'gd', 158: 'gd', 159: 'tb', 160: 'gd'} # long-lived chain members that hold the chain for the first ~1e5 yr ("young" waste): # Zr-93 (1.6 Myr), Tc-99 (0.21 Myr), Pd-107 (6.5 Myr), I-129 (16 Myr), Cs-135 (2.3 Myr) YOUNG_HOLD = {93: 'zr', 99: 'tc', 107: 'pd', 129: 'i', 135: 'cs'} def _endf_float(s): s = s.strip() if not s: return 0.0 m = re.match(r'^([+-]?\d*\.?\d+)([+-]\d+)$', s) if m: return float(m.group(1) + 'e' + m.group(2)) return float(s) def chain_yields(nuclide='U_235', energy=0.0253): p = RAW + 'fy/ENDF-B-VIII.0_nfy.zip' os.makedirs(RAW + 'fy', exist_ok=True) if not os.path.exists(p): with urllib.request.urlopen(NFY_URL, timeout=120) as r: open(p, 'wb').write(r.read()) z = zipfile.ZipFile(p) name = [n for n in z.namelist() if n.endswith('_%s.endf' % nuclide)][0] lines = [l for l in z.read(name).decode().splitlines() if l[70:75] == ' 8459'] vals = lambda l: [_endf_float(l[k * 11:(k + 1) * 11]) for k in range(6)] le = int(vals(lines[0])[2]); pos = 1; best = None for _ in range(le): h = vals(lines[pos]); pos += 1 E = h[0]; nn = int(h[4]); nfp = int(h[5]) nl = (nn + 5) // 6 data = [] for l in lines[pos:pos + nl]: data += vals(l) pos += nl data = np.array(data[:nn]).reshape(nfp, 4) if best is None or abs(E - energy) < abs(best[0] - energy): best = (E, data) E, data = best chain = {} for zafp, fps, yc, dyc in data: A = int(round(zafp)) % 1000 chain[A] = max(chain.get(A, 0.0), yc) # chain yield = largest cumulative yield in the chain return E, chain def elemental_yields(chain, young=False): el = {} for A, y in chain.items(): e = CHAIN_END.get(A) if young and A in YOUNG_HOLD: e = YOUNG_HOLD[A] if e is None: continue el[e] = el.get(e, 0.0) + y return el def fission_pattern(el_y, ref='nd', ref_dex=1.0): """[X/H] enhancement of each element when enough fission waste is added to raise `ref` by ref_dex.""" nsun = {e: 10 ** (A_SUN[e] - 12.0) for e in A_SUN if A_SUN[e] is not None} k = (10 ** ref_dex - 1.0) * nsun[ref] / el_y[ref] # fissions per H atom return {e: float(np.log10(1.0 + k * el_y.get(e, 0.0) / nsun[e])) for e in nsun}, k def mcz_table(): os.makedirs(RAW + 'tracks', exist_ok=True) rows = [] for m in MASSES: p = RAW + 'tracks/M00%sZ13V00.dat' % m if not os.path.exists(p) or os.path.getsize(p) < 1000: with urllib.request.urlopen(TRACK_URL % m, timeout=120) as r: open(p, 'wb').write(r.read()) f = open(p); hdr = f.readline().lstrip('#').split() tr = [] for line in f: q = line.split() if len(q) < 30: continue dd = dict(zip(hdr, q)) try: tr.append((float(dd['time']), 10 ** float(dd['logTeff']), float(dd['M']), float(dd['logg']), float(dd['M']) - float(dd['Mr_b_CE']), float(dd['tau_max']))) except ValueError: pass tr = np.array(tr) ms = tr[(tr[:, 0] > 2e8) & (tr[:, 3] > 3.9)] # main sequence for age in (1e9, 2e9, 4.57e9): if ms[:, 0].max() < age * 0.9: continue i = np.argmin(np.abs(ms[:, 0] - age)) rows.append((float(m.replace('p', '.')), ms[i, 0] / 1e9, ms[i, 1], ms[i, 3], ms[i, 4], ms[i, 5])) return rows def mass_for_dex(element, dex, mcz_msun): """kg of `element` needed to raise its abundance in a convective envelope of mcz_msun by `dex`.""" n_h = X_H * mcz_msun * MSUN / M_H n_x = n_h * 10 ** (A_SUN[element] - 12.0) return (10 ** dex - 1.0) * n_x * ATOMIC_MASS[element] * 1.66054e-27 def main(verbose=True): out = {} with open(DATA + 'fission_pattern.csv', 'w', newline='') as fh: w = csv.writer(fh) w.writerow(['fissile', 'waste_age', 'element', 'atoms_per_100_fissions', 'solar_logeps', 'dex_when_Nd_plus1.0', 'dex_when_Nd_plus1.5', 'dex_when_Nd_plus2.0']) for nuc in ('U_235', 'Pu_239'): E, chain = chain_yields(nuc) tot = sum(chain.values()) for young in (False, True): ey = elemental_yields(chain, young) pats = [fission_pattern(ey, 'nd', dx)[0] for dx in (1.0, 1.5, 2.0)] out[(nuc, young)] = (ey, pats) for e in ['rb', 'sr', 'y', 'zr', 'nb', 'mo', 'ru', 'rh', 'pd', 'ba', 'la', 'ce', 'pr', 'nd', 'sm', 'eu', 'gd']: w.writerow([nuc, 'young(<1e5yr)' if young else 'old(>1e8yr)', e, f'{100 * ey.get(e, 0):.3f}', A_SUN[e], f'{pats[0][e]:.3f}', f'{pats[1][e]:.3f}', f'{pats[2][e]:.3f}']) if young: w.writerow([nuc, 'young(<1e5yr)', 'tc', f'{100 * ey.get("tc", 0):.3f}', '', '', '', '']) if verbose: ey, pats = out[(nuc, False)] print(f'{nuc} thermal (E={E} eV): sum of chain yields = {tot:.3f} per fission') print(' atoms/100 fissions:', ' '.join(f'{e}:{100 * ey.get(e, 0):.1f}' for e in ['sr', 'y', 'zr', 'mo', 'ru', 'ba', 'la', 'ce', 'pr', 'nd', 'sm', 'eu'])) print(' pattern at Nd +1.0:', ' '.join(f'{e}:{pats[0][e]:+.2f}' for e in ['rb', 'sr', 'y', 'zr', 'mo', 'ru', 'ba', 'la', 'ce', 'pr', 'nd', 'sm', 'eu'])) print(' Tc-99 per 100 fissions (young waste):', 100 * out[(nuc, True)][0].get('tc', 0)) rows = mcz_table() with open(DATA + 'convective_envelope_mass.csv', 'w', newline='') as fh: w = csv.writer(fh) w.writerow(['mass_msun', 'age_gyr', 'teff_K', 'logg', 'mcz_msun', 'tau_conv_max_yr', 'kg_Nd_for_+1dex', 'kg_Pr_for_+1dex', 'kg_Nd_for_+2dex']) for r in rows: w.writerow([r[0], f'{r[1]:.2f}', f'{r[2]:.0f}', f'{r[3]:.2f}', f'{r[4]:.3e}', f'{r[5]:.3e}', f'{mass_for_dex("nd", 1.0, r[4]):.3e}', f'{mass_for_dex("pr", 1.0, r[4]):.3e}', f'{mass_for_dex("nd", 2.0, r[4]):.3e}']) if verbose: for r in rows: if abs(r[1] - 2.0) < 0.3 or (r[0] == 1.0): print(f'M={r[0]:.1f} age={r[1]:.1f} Gyr Teff={r[2]:.0f} Mcz={r[4]:.2e} Msun Nd(+1dex)={mass_for_dex("nd", 1.0, r[4]):.2e} kg Pr(+1dex)={mass_for_dex("pr", 1.0, r[4]):.2e} kg') for mcz, lab in [(0.025, 'Sun, helioseismic CZ 0.025 Msun'), (0.0156, 'STAREVOL 1 Msun 4.6 Gyr')]: print(lab, {e: f'{mass_for_dex(e, 1.0, mcz):.2e}' for e in ['pr', 'nd', 'eu', 'la', 'ce', 'ba', 'zr', 'sm']}) write_fission_excess() return out, rows def fission_excess_table(nuclide='U_235', young=False): """Observed excess statistic (Ba-anchored, as in s03) produced by adding fission waste to a star with a scaled-solar pattern, as a function of the resulting Nd enhancement.""" E, chain = chain_yields(nuclide); ey = elemental_yields(chain, young) S = S_FRAC; rows = [] for dx in (0.5, 0.75, 1.0, 1.25, 1.5, 2.0, 2.5, 3.0): p = fission_pattern(ey, 'nd', dx)[0] fs = max((10 ** p['ba'] - (1 - S['ba'])) / S['ba'], 0.05) row = dict(nuclide=nuclide, nd_enhancement=dx, ba=p['ba'], y=p['y'], zr=p['zr'], mo=p['mo'], ru=p['ru'], pr=p['pr']) for e in ('la', 'ce', 'nd', 'sm', 'eu'): pred = np.log10(S[e] * fs + (1 - S[e]) * 1.0) row[e] = p[e]; row['ex_' + e] = p[e] - pred rows.append(row) return rows def write_fission_excess(): rows = fission_excess_table('U_235') + fission_excess_table('Pu_239') cols = list(rows[0].keys()) with open(DATA + 'fission_salt_observed_excess.csv', 'w', newline='') as fh: w = csv.DictWriter(fh, fieldnames=cols); w.writeheader() for r in rows: w.writerow({k: (f'{v:.3f}' if isinstance(v, float) else v) for k, v in r.items()}) return rows if __name__ == '__main__': main()