bt_maxwellian_q_n_m Subroutine

  1. atomic_tables.f90
  2. atomic_tables
  3. bt_maxwellian_q_n_m

public subroutine bt_maxwellian_q_n_m(fqnm, q, T, eb, am, ab, n, m, rate, deexcit)

Calculates Maxwellian reaction rate for a n(\rightarrow)m transition due to a beam with atomic mass ab and energy eb firing into a target with atomic mass am, temperature T, and charge q which has a cross section given by the function fqnm

Arguments

TypeIntentOptionalAttributesName
public function fqnm(a, b, c, d)

Cross section function

Arguments
TypeIntentOptionalAttributesName
real(kind=8), intent(in) :: a
integer, intent(in) :: b
integer, intent(in) :: c
integer, intent(in) :: d
Return Value real(kind=8)
integer, intent(in) :: q

Target charge
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
integer, intent(in) :: m

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

Reaction Rate []
logical, intent(in), optional :: deexcit

Calculate de-excitation reaction rate

Calls

proc~~bt_maxwellian_q_n_m~~CallsGraph proc~bt_maxwellian_q_n_m bt_maxwellian_q_n_m proc~simpsons_rule simpsons_rule proc~bt_maxwellian_q_n_m->proc~simpsons_rule

Called by

proc~~bt_maxwellian_q_n_m~~CalledByGraph proc~bt_maxwellian_q_n_m bt_maxwellian_q_n_m interface~bt_maxwellian bt_maxwellian interface~bt_maxwellian->proc~bt_maxwellian_q_n_m 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_q_n_m

Source Code

subroutine bt_maxwellian_q_n_m(fqnm, q, T, eb, am, ab, n, m, rate, deexcit)
    !+ Calculates Maxwellian reaction rate for a `n`\(\rightarrow)`m` transition due to a beam with atomic mass `ab` and energy `eb`
    !+firing into a target with atomic mass `am`, temperature `T`, and charge `q` which has a cross section given by the function `fqnm`
    interface
        function fqnm(a, b, c, d)
            !+Cross section function
            real(8)             :: fqnm !sigma
            real(8), intent(in) :: a !eb
            integer, intent(in) :: b !q
            integer, intent(in) :: c !n
            integer, intent(in) :: d !m
        end function fqnm
    end interface
    integer, intent(in)           :: q
        !+Target charge
    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
    integer, intent(in)           :: m
        !+Final atomic energy level/state
    real(Float64), intent(out)    :: rate
        !+Reaction Rate [\(cm^3/s\)]
    logical, intent(in), optional :: deexcit
        !+Calculate de-excitation reaction rate

    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, dE, factor, En, Em, v_therm

    integer :: i, j

    if(present(deexcit)) then
        dxc = deexcit
    else
        dxc = .False.
    endif

    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

    En = (13.6d-3)*(1.0 - (1.d0/n)**2.0)
    Em = (13.6d-3)*(1.0 - (1.d0/m)**2.0)
    dE = Em - En

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

    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(dxc) then
                factor = (erel + dE)/erel
                erel = erel + dE
            else
                factor = 1.0
            endif
            if(erel.ge.dE) then
                sig = fqnm(erel/ared, q, n, m)
            else
                sig = 0.d0
            endif
            fr(j) = factor*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
    if(dxc) rate = rate*(real(n)/real(m))**2.0

end subroutine bt_maxwellian_q_n_m