getting reactivity from cross sections

[shimwell]

we can get reactivity from cross sections in addition to Bosch-Hale

DD example where the DD-_3He_n.txt file has been downloaded from https://github.com/shimwell/fusion_cross_sections

from fusion_neutron_utils import reactivity
import numpy as np

import matplotlib.pyplot as plt

import numpy as np
from matplotlib import rc
import matplotlib.pyplot as plt
from scipy.constants import e, k as kB

# Reactant masses in atomic mass units (u).
masses = {'D': 2.014, 'T': 3.016, '3He': 3.016, 
          'O16':15.9949, 'O18': 17.999, 'p': 1.0072, '17F': 17.0020, }

# Energy grid, 1 – 1000 keV, evenly spaced in log-space.
# Egrid = np.logspace(0, 4, 100)

def read_xsec(filename,
              collider_particle_mass,
              target_particle_mass,
              energy_scaling_factor=1.e3,
              xs_scaling_factor=1.e-28):
    """Read in cross section from filename and interpolate to energy grid."""

    E, xs = np.genfromtxt(filename, comments='#', unpack=True)

    m1, m2 = collider_particle_mass, target_particle_mass
    E *= m1 / (m1 + m2)

    E = E*energy_scaling_factor
    xs = xs*xs_scaling_factor

    return E,xs

# from endf site? units are in barns and MeV
DD_E, DD_xs = read_xsec(filename='D_D_-_3He_n.txt',
                  collider_particle_mass=masses['D'],
                  target_particle_mass=masses['D'],
                  xs_scaling_factor=1e-50,
                  energy_scaling_factor=1e3)

egrid=DD_E#np.linspace(1,1000, 20)
b_and_H = []
for temp in egrid: 
    react = reactivity(
        ion_temperature=temp,
        temperature_units='keV',
        reactivity_units='m^3/s',
        reaction='D+D=p+T',
        equation='Bosch-Hale'
    )
    # print(react)
    b_and_H.append(react)

fig, ax = plt.subplots()
# ax.loglog(DD_E, DD_xs, lw=2, label='D,T')

# ax.loglog(egrid, b_and_H, label='b and h')

print(b_and_H/DD_xs)
plt.plot(b_and_H/DD_xs)
plt.ylim(-1,1)
ax.legend()
plt.show()

DT example where

from fusion_neutron_utils import reactivity
import numpy as np

import matplotlib.pyplot as plt

import numpy as np
from matplotlib import rc
import matplotlib.pyplot as plt
from scipy.constants import e, k as kB
# rc('font', **{'family': 'serif', 'serif': ['Computer Modern'], 'size': 14})
# rc('text', usetex=True)

# Reactant masses in atomic mass units (u).
masses = {'D': 2.014, 'T': 3.016, '3He': 3.016, 
          'O16':15.9949, 'O18': 17.999, 'p': 1.0072, '17F': 17.0020, }

# Energy grid, 1 – 1000 keV, evenly spaced in log-space.
# Egrid = np.logspace(0, 4, 100)

def read_xsec(filename,
              collider_particle_mass,
              target_particle_mass,
              energy_scaling_factor=1.e3,
              xs_scaling_factor=1.e-28):
    """Read in cross section from filename and interpolate to energy grid."""

    E, xs = np.genfromtxt(filename, comments='#', unpack=True)

    m1, m2 = collider_particle_mass, target_particle_mass
    E *= m1 / (m1 + m2)

    E = E*energy_scaling_factor
    xs = xs*xs_scaling_factor

    return E,xs

# # D + T -> α + n
# from endf site? units are in barns and eV
DT_E, DT_xs = read_xsec(filename='D_T_-_a_n.txt',
                  collider_particle_mass=masses['D'],
                  target_particle_mass=masses['T'],
                  xs_scaling_factor=1e-50,
                  energy_scaling_factor=1e3)

egrid=DT_E#np.linspace(1,1000, 20)
b_and_H = []
for temp in egrid: 
    react = reactivity(
        ion_temperature=temp,
        temperature_units='keV',
        reactivity_units='m^3/s',
        reaction='D+T=n+a',
        equation='Bosch-Hale'
    )
    print(react)
    b_and_H.append(react)

fig, ax = plt.subplots()
ax.loglog(DT_E, DT_xs, lw=2, label='D,T')

# plt.plot(DT_E, DT_xs)

ax.loglog(egrid, b_and_H, label='b and h')
ax.legend()
plt.show()

[shimwell]

Using data downloaded from deuteron sublibrary https://www.nndc.bnl.gov/endf-releases/?version=B-VIII.1

import endf
>>> h2 = endf.IncidentNeutron.from_endf('d-001_H_002.endf')
>>> import endf
>>> h2.reactions
{2: <Reaction: MT=2 (n,elastic)>, 50: <Reaction: MT=50>, 600: <Reaction: MT=600 (n,p0)>}
>>> h2[2]
<Reaction: MT=2 (n,elastic)>
>>> b=h2[2] 
>>> b.xs
{'0K': <Tabulated1D: 22 points, 1 regions>}
>>> b.xs['0K']
<Tabulated1D: 22 points, 1 regions>
>>> t=b.xs['0K']       
>>> t.x
array([   93000.,    96078.,    99157.,   105310.,   111470.,   120710.,
         139180.,   157650.,   200750.,   250000.,   500000.,   686270.,
        1000000.,  1236100.,  1585900.,  2271800.,  3692900.,  5632400.,
        8078400., 11059000., 15000000., 20000000.])
>>> t.y

[jon-proximafusion]

Not sure how JANIS does this but they appear to have finer grain values for ENDF data than the original ENDF data itself. https://www.oecd-nea.org/janisweb/book/deuterons/H2/MT4/renderer/838