!----------------------------------------------------------------------
!  sigfacx -- quick estimate of excited states cross section
!             enhancement factor
!
!  use the janev/boley/post fit as coded in sigmas.for.  this fit
!  returns the cross section (cm**2) itself;  here i take the fit
!  cross-section at ne=1e12 / a small correction factor to define
!  the "ground-state-cross-section".
!
!  then
!                         total stopping cross section
!    enhancement =      --------------------------------
!                          ground state cross section
!
!  the small correction factor is deduced by reading the value at 1e12
!  of the plots of "delta" on page 2134 (right hand column) of the
!  janev/boley/post article, nucl. fusion vol.29, no. 12 (1989).  the
!  dependence of the correction factor on zeff is small and is ignored;
!  the correction factor is taken to depend on energy only and ranges
!  from 1.0 (no correction at low energy) to 1.1 (at 1000 kev/amu).
!
!  the fit used in the subfunction sigmas.for has been extended to
!  be valid btw 10 kev/amu and 100 kev/amu -- viz. priv. communication
!  w/boley.
!
real(fp) function sigfacx(e1ev, amu, ne, te_ev, zf, lambdgs, &
     lambdexc)
  !
  !============
  use iso_c_binding, only: fp => c_double
  implicit none
  integer :: inumeb,iset,i
  !============
  real(fp) :: enorm,fac2,fact,faci,zsigs,sigmas,zsign,zfacx
  !============
  real(fp) :: e1ev, amu, ne, te_ev, zf(4)
  !
  !  the first five arguments are as in "sigmas.for" except that e1ev and
  !  te_ev are in ev instead of kev.
  !
  !  the last two arguments are output:
  !
  real(fp) :: lambdgs,lambdexc
  !
  !    lambdgs = estimated mean free path (ground state model)
  !    lambdexc= estimated mean free path (accounting for excited states)
  !
  !  (caution-- if e1ev/amu < 10kev/amu, the lambd* outputs are less accurate)
  !
  !------------------------------
  !  correction factor at 1e12 density ...
  !
  parameter (inumeb=7)
  real(fp) :: ebase(inumeb),cfact(inumeb),cflog(inumeb-1)
  !
  data ebase/ &
       10.0e3_fp,50.0e3_fp,200.0e3_fp,500.0e3_fp, &
       1000.0e3_fp,2000.0e3_fp,10.0e6_fp/
  data cfact/ &
       1.00_fp, 1.04_fp, 1.06_fp, 1.08_fp, 1.10_fp, 1.12_fp, 1.16_fp/
  !
  data iset/0/
  !
  save cflog,iset
  !
  !------------------------------
  !  do once on first call ...
  !
  if(iset.eq.0) then
     iset=1
     do i=1,inumeb-1
        cflog(i)=log(ebase(i+1)/ebase(i))
     end do
  end if
  !
  !  interpolate correction factor
  !
  enorm=e1ev/amu
  fac2=1.0_fp
  if(enorm.le.ebase(1)) then
     fact=cfact(1)
     fac2=max(0.0_fp,enorm)/ebase(1)
     enorm=ebase(1)
  else if(enorm.ge.ebase(inumeb)) then
     fact=cfact(inumeb)
  else
     do i=1,inumeb-1
        if((ebase(i).le.enorm).and.(enorm.le.ebase(i+1))) exit
     end do
     faci = log(enorm/ebase(i))/cflog(i)
     fact=cfact(i)+faci*(cfact(i+1)-cfact(i))
  end if
  !
  if(ne.lt.1.0e12_fp) then
     fact=1.0_fp+(fact-1.0_fp)*max(0.0_fp,ne)/1.0e12_fp
  end if
  !
  !  evaluate enhancement factor
  !
  zsigs=sigmas(enorm*amu*0.001_fp, amu, &
       max(1.0e12_fp,ne), te_ev*0.001_fp,zf)
  zsign=sigmas(enorm*amu*0.001_fp, amu, &
       1.0e12_fp, te_ev*0.001_fp,zf)
  !
  zfacx= (fact * zsigs)/zsign
  if(fac2.lt.1.0_fp) then
     zfacx=1.0_fp+(zfacx-1.0_fp)*fac2
  end if
  !
  lambdexc=1/(ne*zsigs)
  lambdgs=fact/(ne*zsign)
  if(fac2.lt.1.0_fp) then
     lambdexc=sqrt(fac2)
     lambdgs=sqrt(fac2)
  end if
  !
  sigfacx = zfacx
  !
  return
end function sigfacx
