bt_maxwellian_n Subroutine

subroutine bt_maxwellian_n(fn, T, eb, am, ab, n, rate)
    !+ Calculates Maxwellian reaction rate for a beam with atomic mass `ab`, energy `eb`, and energy level `n`
    !+firing into a target with atomic mass `am` and temperature `T` which has a cross section given by the function `fn`
    interface
        function fn(a, b)
            !+Cross section function
            real(8)              :: fn !sigma
            real(8), intent(in)  :: a !eb
            integer, intent(in)  :: b !n
        end function fn
    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
    real(Float64), intent(out) :: rate
        !+Reaction Rate [\(cm^3/s\)]

    logical :: dxc
    integer :: n_vr
    real(Float64) :: vr_max, dvr
    real(Float64), dimension(32) :: vr
    real(Float64), dimension(32) :: fr
    integer :: n_vz
    real(Float64) :: vz_max,dvz
    real(Float64), dimension(62) :: vz
    real(Float64), dimension(62) :: fz
    real(Float64) :: T_per_amu, eb_per_amu, ared, sig, sig_eff
    real(Float64) :: zb, u2_to_erel, u2, erel, v_therm, dE

    integer :: i, j

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

    n_vz = 62
    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

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

    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 !for electron interactions

    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(erel.ge.dE) then
                sig = fn(erel/ared,n)
            else
                sig = 0.d0
            endif
            fr(j) = 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

end subroutine bt_maxwellian_n