fsv Function

  function fsv(u)

    !! Integrand function for the hot beam fusion reaction rate
    !! author: P J Knight, CCFE, Culham Science Centre
    !! u : input real : abscissa of integration, = ratio of beam velocity
    !! to the critical velocity
    !! This is the integrand function for the hot beam fusion reaction rate.
    !!     !
    ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

    use constants, only: mproton, echarge

    implicit none

    real(dp) :: fsv

    !  Arguments

    real(dp), intent(in) :: u

    !  Local variables

    real(dp) :: t1,t2,xvc,xvcs

    ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

    t1 = (u**3)/(1.0D0+u**3)

    !  vcritx : critical velocity for electron/ion slowing down of beam ion (m/s)

    xvc = vcritx*u
    xvcs = xvc * xvc * mproton/(echarge * 1000.0D0)
    t2 = sigbmfus(xvcs)

    fsv = t1 * t2

  contains

    ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

    function sigbmfus(vrelsq)

      !! Fusion reaction cross-section
      !! author: P J Knight, CCFE, Culham Science Centre
      !! vrelsq : input real :  square of the speed of the beam ion (keV/amu)
      !! This function evaluates the fusion reaction cross-section as a
      !! function of beam ion velocity (squared).
      !! The functional form of the cross-section is in terms of the equivalent
      !! deuterium energy, i.e. for a tritium beam at 500 keV the energy
      !! used in the cross-section function is 333 keV.
      !!       !
      ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

      implicit none

      real(dp) :: sigbmfus

      !  Arguments

      real(dp), intent(in) :: vrelsq

      !  Local variables

      real(dp) :: a1,a2,a3,a4,a5,atmd,ebm,t1,t2

      ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

      a1 = 45.95D0
      a2 = 5.02D4
      a3 = 1.368D-2
      a4 = 1.076D0
      a5 = 4.09D2

      !  Deuterium atomic mass

      atmd = 2.0D0

      !  Beam kinetic energy

      ebm = 0.5D0 * atmd * vrelsq

      !  Set limits on cross-section at low and high beam energies

      if (ebm < 10.0D0) then
         sigbmfus = 1.0D-27
      else if (ebm > 1.0D4) then
         sigbmfus = 8.0D-26
      else
         t1 = a2/(1.0D0 + (a3 * ebm - a4)**2) + a5
         t2 = ebm * (exp (a1/sqrt(ebm)) - 1.0D0)
         sigbmfus = 1.0D-24 * t1/t2
      end if

    end function sigbmfus

  end function fsv