c
c###################################################################
c
      subroutine pexch4(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 excitation
c     in the dipole-approximation close-coupling method (dacc).
c
c     pe = collision energy in ev
c
c     the coefficient data array passed should contain
c
c     pcf(1)  =   zt, the atomic number of target
c
c     pcf(2)  =   zp, the atomic number of projectile
c
c     pcf(3)  =  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(4) =   n, the principal quantum number of the initial state
c
c     pcf(5) =   m, the principal quantum number of the final state
c
c     if itrans=2, the remainder of the coefficient array should
c                 contain
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
c                total spin is returned
c
c
c     pcf(7) =   lin, the orbital angular monentum of the inital state
c
c     pcf(8) =   lfin, the orbital angular monentum of the final state
c
c     pcf(9) =   mult, spin multiplicity (2s+1) of the inital state
c
c     pcf(10) =   j, the total angular momentum of the initial state
c                only required for fine structure transitions
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(11) = eth,threshold excitation energy for the transition (ev)
c
c     pcf(12) = wnm, energy difference between initial and final states
c
c     pcf(13) = parameter lambda
c
c     kermsg = blank if no errors
c
c     pxs = cross section in cm[2]
c
c      this evaluation function can only derive cross sections for those
c    transitions for which the oscillator strengths have been included
c    in subroutines oscsth1 and oscthe. oscsth1 contains the oscillator
c    strengths for allowed transitions for hydrogen from table a.1. of
c    janev et al. oscthe contains the oscillator strengths for allowed
c    transitions for helium  from table a.4 of janev et al.
c      for transitions beteen excited states of helium the excitation
c    energies are returned by subroutine heexcen and the ionization
c    energies are returned by subroutine heionen.
c
c------------------------------------------------------------------------
c
      double precision pe, pcf, pxs
c
      integer zt, zp, n, m, sumen, itrans
      real lambda
      dimension pcf(13)
      character*(*) kermsg
c
      data ry /1.358e+01/
c
      zt = pcf(1)
      zp = pcf(2)
      rzt = zt
      rzp = zp
      itrans  = pcf(3)
      n = pcf(4)
      m = pcf(5)
      an = n
      am = m
c
      if (m .le. n) then
        kermsg = 'error  principal quantum number m <= n '
        return
      endif
c
      if (kncf .lt. 11) then
c
c        first call to pexch4 determine energy independent
c        parameters and place in pcf for further use
c
c        (n,m) transitions - determine energies, wnm and oscillator
c                            strength.
c
        if (itrans .eq. 1) then
          if (zt .eq. 1) then
            c= 1.0/(an*an) - 1.0/(am*am)
            eth = ry * c
            wnm = 0.5 * c
          else if (zt .eq. 2) then
            call heionen(n, 0, 0, 1, eionn, kermsg)
            call heionen(m, 0, 0, 1, eionm, kermsg)
            eth = eionn - eionm
            wnm = eth /(2.0*ry)
          else
            kermsg = 'target not covered by pexch4 '
            return
          endif
c
c     calculate the oscillator strength, fnm, for the transition from
c     the formula of johnson
c
          if (zt .eq. 1) then
            call oscsth2(n, m, fnm)
          else
            kermsg = 'average oscillator strength missing in pexch4'
            return
          endif
c
c        (nl,ml') transitions - determine energies, wnm and oscillator
c                               strength.
c
        else if (itrans .eq. 2 .and. zt .eq. 2) then
          if (kncf .ne. 10) then
            kermsg =
     1           ' invalid number of coefficients passed to pexch4'
            return
          endif
c
          sumen =  pcf(6)
          lin =  pcf(7)
          lfin = pcf(8)
          mult = pcf(9)
          j = pcf(10)
c
          if (n .eq .1) then
             call heexcen(m, lfin, mult, sumen, eth, kermsg)
             wnm = eth /(2.0*ry)
          else
              call heionen(n, lin, mult, sumen, ethn, kermsg)
              call heionen(m, lfin, mult, sumen, ethm, kermsg)
              eth = ethn - ethm
              wnm = eth /(2.0*ry)
          endif
        endif
c
        if (itrans .eq. 2) then
          call oscsthe(n, lin, m, lfin, mult, j, fnm, kermsg)
        endif
c
c
c     determine the potential strength, lambda
c
      lambda= sqrt(fnm/(2.0*wnm))
c
c        place energy independent parameters in coefficient array and
c        update kncf
c
        pcf(11) = eth
        pcf(12) = wnm
        pcf(13) = lambda
        kncf = 13
      else if (kncf .eq. 13) then
          eth  = pcf(11)
          wnm  = pcf(12)
          lambda = pcf(13)
c
      else
          kermsg = ' incorrect number of coefficients passed to pexch4'
          return
      endif
c
      if(pe .lt. eth) then
        pxs=0.0
        return
      endif
c
c     determine the velocity, v
c     corrected 1/zp factor -- dh -- 13 March 2008
        v = 6.3246e-03* sqrt(pe) / zp
c
c     determine the value of beta
c
      beta = zp* lambda *  wnm / (v * v)
c
c     determine the value of the cross section  pxs
c
      pxs = 1.76e-16 *  zp * lambda * dbeta(beta) / wnm
c
      return
      end