!/ module elossr
!  version 1.0 26-nov-91 from mios:ionbem.for(setsig),bounce.for(aioniz,fsig)
!
!.              %%%%%%%%%%%%%%%%%%%%%%%%%
!.              %                       %
!.              %       elossr          %
!.              %                       %
!.              %%%%%%%%%%%%%%%%%%%%%%%%%
!
function elossr(pebeam,pabeam,kznuc,kzion, &
     kspec,pdens,ptemp,paspec,pzspec, kerror)
  !  calculates the electron loss rate (in /sec) for a beam atom (or ion)
  !  summed over a number of thermal species.
  !
  !. begin index of input arguments:
  !r pebeam            energy of beam particle (ev).
  !r pabeam            atomic weight of beam particle (amu).
  !i kznuc             nuclear charge of beam particle (1,2,...).
  !i kzion             ionization state of beam particle (0,1,...,kznuc-1).
  !i kspec             number of background thermal species.
  !r pdens(kspec)      background species particle density (1/cm^3).
  !r ptemp(kspec)      background species temperature (ev).
  !r paspec(kspec)     background species atomic weight (amu).
  !r pzspec(kspec)     background species charge (-1,0,1,2,... noninteger :: is ok).
  !.    also used to determine identity of species (electron, impurity)
  !. output argument
  !i kerror            error code for problematic input:
  !.         1   kznuc unusable
  !.         2   kzion impossible
  !.         3   kspec .le. 0
  !. end of variable index.
  !. begin code generated by vndx.
  !============
  use iso_c_binding, only: fp => c_double
  use physconst_mod, only : kb_ev, amucgs=>amu_g !amucgs is one atomic mass unit in grams
  implicit none
  real(fp), parameter :: evs = kb_ev*1.0e+7_fp  !evs is one ev in ergs.
  real(fp), parameter :: evsinv = 1.0e-7_fp/kb_ev !evsinv is one erg in ev.
  integer :: mint,mgl,igl,katom,jspec,ipel
  !============
  real(fp) :: elossr,zvrel,zvrel2,zvsig,zv,zsvei
  real(fp) :: zsvpi,zsvcx,zelimp,zteev,aioniz,zsvavg,sgvbt1,zzi,zefq
  !============
  integer :: kznuc,kzion,kspec,kerror,int
  real(fp) :: pebeam,pabeam,pdens(kspec),ptemp(kspec),paspec(kspec), &
       pzspec(kspec)
  !. end code generated by vndx.
  parameter (mint=5, mgl=8)
  data int,igl/ mint, mgl/
  dimension zvrel(mint*mgl),zvrel2(mint*mgl),zvsig(mint*mgl)
  logical ilnucl(54)
  !  ilnucl is .true. for z=1-6, 8, 10, 11, 12, 14, 18, 36, 54.
  data ilnucl/6*.true., .false., .true.,  .false., .true., &
       .true.,  .true., .false., .true., 3*.false., .true., &
       17*.false., .true., 17*.false., .true./
!
  elossr=0._fp
  !  check the input for fatal problems
  if(kznuc.lt.1 .or. kznuc.gt.54 .or. .not.ilnucl(kznuc)) then
     kerror=1
     return
  end if
  if(kzion.lt.0 .or. kzion.gt.kznuc-1) then
     kerror=2
     return
  end if
  if(kspec.le.0) then
     kerror=3
     return
  end if
  !  the speed of the beam particles is used with cold targets.
  zv=sqrt(2._fp*evs*pebeam/(pabeam*amucgs))
  katom=kzion+1
  !  sum the contributions from each background species.
  do jspec=1,kspec
     if(pzspec(jspec).eq.-1._fp) then
        !          electron ionization.
        zteev=ptemp(jspec)
        zsvei=aioniz(kznuc,katom,zteev,pabeam,3)
        zsvavg=zsvei
     elseif(pzspec(jspec).eq.1._fp) then
        !          proton ionization and charge exchange.
        !          averages over a maxwellian target distribution.
        if(kznuc.eq.1) then
           !            hydrogen beam so include charge exchange with background.
           zsvpi=sgvbt1(pebeam,pabeam,ptemp(jspec),paspec(jspec), &
                kznuc,katom,1)
           zsvcx=sgvbt1(pebeam,pabeam,ptemp(jspec),paspec(jspec), &
                kznuc,katom,4)
           zsvavg=zsvpi+zsvcx
        else
           !            non-hydrogenic beam so ignore charge exchange.
           zsvpi=sgvbt1(pebeam,pabeam,ptemp(jspec),paspec(jspec), &
                kznuc,katom,1)
           zsvavg=zsvpi
           zsvcx=0._fp
        end if
     elseif(pzspec(jspec).gt.1._fp) then
        !          impurity collision electron loss (ionization and charge transfer).
        zzi=pzspec(jspec)
        zefq=pebeam/zzi
        !          normally use olson and salop universal formula:
        ipel=5
        !          for very high beam energies use gryzinski formula:
        if(zefq.gt.pabeam*2.e5_fp) ipel=1
        zelimp=zzi*aioniz(kznuc,katom,zefq,pabeam,ipel)
        zsvavg=zv*zelimp
     else
        zsvavg=0._fp
     end if
     elossr=elossr+pdens(jspec)*zsvavg
  end do
  return
end function elossr
