d_t_fusion Function

function d_t_fusion(eb) result(sigma)
    !+Calculates total cross section at a given deuterium energy, `eb`,
    !+for Deuterium-Tritium nuclear reactions in the range [0.5-550 keV]
    !+
    !+###Equation
    !+$$ D + T \rightarrow He^4(3.5 MeV) + n(14.1 MeV)$$
    !+
    !+###References
    !+* Equations 8-9
    !+* Table IV, VI in Ref. 7 [[atomic_tables(module)]]
    real(Float64), intent(in) :: eb
        !+ Deuterium energy [keV]
    real(Float64)             :: sigma
        !+ Cross Section [\(cm^2\)]

    real(Float64), dimension(5), parameter :: A1 = [ 6.927d4,  7.454d8, &
                                                    2.050d6, 5.2002d4,  &
                                                    0.d0 ]
    real(Float64), dimension(4), parameter :: B1 = [ 6.38d1,  -9.95d-1, &
                                                    6.981d-5, 1.728d-4  ]

    real(Float64), dimension(5), parameter :: A2 = [-1.4714d6, 0.d0,  &
                                                     0.d0, 0.d0, 0.d0 ]
    real(Float64), dimension(4), parameter :: B2 = [-8.4127d-3, 4.7983d-6, &
                                                    -1.0748d-9, 8.5184d-14 ]
    real(Float64), parameter :: Bg = 34.3827

    real(Float64), dimension(5) :: A
    real(Float64), dimension(4) :: B
    real(Float64) :: S, E

    E = min(max(eb,0.5),4700.0)

    if(E.le.530.0) then
        A = A1
        B = B1
    else
        A = A2
        B = B2
    endif

    S = (A(1) + E*(A(2) + E*(A(3) + E*(A(4) + E*A(5))))) / &
        (1    + E*(B(1) + E*(B(2) + E*(B(3) + E*B(4)))))

    sigma = (1.0d-27)*(S/(E*exp(Bg/sqrt(E))))

end function d_t_fusion