subroutine coulog(et,at,zt,denb,tempb,ab,zb,cln,nb,Bmag,
     >                    llstix)
C
C (See the warnings below about some (usually small) uncertainties in
C certain parameter regimes.)
C
C calculates the Coulomb logarithm for a test particle of energy et
C (kev), atomic mass at and atomic charge zt, colliding with NB
C maxwellian background species with densities denb (/m**3),
C temperatures tempb (keV), atomic masses ab, atomic charges zb.  The
C electrons should be included as one of the background species, with
C Z=-1 and A=(1./1836.1).  The list of Coulomb logarithms is returned
C in cln.  The magnetic field should be input as bmag (Tesla).
C
C Note that the test particle loglambda's produced by this subroutine
c are not symmetric, i.e. it is not true that
c loglam_alpha_beta=loglam_beta_alpha.  If such symmetry is desired,
c one should set Et=(3/2)T_t and call this routine twice with species
c alpha and beta switched, and take the average value of the log_lambda's
c that are returned.
C
C Written by Greg Hammett, 16-May-1989.
C
C The exact details of the Coulomb logarithm are probably not worth
C worrying about too much.  It has the form
C
C    log(Lambda) = log(r_max/r_min)
C
C For typical values of log(Lambda) ~ 18, a factor of 2 difference in 
C r_max/r_min gives only a 4% correction to log(Lambda).  Since the whole
C assumption of small-angle scattering in the Fokker-Planck approximation
C to the collision operator is only accurate to ~O(1/log(Lambda)), this is
C good enough for most purposes.  Often the most important correction is
C the difference between the classical and quantum r_min, which can give about
C an 18% difference in log(Lambda) at T_e=10 keV and n_e=1.e20/m^3 (using the
C simplified thermal expressions in the NRL formulary).
C
C The recipe I use here is similiar to the
C one given in the NRL plasma formulary.  So it includes the quantum
C mechanics correction which is usually important for collisions with
C electrons.  Rather than completely ignoring the Debye shielding by
C particles with vthermal < v, it contains a smooth transition which I
C have rigorously shown (so I have thought) from the Balescu-Lenard
C to be correct in the limits vthermal>>v and vthermal<<v.
C
C In 1989 I wrote the following, which I now think overstates an issue:
C (However, the Balescu-Lenard
C operator predicts a substantial amount of Cerenkov radiation of
C plasma waves if the particle speed exceeds the electron thermal
C speed.  I am ignoring this Cerenkov radiation (and its associated
C drag).  To look at it properly, Nat Fisch claims that one would have
C to consider the reabsorption of the Cerenkov radiation as well.)
C
C I think it was first in 1991 that I found more accurate ways of
C evaluating the integrals in the Balescu-Lenard operator for energetic
C particles, leading to the comment below:
C
C Warning (usually small): GWH 3/1/91 (updated 4/7/2016): 
C
C The formulas I am using here reduce Debye shielding somewhat and thus
C increase the drag rate somewhat (via a slow logarithmic factor) if the
C test particle is much faster than the *thermal electrons* (like an
C electron tail).  This corresponds to a (usually small) enhancement of
C the effective log(Lambda) factor by Cerenkov radiation.  This fast
C particle limit is not discussed much in the literature or in textbooks,
C but I think the formulas here are at least roughly okay (usually good to
C a few percent), and are fairly consistent with what one gets from the
C Balescu-Lenard collision operator for high velocity particles colliding
C with thermal particles.  At one time I thought that the Balescu-Lenard
C operator blew up in this limit, but I later found how to do the
C complicated asymptotics of the integrals in this regime properly, and
C it asymptotically has the form given here (if I did the calculations
C right).  Comparing with the Perkins 1965 article cited below also
C supports the formulas here, for the case of very fast ions (faster
C than thermal electrons) slowing down on thermal electrons.  (Compare
C with Eqs. 26 and 28 in Perkins 1965.)  The formulas here also agree
C with Perkins 1965 on the energy drag rate for fast electrons slowing
C down on thermal electrons, but there are some differences for the
C electron drag rate (the rate of loss of directed momentum, sometimes
C called dynamical friction), and for slowing down of fast ions or
C electrons on thermal ions (but for high velocity ions, drag on thermal
C ions is weak compared to drag on thermal electrons anyway).
C
C The biggest difference between the formulas here and in Perkins might
C be for the dynamical friction drag rate for runaway electrons in a
C cold plasma.  The formula here says
C
C    b_max ~ (v/v_t) b_max_thermal, 
C
C while Perkins Eq. 28 appears to have
C
C    b_max ~ (v/v_t)^(1/2) b_max_thermal.
C
C So the Perkins formula would reduce the log(Lambda) here by 
C log((v_t/v)^(1/2)).  For a 100 keV electron colliding with 10 eV
C thermal electrons, that would reduce ln(Lambda) by about 2.3, or
C about 20% at ln(Lambda) ~ 10.
C
C End of 2016 comment/warning.

C
C The main references for this stuff are D.V. Sivukhin, Reviews of
C Plasma Physics, Vol. 4 (1966) (especially good on quantum corrections,
C physical insight), Trubnikov, Reviews of Plasma Physics, Vol. 1
C (1963) (simple, complete derivation of Fokker-Planck equation,
C complete with justification for the Coulomb Logarithm in terms of
C deriving the collision operator from the dielectric response of the
C plasma), Rob Goldston's Ph.D. thesis (references on strong magnetic
C field corrections, collisions with multi-electron impurities),
C Krommes' class notes on the derivation of the Landau operator from the
C Balescu-Lenard operator, and my notes on Krommes' notes relating to
C the problem when the test particle is much faster than the thermal
C electron velocity.
C
	real et,at,zt,Bmag
	real denb(nb),tempb(nb),ab(nb),zb(nb),cln(nb)
 
	Logical llstix       ! tbt 7/10/95
 
c
c first calculate the maximum impact parameter, rmax:
c
c rmax=( sum_over_species omega2/vrel2 )**(-1/2)
c
c where,
c omega2 = omega_p**2 + omega_c**2
c vrel2=T/m+2.*E_t/m_t
c
	sum=0.0
	do 100 i=1,nb
 	omega2=1.74*zb(i)**2/ab(i)*denb(i)
	1 +9.18e15*zb(i)**2/ab(i)**2*Bmag**2
	vrel2=9.58e10*(tempb(i)/ab(i) + 2.*et/at)
100	sum=sum+omega2/vrel2
	rmax=sqrt(1./sum)
 
c next calculate rmin, including quantum corrections.  The classical
C rmin is:
c
C rmincl = e_alpha e_beta / (m_ab vrel**2)
C
C where m_ab = m_a m_b / (m_a+m_b) is the reduced mass.
c vrel**2 = 3 T_b/m_b + 2 E_a / m_a
c (Note:  the two different definitions of vrel2 used in this code
c are each correct for their application.)
c
C The quantum rmin is:
C
C rminqu = hbar/( 2 exp(0.5) m_ab vrel)
C
C and the proper rmin is the larger of rmincl and rminqu
C
	do 200 i=1,nb
	vrel2=9.58e10*(3*tempb(i)/ab(i)+2*et/at)
	rmincl=0.13793*abs(zb(i)*zt)*(ab(i)+at)/ab(i)/at/vrel2
	rminqu=1.9121e-8*(ab(i)+at)/ab(i)/at/sqrt(vrel2)
	rmin=max(rmincl,rminqu)
	cln(i)=log(rmax/rmin)
	if(cln(i) .lt. 1.0) then
	    write(6,*) 'warning from COULOG: coulomb logarithms < 1!'
	    cln(i)=1.0
	endif
 
200	continue

! Another 2016 comment:
!
! Comparing with Perkins, 1965,  Phys. Fluids 8, 1361, for v>>v_te, 
! and focussing on just the log(Lambda) for the fast particle energy
! loss rate on thermal electrons (the g_1(v) term in his Eq. 26)
! suggests that the above should be modified from:
!
!    ln(Lambda) = log(r_max/r_min)
!
! to
!
!    ln(Lambda) = ln(exp(0.5)/1.467*r_max/r_min)
!               = ln(1.12*r_max/r_min)
!
! This is only a very minor correction!!! (If I did the calculation right.)
! One could eventually also check other regimes and thus the Pade
! interpolation formulas. 
!
! Could further modify the above formula to be
!
!    ln(Lambda) => ln(sqrt(1+(1.12*r_max/r_min)**2))
! 
! This can be derived from the calculation for the dynamical friction
! drag rate, and gives a positive value for ln(Lambda) for all possible
! r_max/r_min, but it is not correct compared to molecular dynamics
! simulations for small r_max/r_min, where certain correlation effects
! become important.
!
! Expanding for large Lambda, we get
!
!    ln(sqrt(1+(C*Lambda)^2)) = ln(C*Lambda)*ln(sqrt(1+1/(C*Lambda)^2))
!                             ~ ln(C*Lambda) + 0.5/(C*Lambda)^2
! 
! and it is known that there are other corrections of order
! 1/sqrt(Lambda_eff) (see Perkins 1965, if I understand him right).
! Perkins calculates the first order correction C, but no one (I think)
! has tried to go beyond that.  These other terms would dominate over
! the above 0.5/Lambda_eff^2 term.  Also, although one could in 
! principle calculate the dynamical drag rate or the mean energy loss
! rate to arbitrary order in 1/Lambda, the whole assumption that
! small-angle scattering dominates breaks down so that a Fokker-Planck
! form for the collision operator is no longer valid and one should
! use a Boltzmann collision operator and include large-angle
! scattering.
!
! See papers by Baalrud et al. for a better transition to the strongly coupled 
! regime.  An improved fit is if r_max is limited to be no smaller than
! 1/n^(1/3), the average interparticle spacing.  But I think even that does 
! not do as well as empirical fits to molecular dynamics (MD) simulations of
! the strongly-coupled regime.

c debug section:  set all log(lambda)'s to a simple formula for
c log(lambda_e) to benchmark with Stix's solutions:
 
	if(llstix) then
	do 300 i=1,nb
300	cln(i)=24.0-alog(sqrt(denb(1)/1.e6)/(tempb(1)*1000.))	
	endif
 
	return
	end