Get beam-beam neutralization/cx rates
Type | Intent | Optional | Attributes | Name | ||
---|---|---|---|---|---|---|
real(kind=Float64), | intent(in), | dimension(nlevs) | :: | denn | Neutral density [cm^-3] |
|
real(kind=Float64), | intent(in), | dimension(3) | :: | vi | Ion velocity [cm/s] |
|
real(kind=Float64), | intent(in), | dimension(3) | :: | vn | Neutral velocity [cm/s] |
|
real(kind=Float64), | intent(out), | dimension(nlevs) | :: | rates | Reaction rates [1/s] |
subroutine bb_cx_rates(denn, vi, vn, rates)
!+ Get beam-beam neutralization/cx rates
real(Float64), dimension(nlevs), intent(in) :: denn
!+ Neutral density [cm^-3]
real(Float64), dimension(3), intent(in) :: vi
!+ Ion velocity [cm/s]
real(Float64), dimension(3), intent(in) :: vn
!+ Neutral velocity [cm/s]
real(Float64), dimension(nlevs), intent(out) :: rates
!+ Reaction rates [1/s]
real(Float64), dimension(nlevs,nlevs) :: neut !!rate coeff
real(Float64) :: eb !! relative Energy
type(InterpolCoeffs1D) :: c
real(Float64) :: dlogE, logEmin, logeb
real(Float64) :: vrel !! relative velocity
integer :: ebi, neb, err
!Eeff
vrel=norm2(vi-vn)
eb=v2_to_E_per_amu*vrel**2 ! [kev/amu]
logeb = log10(eb)
logEmin = tables%H_H_cx_cross%logemin
dlogE = tables%H_H_cx_cross%dlogE
neb = tables%H_H_cx_cross%nenergy
call interpol_coeff(logEmin,dlogE,neb,logeb,c,err)
ebi = c%i
if(err.eq.1) then
if(inputs%verbose.ge.0) then
write(*,'(a)') "BB_CX_RATES: Eb out of range of H_H_cx table. Using nearest energy value."
write(*,'("eb = ",ES10.3," [keV]")') eb
endif
if(ebi.lt.1) then
ebi=1
c%b1=1.0 ; c%b2=0.0
else
ebi=neb-1
c%b1=0.0 ; c%b2=1.0
endif
endif
neut(:,:) = (c%b1*tables%H_H_cx_cross%log_cross(:,:,ebi) + &
c%b2*tables%H_H_cx_cross%log_cross(:,:,ebi+1))
where (neut.lt.tables%H_H_cx_cross%minlog_cross)
neut = 0.d0
elsewhere
neut = exp(neut*log_10)
end where
rates=matmul(neut,denn)*vrel
end subroutine bb_cx_rates