bt_maxwellian_n_m Subroutine

subroutine bt_maxwellian_n_m(fnm, T, eb, am, ab, n, m, rate, deexcit)
    !+ Calculates Maxwellian reaction rate for a `n`\(\rightarrow)`m` transition due to a beam with atomic mass `ab` and energy `eb`
    !+firing into a target with atomic mass `am` and temperature `T` which has a cross section given by the function `fnm`
    interface
        function fnm(a, b, c)
            !+Cross section function
            real(8)             :: fnm !sigma
            real(8), intent(in) :: a !eb
            integer, intent(in) :: b !n
            integer, intent(in) :: c !m
        end function fnm
    end interface
    real(Float64), intent(in)     :: T
        !+Target temperature [keV]
    real(Float64), intent(in)     :: eb
        !+Beam energy [keV]
    real(Float64), intent(in)     :: am
        !+Target atomic mass [amu]
    real(Float64), intent(in)     :: ab
        !+Beam atomic mass [amu]
    integer, intent(in)           :: n
        !+Initial atomic energy level/state
    integer, intent(in)           :: m
        !+Final atomic energy level/state
    real(Float64), intent(out)    :: rate
        !+Reaction Rate [\(cm^3/s\)]
    logical, intent(in), optional :: deexcit
        !+Calculate de-excitation reaction rate

    logical :: dxc
    integer, parameter :: n_vr = 30
    real(Float64) :: vr_max, dvr
    real(Float64), dimension(n_vr) :: vr
    real(Float64), dimension(n_vr) :: fr
    integer, parameter :: n_vz = 60
    real(Float64) :: vz_max,dvz
    real(Float64), dimension(n_vz) :: vz
    real(Float64), dimension(n_vz) :: fz
    real(Float64) :: T_per_amu, eb_per_amu, ared, sig, sig_eff
    real(Float64) :: zb, u2_to_erel, u2, erel, dE, factor, En, Em, v_therm

    integer :: i, j

    if(present(deexcit)) then
        dxc = deexcit
    else
        dxc = .False.
    endif

    vr_max = 4.d0
    dvr = vr_max/(n_vr - 1.d0)
    do i=1,n_vr
        vr(i) = (i-1)*dvr
    enddo

    vz_max = 4.d0
    dvz = 2.0*vz_max/(n_vz - 1.d0)
    do i=1,n_vz
        vz(i) = (i-1)*dvz - vz_max
    enddo

    En = (13.6d-3)*(1.0 - (1.d0/n)**2.0)
    Em = (13.6d-3)*(1.0 - (1.d0/m)**2.0)
    dE = Em - En

    T_per_amu = max(T, 1.d-6)/am
    eb_per_amu = eb/ab
    ared = am*ab/(am + ab)

    v_therm = 1.384d6 * sqrt(T_per_amu*1.d3)
    zb = sqrt(eb_per_amu/T_per_amu)
    u2_to_erel = ared*T_per_amu
    if(ared.lt.0.5) ared = 1.0

    fz = 0.d0
    fr = 0.d0

    do i=1,n_vz
        do j=1,n_vr
            u2 = (zb - vz(i))**2.0 + vr(j)**2.0
            erel = u2_to_erel*u2
            if(dxc) then
                factor = (erel + dE)/erel
                erel = erel + dE
            else
                factor = 1.0
            endif
            if(erel.ge.dE) then
                sig = fnm(erel/ared, n, m)
            else
                sig = 0.d0
            endif
            fr(j) = factor*sig*sqrt(u2)*exp(-(vz(i)**2.0 + vr(j)**2.0))*vr(j)
        enddo
        fz(i) = simpsons_rule(fr, dvr)
    enddo

    sig_eff = (2.0/sqrt(PI))*simpsons_rule(fz, dvz)
    rate = sig_eff*v_therm
    if(dxc) rate = rate*(real(n)/real(m))**2.0

end subroutine bt_maxwellian_n_m