p_excit_n Function

function p_excit_n(eb, n, m_max) result(sigma)
    !+Calculates an array of cross sections for a proton-Hydrogen impact excitation transitions from
    !+the `n` state to the \(m=1..{m_max}\) state at energy `eb`
    !+
    !+###Equation
    !+$$ H^+ + H(n) \rightarrow H^+ + H(m=1..m_{max}), m \gt n $$
    !+
    !+###References
    !+* Eq. 29.b and Table 4 in Ref. 2 for \(n = 1\) and \(m = 2\) [[atomic_tables(module)]]
    !+* Eq. 30 and Table 5 in Ref. 2 for \(n = 1\) and \(m = 3-6\) [[atomic_tables(module)]]
    !+* Eq. 31 and Table 5 in Ref. 2 for \(n = 1\) and \(m \gt 6\) [[atomic_tables(module)]]
    !+* Eq. 32 and Table 6 in Ref. 2 for \(n = 2\) and \(m \le 5\) [[atomic_tables(module)]]
    !+* Eq. 33 and Table 6 in Ref. 2 for \(n = 2\) and \(m = 6-10\) [[atomic_tables(module)]]
    !+* Eq. 34 and Table 6 in Ref. 2 for \(n = 2\) and \(m \gt 10\) [[atomic_tables(module)]]
    !+* Eq. 35 and Table 7 in Ref. 2 for \(n = 3\) and \(m \le 6\) [[atomic_tables(module)]]
    !+* Eq. 36 and Table 7 in Ref. 2 for \(n = 3\) and \(m = 7-10\) [[atomic_tables(module)]]
    !+* Eq. 37 and Table 7 in Ref. 2 for \(n = 3\) and \(m \gt 10\) [[atomic_tables(module)]]
    !+* Eq. 38-39 in Ref. 2 for \(n \gt 3\) and \(m \gt 4\) [[atomic_tables(module)]]
    !+
    real(Float64), intent(in)       :: eb
        !+ Relative collision energy [keV/amu]
    integer, intent(in)             :: n
        !+ Initial atomic energy level/state
    integer, intent(in)             :: m_max
        !+ Number of m states to calculate
    real(Float64), dimension(m_max) :: sigma
        !+ Array of cross sections where the m'th index refer to the transition
        !+ from `n` to m [\(cm^2\)]

    integer :: m
    real(Float64) :: nf, mf, Etil, s, D, A, G, L, F
    real(Float64) :: y, zpl, zmi, C2pl, C2mi, H
    sigma = 0.d0

    select case (n)
        case (0)
            stop
        case (1)
            sigma = p_excit_1_janev(eb,m_max)
        case (2)
            sigma = p_excit_2_janev(eb,m_max)
        case (3)
            sigma = p_excit_3_janev(eb,m_max)
        case DEFAULT
            nf = real(n)
            m_loop: do m=1,m_max
                if(n.ge.m) then
                    sigma(m) = 0.d0
                    cycle m_loop
                endif
                mf = real(m)
                Etil = Eb/25.0
                s = (mf-nf)

                D = exp(-1.0/(nf*mf*Etil**2.0))
                A = 8.0/(3.0*s)*(mf/(s*nf))**3*(0.184-0.04/s**(2.0/3.0)) * &
                    (1.0 - 0.2*s/(nf*mf))**(1.0 + 2.0*s)
                G = 0.5*( Etil*nf**2.0/(mf - 1.0/mf) )**3.
                L = log(1.0 + 0.53*Etil**2.0 * nf*(mf - 2.0/mf)/(1.0 + 0.4*Etil))
                F = ( 1.0 - 0.3*s*D/(nf*mf) )**(1.0 + 2.0*s)

                y = 1.0/( 1.0 - D*log(18*s)/(4.0*s) )
                zpl = 2.0/(Etil * nf**2 * ( (2.0 - (nf/mf)**2)**0.5 + 1.0))
                zmi = 2.0/(Etil * nf**2 * ( (2.0 - (nf/mf)**2)**0.5 - 1.0))
                C2pl = zpl**2 * log(1.0 + 2.0*zpl/3.0)/(2.0*y + 3.0*zpl/2.0)
                C2mi = zmi**2 * log(1.0 + 2.0*zmi/3.0)/(2.0*y + 3.0*zmi/2.0)

                H = C2mi - C2pl

                sigma(m) = ((8.8d-17*n**4)/Etil)*(A*L*D + F*G*H)
            enddo m_loop
    end select

end function p_excit_n