neoclassics_calc_nu_star_fromT Function

    function neoclassics_calc_nu_star_fromT(iota)
        !! Calculates the collision frequency
        !
        ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
        use const_and_precisions, only: pi, me_, mp_, eps0_,e_, keV_
        use physics_variables, only: rmajor, te,ti, dene,deni, dnalp, fdeut

        real(dp),dimension(4) :: neoclassics_calc_nu_star_fromT
        real(dp) :: t,erfn,phixmgx,expxk,xk, lnlambda,x,v
        real(dp),dimension(4) :: temp, mass,density,z
        real(dp) :: iota

        integer :: jj,kk

        temp = (/te,ti,ti,ti /) * keV_
        density = (/dene,deni * fdeut,deni*(1-fdeut),dnalp /)

        !          e      D      T         a (He)
        mass = (/me_,mp_*2.0d0,mp_*3.0d0,mp_*4.0d0/)
        z = (/-1.0d0,1.0d0,1.0d0,2.0d0/) * e_

        ! transform the temperature back in eV
        ! Formula from L. Spitzer.Physics of fully ionized gases.  Interscience, New York, 1962
        lnlambda = 32.2d0 - 1.15d0*log10(density(1)) + 2.3d0*log10(temp(1)/e_)

        neoclassics_calc_nu_star_fromT(:) = 0.0d0

        do jj = 1, 4
            v = sqrt(2d0 * temp(jj)/mass(jj))
            do kk = 1,4
                xk = (mass(kk)/mass(jj))*(temp(jj)/temp(kk))

                if (xk < 200.d0) then
                    expxk = exp(-xk)
                else
                    expxk = 0.0d0
                endif

                t = 1.0d0/(1.0d0+0.3275911d0*sqrt(xk))
                erfn = 1.0d0-t*(.254829592d0 + t*(-.284496736d0 + t*(1.421413741d0       &
                        + t*(-1.453152027d0 +t*1.061405429d0))))*expxk
                phixmgx = (1.0d0-0.5d0/xk)*erfn + expxk/sqrt(pi*xk)
                neoclassics_calc_nu_star_fromT(jj) = neoclassics_calc_nu_star_fromT(jj) + density(kk)*(z(jj)*z(kk))**2 &
                            *lnlambda*phixmgx/(4.0d0*pi*eps0_**2*mass(jj)**2*v**4) * rmajor/iota
            enddo
        enddo

    end function neoclassics_calc_nu_star_fromT