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
and temperature T
which has a cross section given by the function fnm
Type | Intent | Optional | Attributes | Name | ||||||||||||||||||||||||||||||
---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
public function fnm(a, b, c)Cross section function Arguments
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 |
||||||||||||||||||||||||||||||
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 |
subroutine bt_maxwellian_n_m(fnm, 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` and temperature `T` which has a cross section given by the function `fnm`
interface
function fnm(a, b, c)
!+Cross section function
real(8) :: fnm !sigma
real(8), intent(in) :: a !eb
integer, intent(in) :: b !n
integer, intent(in) :: c !m
end function fnm
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
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 :: n_vr
real(Float64) :: vr_max, dvr
real(Float64), dimension(32) :: vr
real(Float64), dimension(32) :: fr
integer :: n_vz
real(Float64) :: vz_max, dvz
real(Float64), dimension(62) :: vz
real(Float64), dimension(62) :: 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
n_vr = 32
vr_max = 4.d0
dvr = vr_max/(n_vr - 1.d0)
do i=1,n_vr
vr(i) = (i-1)*dvr
enddo
n_vz = 62
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 = fnm(erel/ared, 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_n_m