Use von Dreele approximation to Lynn & Seeger approximation for resonant elements

[pkienzle]

Resonance parameters from von Dreele (2024) DOI:10.1107/S1600576724005375

103-Rh has 100% natural abundance so natural Rh is a copy of 103-Rh.

113-Cd has 12.22%, which scales Re and Im for natural Cd.

b0 values for 103-Rh, Cd, 113-Cd, 239-Pu, 240-Pu were not provided in the paper. They can be determined by matching the b_c values at 2200 m/s recorded in nsf.py.

from periodictable.nsf import neutron_energy, _4PI_100
import numpy as np

resonant_isotopes = {
#    Iso: (b0      Re        Im        E0         G       A1 E1 A2 E2)
    "Rh": (??,   78.92,  6156.0,   1257.0,    156,      0, 0, 0, 0),
"103-Rh": (??,   78.92,  6156.0,   1257.0,    156,      0, 0, 0, 0),
# [PAK] 113Cd resonance from Volev (2013) DOI:10.1016/j.nimb.2013.01.030
# with b0 adjusted to match the thermal b_c.
# E0=178.70(10), Γ=113.5(1), gΓm=0.4800(15), λ0=(81.787/E0)**0.5, μ=(A+n)/(An)=1.0088
# Re=2*gΓm*λ0*μ/(8*π)*1e4=260.7(8), Im=2*gΓm*Γ*λ0*μ/(16*π)*1e4=14794(48)
    "Cd": (0.6657,  31.86,   1808.,   178.7,    113.5,    0, 0, 0, 0),
"113-Cd": (0.6949, 260.7,   14794.,   178.7,    113.5,    0, 0, 0 ,0),
    "Sm": (0.4099,  30.10,    949,     97.75,    65.14,   0, 0, 0, 0),
"149-Sm": (0.470,  216.4,    6830,     97.74,    65.15,   0, 0, 0, 0),
    "Eu": (0.7550,   8.55,    385.5,  321.01,    87.32,   7.14, 459.65, 1.25, -31.0),
"151-Eu": (0.731,   17.88,    807,    321.0,     87.35,   7.14, 459.64, 1.22, -30.5),
    "Gd": (0.6794,  72.72,   3866,     30.40,   105.60,   0, 0, 0, 0),
"155-Gd": (0.6820,  88.74,   4672.3,   28.069,  105.135,  0, 0, 0, 0),
"157-Gd": (0.6428, 379.9,   20178,     31.0194, 105.74,   0, 0, 0, 0),
    "Er": (0.8525,  12.08,    523.5,  460.30,    88.08,   0.6522, 584.42, 0, 0),
"167-Er": (0.5635,  52.42,   2281.9,  460.16,    87.93,   0.6481, 584.23, 0, 0),
    "Yb": (1.2273,   0.56,     36.2,  596.97,    66.165,  0, 0, 0, 0),
"168-Yb": (0.7046, 635.179, 21152.9,  596.9696,  66.1647, 0, 0, 0, 0),
"176-Lu": (0.8042,  25.55,    765,    141.36,    60.560,  0, 0, 0, 0),
"239-Pu": (??,   22.66,   1156.0,  296,      102.0,    0, 0, 0, 0),
"240-Pu": (??,   515.8,   8356.0, 1057,       32.4,    0, 0, 0, 0),
}

def init_resonance(table):
    for el in table:
        el.neutron.resonance = resonant_isotopes.get(str(el), None)
        for isonum in el.isotopes:
            iso = el[isonum]
            iso.neutron.resonance = resonant_isotopes.get(str(iso), None)

def scattering_by_wavelength(self, wavelength):
        r"""
        Return scattering length and total cross section for each wavelength.

        For select isotopes this returns the energy-dependent $\mathrm{Re}(b_c)$
        and $\mathrm{Im}(b_c)$ from low energy resonances. Total scattering is
        returned as $4\pi/100 |b_c|^2$ with no contribution for bound incoherent
        scattering. This is incorrect for odd-numbered isotopes.

        :Parameters:
            *wavelength* \: float(s) | |Ang|

        :Returns:
            *b_c* \: complex(s) | fm

            *sigma_s* \: float(s) | barn
        """
        coeff = self.resonance
        if coeff is None:
            ones = 1 if np.isscalar(wavelength) else np.ones_like(wavelength)
            return ones*self.b_c_complex, ones*self.total
        E = neutron_energy(wavelength)
        b0, Re, Im, E0, G, A1, E1, A2, E2 = coeff
        bp, bpp = b0, 0
        for Ak, Ek in ((1, E0), (A1, E1), (A2, E2)):
            if Ak == 0: break
            dE = (E - Ek)
            scale = Ak / (dE**2 + G**2/4)
            bp += Re * dE * scale
            bpp += Im * scale
        b_c = 10*(bp - 1j*bpp)
        # TODO: sigma_s should include an incoherent contribution?
        sigma_s = _4PI_100*abs(b_c)**2 # 1 barn = 1 fm^2 1e-2 barn/fm^2
        return b_c, sigma_s