neoclassics_calc_nu Function

    function neoclassics_calc_nu()
        !! Calculates the collision frequency
        !
        ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
        use const_and_precisions, only: pi, me_, mp_, eps0_,e_

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

        integer :: jj,ii,kk

        temp = temperatures
        density = densities

        !          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(:,:) = 0.0

        do jj = 1, 4
           do ii = 1, no_roots
              x = roots(ii)
              do kk = 1,4
                 xk = (mass(kk)/mass(jj))*(temp(jj)/temp(kk))*x
                 expxk = exp(-xk)
                 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)
                 v = sqrt(2.0d0*x*temp(jj)/mass(jj))
                 neoclassics_calc_nu(jj,ii) = neoclassics_calc_nu(jj,ii) + density(kk)*(z(jj)*z(kk))**2 &
                              *lnlambda *phixmgx/(4.0*pi*eps0_**2*mass(jj)**2*v**3)
              enddo
           enddo
        enddo

    end function neoclassics_calc_nu