! this function calculates the beam stopping cross section using r.k janev's
! formulas in
!     "penetration of energetic neutral beams into fusion plasma", _nuclear
!     fusion_, vol.29, no.12 (1989), pp.2125-2139.
!
! coefficients for the data statements were not taken from this article, but
! through private communications with boley.  they should be more accurate for
! tftr beam energies.
!
!
!
! this function can be called by : sigmas (e1, mass, ne, te, zf(1-4))
!     e1 (real) is the beam energy in kev
!     mass (real) is the mass of the beam particle in amu
!     ne (real) is the electron density  in 1 / cm 3
!     te (real) is the electron temperature in kev
!
!     zf(1) is the fractional z effective for he
!     zf(2)              "                    c
!     zf(3)              "                    o
!     zf(4)              "                    fe
!
!     a zf of 0 will cause this routine to ignore any effects do to the
!     presence of that impurity.
!
!
!
! ----------------------------------------------------------------------
!
real(fp) function sigmas (e1, mass, ne, te, zf)
  !
  ! declaring input variables
  !
  !============
  use iso_c_binding, only: fp => c_double
  implicit none
  real(fp) :: s1
  !============
  real(fp) :: e1, mass, ne, te              !energy, mass, dens e, e temp
  real(fp) :: zf(4)
  !
  !
  ! declaring variables for routine's use
  !
  integer :: i,j,k,q                   !index variables
  !                                         q reffers to impurity:
  !                                         1 = he
  !                                         2 = c
  !                                         3 = o
  !                                         4 = fe
  !
  !
  real(fp) :: e                            ! energy/mass (e1/mass)
  real(fp) :: lne                          ! ln (beam energy)
  real(fp) :: n1, n(4)                     ! ion densities
  real(fp) :: sion                         ! sum of ion densities
  real(fp) :: lnne                         ! ln (electron density/10**13)
  real(fp) :: lnte                         ! ln (electron temperature)
  real(fp) :: z(4)                         ! data variable for zi
  real(fp) :: ifac                         ! ifac = lne  in s1 and sz(q)
  real(fp) :: jfac                         ! jfac = lnne in s1 and sz(q)
  real(fp) :: kfac                         ! kfac = lnte in s1 and sz(q)
  !
  !
  real(fp) :: a(3,3,2)                     !data matrix for no impurities
  real(fp) :: b(4,2,2,4)                   !data matrix for impurities
  !                                        the 4th index is q specifying which
  !                                        impurity.
  real(fp) :: sz(4)                        !sz(q) to be used in equation (26) and
  !                                        computed by equation (25)
  real(fp) :: bigsum                       !argument in brackets in equation (26)
  !
  !
  !              end of variable declaration
  !------------------------------------------------------------------------------
  !
  !     data for coefficient matrices (values are not from janev.  they are
  !                                    good for within tftr energy ranges  )
  !
  data (((  a(i,j,k),  k = 1, 2),  j = 1, 3),  i = 1, 3)/&
       3.95_fp,     1.60e-02_fp,&
       -3.84e-02_fp,    -5.98e-03_fp,&
       -3.10e-03_fp,    -1.09e-03_fp,&
       3.67e-01_fp,    -2.15e-02_fp,&
       3.07e-02_fp,     1.78e-03_fp,&
       3.16e-03_fp,     3.47e-04_fp,&
       -9.95e-03_fp,     6.19e-04_fp,&
       -2.36e-03_fp,    -1.67e-04_fp,&
       -1.31e-04_fp,    -2.28e-05_fp/
  !
  !
  !     data statemenoats for impurity coeficients
  !       note:there is a conditional statement so that if zeff fractional for
  !            a given impurity is 0, the data will be ignored since it is
  !            useless.
  !
  !     ------he------------------
  !
  data ((( b(i,j,k,1),      k = 1, 2), j = 1, 2), i = 1, 4)/&
       -1.76_fp,     -2.90e-01_fp,&
       5.43e-02_fp,      1.04e-02_fp,&
       7.49e-01_fp,      1.57e-01_fp,&
       -3.74e-02_fp,     -6.70e-03_fp,&
       -7.28e-02_fp,     -2.45e-02_fp,&
       6.11e-03_fp,      1.14e-03_fp,&
       2.05e-03_fp,      1.34e-03_fp,&
       -2.82e-04_fp,     -5.42e-05_fp/
  !
  !--------c----------------------
  !
  data ((( b(i,j,k,2),     k = 1,2),  j = 1, 2),  i = 1, 4)/&
       -1.89e-01_fp,     -3.22e-02_fp,&
       5.43e-02_fp,      3.83e-03_fp,&
       -1.41e-03_fp,      8.98e-03_fp,&
       -3.34e-02_fp,     -1.97e-03_fp,&
       3.10e-02_fp,      7.43e-04_fp,&
       5.08e-03_fp,      1.96e-04_fp,&
       -2.54e-03_fp,     -4.21e-05_fp,&
       -2.20e-04_fp,      8.32e-07_fp/
  !
  !--------o----------------------
  !
  data ((( b(i,j,k,3),    k = 1, 2),  j = 1, 2),  i = 1, 4)/&
       -1.07e-01_fp,     -3.36e-02_fp,&
       4.90e-02_fp,      3.77e-03_fp,&
       -3.72e-02_fp,      1.11e-02_fp,&
       -2.89e-02_fp,     -1.84e-03_fp,&
       3.34e-02_fp,      7.20e-06_fp,&
       4.14e-03_fp,      1.64e-04_fp,&
       -2.50e-03_fp,      1.08e-05_fp,&
       -1.65e-04_fp,      2.45e-06_fp/
  !
  !-----------fe--------------------
  !
  data ((( b(i,j,k,4),     k = 1, 2),  j = 1, 2),  i = 1, 4)/&
       5.46e-02_fp,     -3.89e-02_fp,&
       1.71e-02_fp,     -7.35e-04_fp,&
       -8.97e-02_fp,      2.36e-02_fp,&
       -7.46e-03_fp,      9.94e-04_fp,&
       2.96e-02_fp,     -4.71e-03_fp,&
       2.52e-04_fp,     -3.05e-04_fp,&
       -1.75e-03_fp,      3.63e-04_fp,&
       4.10e-05_fp,      2.37e-05_fp/
  !
  data z/2.0_fp,6.0_fp,8.0_fp,26.0_fp/
  !
  !------------------------------------------------------------
  !     end of data
  !------------------------------------------------------------
  !
  !     compute : e, lne, lnte, lnne
  !
  e = e1/mass
  lne = log(e)
  lnte = log(te)
  lnne = log (ne/1e13)
  !
  do q = 1, 4
     n(q) = ne * zf(q) / z(q)**2
  enddo
     !
     !
     !    ______ sum s1 _________________________________
     !
  s1 = 0
  do i = 1, 3
     do j = 1, 3
        do k = 1, 2
           !
           ifac = lne
           if ((i .eq. 1).and.(lne  .eq. 0.0_fp)) &
                ifac = 1.0_fp
           jfac = lnne
           if ((j .eq. 1).and.(lnne .eq. 0.0_fp)) &
                jfac = 1.0_fp
           kfac = lnte
           if ((k .eq. 1).and.(lnte .eq. 0.0_fp)) &
                kfac = 1.0_fp
           !
           s1 = s1 + (a(i,j,k) * ifac**(i-1) * jfac**(j-1) *&
                kfac**(k-1))
        enddo
     enddo
  enddo
              !
              !
              !     ______sum sz(q)________________________________
              !
  do q = 1, 4
     !
     sz(q) = 0
     if (zf(q) .eq. 0.0_fp) cycle
     !
     do i = 1, 4
        do j = 1, 2
           do k = 1, 2
              !
              ifac = lne
              if ((i .eq. 1).and.(lne  .eq. 0.0_fp)) &
                   ifac = 1.0_fp
              jfac = lnne
              if ((j .eq. 1).and.(lnne .eq. 0.0_fp)) &
                   jfac = 1.0_fp
              kfac = lnte
              if ((k .eq. 1).and.(lnte .eq. 0.0_fp)) &
                   kfac = 1.0_fp
              !
              sz(q) = sz(q) + (b(i,j,k,q) * ifac**(i-1) *&
                   jfac**(j-1) * kfac**(k-1))
              !
           enddo
        enddo
     enddo
  enddo                          !
  !
                          !     ______sum bigsum________________________________
                          !
  bigsum = 0
  do q = 1, 4
     bigsum = bigsum + (n(q) * z(q) * (z(q) - 1) * sz(q))
  enddo
                             !
                             !
  !     ______compute sigmas____________________________
                             !
  sigmas = exp(s1)/e   *   (1 + bigsum/ne)    *    1.0e-16_fp
  !
  return
end function sigmas
