bosch_hale Function

  function bosch_hale(t,reaction)

    !! Routine to calculate the fusion reaction rate
    !! author: R Kemp, CCFE, Culham Science Centre
    !! author: P J Knight, CCFE, Culham Science Centre
    !! t : input real : Maxwellian density-weighted ion temperature (keV)
    !! reaction : input integer : flag for fusion reaction to use:
    !! 1 : D-T reaction
    !! 2 : D-3He reaction
    !! 3 : D-D 1st reaction (50% probability)
    !! 4 : D-D 2nd reaction (50% probability)
    !! This routine calculates the volumetric fusion reaction rate
    !! <I>&lt;sigma v&gt;</I> in m3/s for one of four nuclear reactions,
    !! using the Bosch-Hale parametrization.
    !! <P>The valid range of the fit is 0.2 keV < t < 100 keV
    !! Bosch and Hale, Nuclear Fusion 32 (1992) 611-631
    !
    ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

    implicit none

    real(dp) :: bosch_hale

    !  Arguments

    real(dp), intent(in) :: t
    integer, intent(in) :: reaction

    !  Local variables

    integer, parameter :: DT=1, DHE3=2, DD1=3, DD2=4
    real(dp) :: theta1, theta, xi
    real(dp), dimension(4) :: bg, mrc2
    real(dp), dimension(4,7) :: cc

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

    if  (t == 0.0D0) then
       bosch_hale = 0.0D0
       return
    end if

    !  Gamov constant, BG

    bg(DT)   = 34.3827D0  !  D + T --> 4He + n reaction
    bg(DHE3) = 68.7508D0  !  D + 3He --> 4He + p reaction
    bg(DD1)  = 31.3970D0  !  D + D --> 3He + n reaction
    bg(DD2)  = 31.3970D0  !  D + D --> T + p reaction

    !  Reduced mass of the particles, keV

    mrc2(DT)   = 1.124656D6
    mrc2(DHE3) = 1.124572D6
    mrc2(DD1)  = 0.937814D6
    mrc2(DD2)  = 0.937814D6

    !  Parametrization coefficients

    cc(DT,1) =  1.17302D-9
    cc(DT,2) =  1.51361D-2
    cc(DT,3) =  7.51886D-2
    cc(DT,4) =  4.60643D-3
    cc(DT,5) =  1.35000D-2
    cc(DT,6) = -1.06750D-4
    cc(DT,7) =  1.36600D-5

    cc(DHE3,1) =  5.51036D-10
    cc(DHE3,2) =  6.41918D-3
    cc(DHE3,3) = -2.02896D-3
    cc(DHE3,4) = -1.91080D-5
    cc(DHE3,5) =  1.35776D-4
    cc(DHE3,6) =  0.00000D0
    cc(DHE3,7) =  0.00000D0

    cc(DD1,1) =  5.43360D-12
    cc(DD1,2) =  5.85778D-3
    cc(DD1,3) =  7.68222D-3
    cc(DD1,4) =  0.00000D0
    cc(DD1,5) = -2.96400D-6
    cc(DD1,6) =  0.00000D0
    cc(DD1,7) =  0.00000D0

    cc(DD2,1) =  5.65718D-12
    cc(DD2,2) =  3.41267D-3
    cc(DD2,3) =  1.99167D-3
    cc(DD2,4) =  0.00000D0
    cc(DD2,5) =  1.05060D-5
    cc(DD2,6) =  0.00000D0
    cc(DD2,7) =  0.00000D0

    theta1 = t*(cc(reaction,2) + t*(cc(reaction,4) + t*cc(reaction,6))) / &
         (1.0D0 + t*(cc(reaction,3) + t*(cc(reaction,5) + t*cc(reaction,7))))
    theta = t/(1.0D0 - theta1)

    xi = ((bg(reaction)**2)/(4.0D0*theta))**0.3333333333D0

    !  Volumetric reaction rate <sigma v> (m3/s)

    bosch_hale = 1.0D-6 * cc(reaction,1) * theta * &
         sqrt( xi/(mrc2(reaction)*t**3) ) * exp(-3.0D0*xi)

  end function bosch_hale