c
c###################################################################
c
      subroutine ionbea(pe, pcf, kncf, pxs, kermsg)
c
c     this is a subroutine to calculate cross sections (cm[2])
c     versus energy (ev) for heavy particle impact ionization
c     in the binary encounter approximation (bea). this routine is
c     only valid for a hydrogenic core.
c
c     pe = collision energy in ev
c
c     the coefficient data array passed should contain
c
c     pcf(1)  =   ant , the atomic number of target
c
c     pcf(2)  =   zct, the charge of the core of the target
c
c     pcf(3)  =   zp, the charge of the projectile ion
c
c     pcf(4) =   eion, binding energy of target atom. if this value is
c                zero and z0 = 2 (he) the binding energy is determined
c                from the routine heionen.
c
c     pcf(5)  =  itrans,  integer which defines the type of transition.
c                 for transitions defined only in terms of the initial
c                 and final principal quantum numbers (n,m), itrans=1.
c                 for tranitions between (nl,ml') states itrans=2
c
c     pcf(6) =   sumen, this indicates choice of ionization energy to be
c                returned.  if sumen=0  the ionization energy for the
c                specific state (quantified by n, l amd mult) is
c                returned.  if sumen=1  the ionization energy taken as
c                an average over angular momentum and total spin is
c                returned. this only applies to helium as a target.
c
c     pcf(7) =   n, the principal quantum number of the initial state
c
c     pcf(8) =   lin, the orbital angular monentum of the inital state
c
c     pcf(9) =   mult, spin multiplicity of the initial state
c
c    - warning- .
c
c        the coefficient array pcf is updated by this routine to
c     include energy independent constants. these coefficients can be
c     used in subsequent calls for the same entry. the coefficeients
c     added are:
c
c     pcf(10) = eth,threshold enegry for ionization
c
c     pcf(11) = neff, effective value of the principal quantum number
c               used. (see doc=h-he-plasma for definition)
c
c     kermsg = blank if no errors
c
c     pxs = cross section in cm[2]
c
c
c------------------------------------------------------------------------
c
      double precision pe, pcf, pxs
c
      integer ant, n, zct, zp, itrans, sumen, lin, mult
      real neff
      dimension pcf(11)
      character*(*) kermsg
c
      data ry /1.358e+01/
c
      ant  = pcf(1)
      zct  = pcf(2)
      zp   = pcf(3)
      eion = pcf(4)
      itrans = pcf(5)
      sumen  = pcf(6)
c
      if (itrans .ne. 1 .and. itrans .ne. 2) then
        kermsg = ' error invalid fifth coefficient passed to #ionbea'
        return
      endif
c
      if (sumen .ne. 0 .and. sumen .ne. 1) then
        kermsg = ' error invalid sixth coefficient passed to #ionbea'
        return
      endif
c
      n = pcf(7)
      an = n
      kermsg = ' '
      if (kncf .le. 9) then
c
c        first call to ionbea determine energy independent
c        parameters and place in pcf for further use
c
c        determine the binding threshold energy (eion)
c
        if (eion .gt. 0.0 ) then
            enion = eion
            neff=n
        else if (itrans .eq. 1) then
          if (ant .eq. 1) then
            enion = ry / (an*an)
            neff=n
          else if (ant .eq. 2) then
            call heionen(n, 0, 0, 1, enion, kermsg)
            neff = sqrt(ry/enion)
          else
            kermsg = 'target not covered by ionbea '
            return
          endif
c
        else if (itrans .eq. 2 .and. ant .eq. 2) then
          if (kncf .ne. 9) then
            kermsg = ' invalid number of coefficients passed to ionbea'
            return
          endif
c
          lin = pcf(8)
          mult = pcf(9)
          call heionen(n, lin, mult, sumen, enion, kermsg)
          neff = sqrt(ry/enion)
        endif
c
c        place energy independent parameters in coefficient array and
c        update kncf
c
        pcf(10) = enion
        pcf(11) = neff
        kncf = 11
      else if (kncf .eq. 11) then
          enion = pcf(10)
          neff = pcf(11)
c
      else
          kermsg = ' incorrect number of coefficients passed to ionbea'
          return
      endif
c
      if(pe .lt. enion) then
        pxs=0.0
        return
      endif
c
c     determine alpha
c
      alpha= (6.3246e-03/zp) * neff * dsqrt (pe)
c
c     determine from the value of alpha which formula to use to
c     calculate the cross section
c
      alpsq = alpha  * alpha
      if (alpha .lt. 0.207) then
         pxs = 0.0
      else if (alpha .lt. 1.207) then
         pxs = 5.867e-17 * (zp**2) * (neff**4) *
     1         ( alpha - (0.164/alpsq) + 0.1875/(alpsq*(alpha+ 1.0)))
        else
           pxs = 1.467e-16 * (zp**2) * (neff**4) *
     1                        ( 1.0 - (0.15/(alpsq-1.0)))/alpsq
      endif
c
      return
      end