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
Type | Intent | Optional | Attributes | Name | ||||||||||||||||
---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
public function fn(a)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] |
||||||||||||||||
real(kind=Float64), | intent(out) | :: | rate | Reaction Rate [] |
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