OpenMC/openmc/data/reconstruct.pyx

from libc.stdlib cimport malloc, calloc, free
from libc.math cimport cos, sin, sqrt, atan, M_PI

cimport numpy as np
import numpy as np
from numpy.linalg import inv
cimport cython


cdef extern from "complex.h":
    double cabs(double complex)
    double complex conj(double complex)
    double creal(complex double)
    double cimag(complex double)
    double complex cexp(double complex)

# Physical constants are from CODATA 2014
cdef double NEUTRON_MASS_ENERGY = 939.5654133e6  # eV/c^2
cdef double HBAR_C = 197.3269788e5  # eV-b^0.5


@cython.cdivision(True)
def wave_number(double A, double E):
    r"""Neutron wave number in center-of-mass system.

    ENDF-102 defines the neutron wave number in the center-of-mass system in
    Equation D.10 as

    .. math::
        k = \frac{2m_n}{\hbar} \frac{A}{A + 1} \sqrt{|E|}

    Parameters
    ----------
    A : double
        Ratio of target mass to neutron mass
    E : double
        Energy in eV

    Returns
    -------
    double
        Neutron wave number in b^-0.5

    """
    return A/(A + 1)*sqrt(2*NEUTRON_MASS_ENERGY*abs(E))/HBAR_C

@cython.cdivision(True)
cdef double _wave_number(double A, double E):
    return A/(A + 1)*sqrt(2*NEUTRON_MASS_ENERGY*abs(E))/HBAR_C


@cython.cdivision(True)
cdef double phaseshift(int l, double rho):
    """Calculate hardsphere phase shift as given in ENDF-102, Equation D.13

    Parameters
    ----------
    l : int
        Angular momentum quantum number
    rho : float
        Product of the wave number and the channel radius

    Returns
    -------
    double
        Hardsphere phase shift

    """
    if l == 0:
        return rho
    elif l == 1:
        return rho - atan(rho)
    elif l == 2:
        return rho - atan(3*rho/(3 - rho**2))
    elif l == 3:
        return rho - atan((15*rho - rho**3)/(15 - 6*rho**2))
    elif l == 4:
        return rho - atan((105*rho - 10*rho**3)/(105 - 45*rho**2 + rho**4))


@cython.cdivision(True)
def penetration_shift(int l, double rho):
    r"""Calculate shift and penetration factors as given in ENDF-102, Equations D.11
    and D.12.

    Parameters
    ----------
    l : int
        Angular momentum quantum number
    rho : float
        Product of the wave number and the channel radius

    Returns
    -------
    double
        Penetration factor for given :math:`l`
    double
        Shift factor for given :math:`l`

    """
    cdef double den

    if l == 0:
        return rho, 0.
    elif l == 1:
        den = 1 + rho**2
        return rho**3/den, -1/den
    elif l == 2:
        den = 9 + 3*rho**2 + rho**4
        return rho**5/den, -(18 + 3*rho**2)/den
    elif l == 3:
        den = 225 + 45*rho**2 + 6*rho**4 + rho**6
        return rho**7/den, -(675 + 90*rho**2 + 6*rho**4)/den
    elif l == 4:
        den = 11025 + 1575*rho**2 + 135*rho**4 + 10*rho**6 + rho**8
        return rho**9/den, -(44100 + 4725*rho**2 + 270*rho**4 + 10*rho**6)/den


@cython.boundscheck(False)
@cython.wraparound(False)
@cython.cdivision(True)
def reconstruct_mlbw(mlbw, double E):
    """Evaluate cross section using MLBW data.

    Parameters
    ----------
    mlbw : openmc.data.MultiLevelBreitWigner
        Multi-level Breit-Wigner resonance parameters
    E : double
        Energy in eV at which to evaluate the cross section

    Returns
    -------
    elastic : double
        Elastic scattering cross section in barns
    capture : double
        Radiative capture cross section in barns
    fission : double
        Fission cross section in barns

    """
    cdef int i, nJ, ij, l, n_res, i_res
    cdef double elastic, capture, fission
    cdef double A, k, rho, rhohat, I
    cdef double P, S, phi, cos2phi, sin2phi
    cdef double Ex, Q, rhoc, rhochat, P_c, S_c
    cdef double jmin, jmax, j, Dl
    cdef double E_r, gt, gn, gg, gf, gx, P_r, S_r, P_rx
    cdef double gnE, gtE, Eprime, x, f
    cdef double *g
    cdef double (*s)[2]
    cdef double [:,:] params

    I = mlbw.target_spin
    A = mlbw.atomic_weight_ratio
    k = _wave_number(A, E)

    elastic = 0.
    capture = 0.
    fission = 0.

    for i, l in enumerate(mlbw._l_values):
        params = mlbw._parameter_matrix[l]

        rho = k*mlbw.channel_radius[l](E)
        rhohat = k*mlbw.scattering_radius[l](E)
        P, S = penetration_shift(l, rho)
        phi = phaseshift(l, rhohat)
        cos2phi = cos(2*phi)
        sin2phi = sin(2*phi)

        # Determine shift and penetration at modified energy
        if mlbw._competitive[i]:
            Ex = E + mlbw.q_value[l]*(A + 1)/A
            rhoc = mlbw.channel_radius[l](Ex)
            rhochat = mlbw.scattering_radius[l](Ex)
            P_c, S_c = penetration_shift(l, rhoc)
            if Ex < 0:
                P_c = 0

        # Determine range of total angular momentum values based on equation
        # 41 in LA-UR-12-27079
        jmin = abs(abs(I - l) - 0.5)
        jmax = I + l + 0.5
        nJ = int(jmax - jmin + 1)

        # Determine Dl factor using Equation 43 in LA-UR-12-27079
        Dl = 2*l + 1
        g = <double *> malloc(nJ*sizeof(double))
        for ij in range(nJ):
            j = jmin + ij
            g[ij] = (2*j + 1)/(4*I + 2)
            Dl -= g[ij]

        s = <double (*)[2]> calloc(2*nJ, sizeof(double))
        for i_res in range(params.shape[0]):
            # Copy resonance parameters
            E_r = params[i_res, 0]
            j = params[i_res, 2]
            ij = int(j - jmin)
            gt = params[i_res, 3]
            gn = params[i_res, 4]
            gg = params[i_res, 5]
            gf = params[i_res, 6]
            gx = params[i_res, 7]
            P_r = params[i_res, 8]
            S_r = params[i_res, 9]
            P_rx = params[i_res, 10]

            # Calculate neutron and total width at energy E
            gnE = P*gn/P_r  # ENDF-102, Equation D.7
            gtE = gnE + gg + gf
            if gx > 0:
                gtE += gx*P_c/P_rx

            Eprime = E_r + (S_r - S)/(2*P_r)*gn  # ENDF-102, Equation D.9
            x = 2*(E - Eprime)/gtE    # LA-UR-12-27079, Equation 26
            f = 2*gnE/(gtE*(1 + x*x)) # Common factor in Equation 40
            s[ij][0] += f             # First sum in Equation 40
            s[ij][1] += f*x           # Second sum in Equation 40
            capture += f*g[ij]*gg/gtE
            if gf > 0:
                fission += f*g[ij]*gf/gtE

        for ij in range(nJ):
            # Add all but last term of LA-UR-12-27079, Equation 40
            elastic += g[ij]*((1 - cos2phi - s[ij][0])**2 +
                              (sin2phi + s[ij][1])**2)

        # Add final term with Dl from Equation 40
        elastic += 2*Dl*(1 - cos2phi)

        # Free memory
        free(g)
        free(s)

    capture *= 2*M_PI/(k*k)
    fission *= 2*M_PI/(k*k)
    elastic *= M_PI/(k*k)

    return (elastic, capture, fission)


@cython.boundscheck(False)
@cython.wraparound(False)
@cython.cdivision(True)
def reconstruct_slbw(slbw, double E):
    """Evaluate cross section using SLBW data.

    Parameters
    ----------
    slbw : openmc.data.SingleLevelBreitWigner
        Single-level Breit-Wigner resonance parameters
    E : double
        Energy in eV at which to evaluate the cross section

    Returns
    -------
    elastic : double
        Elastic scattering cross section in barns
    capture : double
        Radiative capture cross section in barns
    fission : double
        Fission cross section in barns

    """
    cdef int i, l, i_res
    cdef double elastic, capture, fission
    cdef double A, k, rho, rhohat, I
    cdef double P, S, phi, cos2phi, sin2phi, sinphi2
    cdef double Ex, rhoc, rhochat, P_c, S_c
    cdef double E_r, J, gt, gn, gg, gf, gx, P_r, S_r, P_rx
    cdef double gnE, gtE, Eprime, f
    cdef double x, theta, psi, chi
    cdef double [:,:] params

    I = slbw.target_spin
    A = slbw.atomic_weight_ratio
    k = _wave_number(A, E)

    elastic = 0.
    capture = 0.
    fission = 0.

    for i, l in enumerate(slbw._l_values):
        params = slbw._parameter_matrix[l]

        rho = k*slbw.channel_radius[l](E)
        rhohat = k*slbw.scattering_radius[l](E)
        P, S = penetration_shift(l, rho)
        phi = phaseshift(l, rhohat)
        cos2phi = cos(2*phi)
        sin2phi = sin(2*phi)
        sinphi2 = sin(phi)**2

        # Add potential scattering -- first term in ENDF-102, Equation D.2
        elastic += 4*M_PI/(k*k)*(2*l + 1)*sinphi2

        # Determine shift and penetration at modified energy
        if slbw._competitive[i]:
            Ex = E + slbw.q_value[l]*(A + 1)/A
            rhoc = k*slbw.channel_radius[l](Ex)
            rhochat = k*slbw.scattering_radius[l](Ex)
            P_c, S_c = penetration_shift(l, rhoc)
            if Ex < 0:
                P_c = 0

        for i_res in range(params.shape[0]):
            # Copy resonance parameters
            E_r = params[i_res, 0]
            J = params[i_res, 2]
            gt = params[i_res, 3]
            gn = params[i_res, 4]
            gg = params[i_res, 5]
            gf = params[i_res, 6]
            gx = params[i_res, 7]
            P_r = params[i_res, 8]
            S_r = params[i_res, 9]
            P_rx = params[i_res, 10]

            # Calculate neutron and total width at energy E
            gnE = P*gn/P_r  # Equation D.7
            gtE = gnE + gg + gf
            if gx > 0:
                gtE += gx*P_c/P_rx

            Eprime = E_r + (S_r - S)/(2*P_r)*gn  # Equation D.9
            gJ = (2*J + 1)/(4*I + 2)  # Mentioned in section D.1.1.4

            # Calculate common factor for elastic, capture, and fission
            # cross sections
            f = M_PI/(k*k)*gJ*gnE/((E - Eprime)**2 + gtE**2/4)

            # Add contribution to elastic per Equation D.2
            elastic += f*(gnE*cos2phi - 2*(gg + gf)*sinphi2
                          + 2*(E - Eprime)*sin2phi)

            # Add contribution to capture per Equation D.3
            capture += f*gg

            # Add contribution to fission per Equation D.6
            if gf > 0:
                fission += f*gf

    return (elastic, capture, fission)


@cython.boundscheck(False)
@cython.wraparound(False)
@cython.cdivision(True)
def reconstruct_rm(rm, double E):
    """Evaluate cross section using Reich-Moore data.

    Parameters
    ----------
    rm : openmc.data.ReichMoore
        Reich-Moore resonance parameters
    E : double
        Energy in eV at which to evaluate the cross section

    Returns
    -------
    elastic : double
        Elastic scattering cross section in barns
    capture : double
        Radiative capture cross section in barns
    fission : double
        Fission cross section in barns

    """
    cdef int i, l, m, n, i_res
    cdef int i_s, num_s, i_J, num_J
    cdef double elastic, capture, fission, total
    cdef double A, k, rho, rhohat, I
    cdef double P, S, phi
    cdef double smin, smax, s, Jmin, Jmax, J, j
    cdef double E_r, gn, gg, gfa, gfb, P_r
    cdef double E_diff, abs_value, gJ
    cdef double Kr, Ki, x
    cdef double complex Ubar, U_, factor
    cdef bint hasfission
    cdef np.ndarray[double, ndim=2] one
    cdef np.ndarray[double complex, ndim=2] K, Imat, U
    cdef double [:,:] params

    # Get nuclear spin
    I = rm.target_spin

    elastic = 0.
    fission = 0.
    total = 0.
    A = rm.atomic_weight_ratio
    k = _wave_number(A, E)
    one = np.eye(3)
    K = np.zeros((3,3), dtype=complex)

    for i, l in enumerate(rm._l_values):
        # Check for l-dependent scattering radius
        rho = k*rm.channel_radius[l](E)
        rhohat = k*rm.scattering_radius[l](E)

        # Calculate shift and penetrability
        P, S = penetration_shift(l, rho)

        # Calculate phase shift
        phi = phaseshift(l, rhohat)

        # Calculate common factor on collision matrix terms (term outside curly
        # braces in ENDF-102, Eq. D.27)
        Ubar = cexp(-2j*phi)

        # The channel spin is the vector sum of the target spin, I, and the
        # neutron spin, 1/2, so can take on values of |I - 1/2| < s < I + 1/2
        smin = abs(I - 0.5)
        smax = I + 0.5
        num_s = int(smax - smin + 1)

        for i_s in range(num_s):
            s = i_s + smin

            # Total angular momentum is the vector sum of l and s and can assume
            # values between |l - s| < J < l + s
            Jmin = abs(l - s)
            Jmax = l + s
            num_J = int(Jmax - Jmin + 1)

            for i_J in range(num_J):
                J = i_J + Jmin

                # Initialize K matrix
                for m in range(3):
                    for n in range(3):
                        K[m,n] = 0.0

                hasfission = False
                if (l, J) in rm._parameter_matrix:
                    params = rm._parameter_matrix[l, J]

                    for i_res in range(params.shape[0]):
                        # Sometimes, the same (l, J) quantum numbers can occur
                        # for different values of the channel spin, s. In this
                        # case, the sign of the channel spin indicates which
                        # spin is to be used. If the spin is negative assume
                        # this resonance comes from the I - 1/2 channel and vice
                        # versa.
                        j = params[i_res, 2]
                        if l > 0:
                            if (j < 0 and s != smin) or (j > 0 and s != smax):
                                continue

                        # Copy resonance parameters
                        E_r = params[i_res, 0]
                        gn = params[i_res, 3]
                        gg = params[i_res, 4]
                        gfa = params[i_res, 5]
                        gfb = params[i_res, 6]
                        P_r = params[i_res, 7]

                        # Calculate neutron width at energy E
                        gn = sqrt(P*gn/P_r)

                        # Calculate j/2 * inverse of denominator of K matrix terms
                        factor = 0.5j/(E_r - E - 0.5j*gg)

                        # Upper triangular portion of K matrix -- see ENDF-102,
                        # Equation D.28
                        K[0,0] = K[0,0] + gn*gn*factor
                        if gfa != 0.0 or gfb != 0.0:
                            # Negate fission widths if necessary
                            gfa = (-1 if gfa < 0 else 1)*sqrt(abs(gfa))
                            gfb = (-1 if gfb < 0 else 1)*sqrt(abs(gfb))

                            K[0,1] = K[0,1] + gn*gfa*factor
                            K[0,2] = K[0,2] + gn*gfb*factor
                            K[1,1] = K[1,1] + gfa*gfa*factor
                            K[1,2] = K[1,2] + gfa*gfb*factor
                            K[2,2] = K[2,2] + gfb*gfb*factor
                            hasfission = True

                # Get collision matrix
                gJ = (2*J + 1)/(4*I + 2)
                if hasfission:
                    # Copy upper triangular portion of K to lower triangular
                    K[1,0] = K[0,1]
                    K[2,0] = K[0,2]
                    K[2,1] = K[1,2]

                    Imat = inv(one - K)
                    U = Ubar*(2*Imat - one)  # ENDF-102, Eq. D.27
                    elastic += gJ*cabs(1 - U[0,0])**2  # ENDF-102, Eq. D.24
                    total += 2*gJ*(1 - creal(U[0,0]))  # ENDF-102, Eq. D.23

                    # Calculate fission from ENDF-102, Eq. D.26
                    fission += 4*gJ*(cabs(Imat[1,0])**2 + cabs(Imat[2,0])**2)
                else:
                    U_ = Ubar*(2/(1 - K[0,0]) - 1)
                    if abs(creal(K[0,0])) < 3e-4 and abs(phi) < 3e-4:
                        # If K and phi are both very small, the calculated cross
                        # sections can lose precision because the real part of U
                        # ends up very close to unity. To get around this, we
                        # use Euler's formula to express Ubar by real and
                        # imaginary parts, expand cos(2phi) = 1 - 2phi^2 +
                        # O(phi^4), and then simplify
                        Kr = creal(K[0,0])
                        Ki = cimag(K[0,0])
                        x = 2*(-Kr + (Kr*Kr + Ki*Ki)*(1 - phi*phi) + phi*phi -
                               sin(2*phi)*Ki)/((1 - Kr)*(1 - Kr) + Ki*Ki)
                        total += 2*gJ*x
                        elastic += gJ*(x*x + cimag(U_)**2)
                    else:
                        total += 2*gJ*(1 - creal(U_))   # ENDF-102, Eq. D.23
                        elastic += gJ*cabs(1 - U_)**2   # ENDF-102, Eq. D.24

    # Calculate capture as difference of other cross sections as per ENDF-102,
    # Equation D.25
    capture = total - elastic - fission

    elastic *= M_PI/(k*k)
    capture *= M_PI/(k*k)
    fission *= M_PI/(k*k)

    return (elastic, capture, fission)