bt_maxwellian_n Subroutine

  1. atomic_tables.f90
  2. atomic_tables
  3. bt_maxwellian_n

public subroutine bt_maxwellian_n(fn, T, eb, am, ab, n, rate)

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

Arguments

TypeIntentOptionalAttributesName
public function fn(a, b)

Cross section function

Arguments
TypeIntentOptionalAttributesName
real(kind=8), intent(in) :: a
integer, intent(in) :: b
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]
integer, intent(in) :: n

Initial atomic energy level/state
real(kind=Float64), intent(out) :: rate

Reaction Rate []

Calls

proc~~bt_maxwellian_n~~CallsGraph proc~bt_maxwellian_n bt_maxwellian_n proc~simpsons_rule simpsons_rule proc~bt_maxwellian_n->proc~simpsons_rule

Called by

proc~~bt_maxwellian_n~~CalledByGraph proc~bt_maxwellian_n bt_maxwellian_n interface~bt_maxwellian bt_maxwellian interface~bt_maxwellian->proc~bt_maxwellian_n proc~write_bt_d_he3 write_bt_D_He3 proc~write_bt_d_he3->interface~bt_maxwellian proc~write_bt_h_aq write_bt_H_Aq proc~write_bt_h_aq->interface~bt_maxwellian proc~write_bt_d_t write_bt_D_T proc~write_bt_d_t->interface~bt_maxwellian proc~write_bt_h_h write_bt_H_H proc~write_bt_h_h->interface~bt_maxwellian proc~write_bt_d_d write_bt_D_D proc~write_bt_d_d->interface~bt_maxwellian proc~write_bt_h_e write_bt_H_e proc~write_bt_h_e->interface~bt_maxwellian program~generate_tables generate_tables program~generate_tables->proc~write_bt_d_he3 program~generate_tables->proc~write_bt_h_aq program~generate_tables->proc~write_bt_h_h program~generate_tables->proc~write_bt_d_d program~generate_tables->proc~write_bt_h_e

Contents

Source Code

bt_maxwellian_n

Source Code

subroutine bt_maxwellian_n(fn, T, eb, am, ab, n, rate)
    !+ Calculates Maxwellian reaction rate for a beam with atomic mass `ab`, energy `eb`, and energy level `n`
    !+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, b)
            !+Cross section function
            real(8)              :: fn !sigma
            real(8), intent(in)  :: a !eb
            integer, intent(in)  :: b !n
        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]
    integer, intent(in)        :: n
        !+Initial atomic energy level/state
    real(Float64), intent(out) :: rate
        !+Reaction Rate [\(cm^3/s\)]

    logical :: dxc
    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)
    dE = (13.6d-3)/(n**2.0)

    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
    if(ared.lt.0.5) ared = 1.0 !for electron interactions

    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
            if(erel.ge.dE) then
                sig = fn(erel/ared,n)
            else
                sig = 0.d0
            endif
            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_n