bt_maxwellian_eb Subroutine

  1. atomic_tables.f90
  2. atomic_tables
  3. bt_maxwellian_eb

public subroutine bt_maxwellian_eb(fn, T, eb, am, ab, rate)

Calculates Maxwellian reaction rate for a beam with atomic mass ab and energy eb firing into a target with atomic mass am and temperature T which has a cross section given by the function fn

Arguments

Type IntentOptional AttributesName
public function fn(a)

Cross section function

Arguments
Type IntentOptional AttributesName
real(kind=8), intent(in) :: a
Return Value real(kind=8)
real(kind=Float64), intent(in) :: T

Target temperature [keV]
real(kind=Float64), intent(in) :: eb

Beam energy [keV]
real(kind=Float64), intent(in) :: am

Target atomic mass [amu]
real(kind=Float64), intent(in) :: ab

Beam atomic mass [amu]
real(kind=Float64), intent(out) :: rate

Reaction Rate []

Calls

proc~~bt_maxwellian_eb~~CallsGraph proc~bt_maxwellian_eb bt_maxwellian_eb proc~simpsons_rule simpsons_rule proc~bt_maxwellian_eb->proc~simpsons_rule

Called by

proc~~bt_maxwellian_eb~~CalledByGraph proc~bt_maxwellian_eb bt_maxwellian_eb interface~bt_maxwellian bt_maxwellian interface~bt_maxwellian->proc~bt_maxwellian_eb

Contents

Source Code

bt_maxwellian_eb

Source Code

subroutine bt_maxwellian_eb(fn, T, eb, am, ab, rate)
    !+ Calculates Maxwellian reaction rate for a beam with atomic mass `ab` and energy `eb`
    !+firing into a target with atomic mass `am` and temperature `T` which has a cross section given by the function `fn`
    interface
        function fn(a)
            !+Cross section function
            real(8)              :: fn !sigma
            real(8), intent(in)  :: a !eb
        end function fn
    end interface
    real(Float64), intent(in)  :: T
        !+Target temperature [keV]
    real(Float64), intent(in)  :: eb
        !+Beam energy [keV]
    real(Float64), intent(in)  :: am
        !+Target atomic mass [amu]
    real(Float64), intent(in)  :: ab
        !+Beam atomic mass [amu]
    real(Float64), intent(out) :: rate
        !+Reaction Rate [\(cm^3/s\)]

    integer, parameter :: n_vr = 30
    real(Float64) :: vr_max, dvr
    real(Float64), dimension(n_vr) :: vr
    real(Float64), dimension(n_vr) :: fr
    integer, parameter :: n_vz = 60
    real(Float64) :: vz_max,dvz
    real(Float64), dimension(n_vz) :: vz
    real(Float64), dimension(n_vz) :: fz
    real(Float64) :: T_per_amu, eb_per_amu, ared, sig, sig_eff
    real(Float64) :: zb, u2_to_erel, u2, erel, v_therm, dE

    integer :: i, j

    vr_max = 4.d0
    dvr = vr_max/(n_vr - 1.d0)
    do i=1,n_vr
        vr(i) = (i-1)*dvr
    enddo

    vz_max = 4.d0
    dvz = 2.0*vz_max/(n_vz - 1.d0)
    do i=1,n_vz
        vz(i) = (i-1)*dvz - vz_max
    enddo

    T_per_amu = max(T, 1.d-6)/am
    eb_per_amu = eb/ab
    ared = am*ab/(am + ab)

    v_therm = 1.384d6 * sqrt(T_per_amu*1.d3)
    zb = sqrt(eb_per_amu/T_per_amu)
    u2_to_erel = ared*T_per_amu

    fz = 0.d0
    fr = 0.d0

    do i=1,n_vz
        do j=1,n_vr
            u2 = (zb - vz(i))**2.0 + vr(j)**2.0
            erel = u2_to_erel*u2
            sig = fn(erel)
            fr(j) = sig*sqrt(u2)*exp(-(vz(i)**2.0 + vr(j)**2.0))*vr(j)
        enddo
        fz(i) = simpsons_rule(fr, dvr)
    enddo

    sig_eff = (2.0/sqrt(PI))*simpsons_rule(fz, dvz)
    rate = sig_eff*v_therm

end subroutine bt_maxwellian_eb