bt_cx_rates Subroutine

  1. fidasim.f90
  2. libfida
  3. bt_cx_rates

public subroutine bt_cx_rates(plasma, denn, an, vi, rates)

Get beam-target neutralization/cx rates

Arguments

TypeIntentOptionalAttributesName
type(LocalProfiles), intent(in) :: plasma

Plasma parameters (for neutral temperature and vrot)
real(kind=Float64), intent(in), dimension(nlevs):: denn

Neutral density [cm^-3]
real(kind=Float64), intent(in) :: an

Neutral mass [amu]
real(kind=Float64), intent(in), dimension(3):: vi

Ion velocity [cm/s]
real(kind=Float64), intent(out), dimension(nlevs):: rates

Reaction rates [1/s]

Calls

proc~~bt_cx_rates~~CallsGraph proc~bt_cx_rates bt_cx_rates interface~interpol_coeff interpol_coeff proc~bt_cx_rates->interface~interpol_coeff

Called by

proc~~bt_cx_rates~~CalledByGraph proc~bt_cx_rates bt_cx_rates proc~pnpa_mc pnpa_mc proc~pnpa_mc->proc~bt_cx_rates proc~pnpa_f pnpa_f proc~pnpa_f->proc~bt_cx_rates proc~pfida_f pfida_f proc~pfida_f->proc~bt_cx_rates proc~pfida_mc pfida_mc proc~pfida_mc->proc~bt_cx_rates proc~fida_weights_los fida_weights_los proc~fida_weights_los->proc~bt_cx_rates program~fidasim fidasim program~fidasim->proc~pnpa_mc program~fidasim->proc~pnpa_f program~fidasim->proc~pfida_f program~fidasim->proc~pfida_mc program~fidasim->proc~fida_weights_los

Contents

Source Code

bt_cx_rates

Source Code

subroutine bt_cx_rates(plasma, denn, an, vi, rates)
    !+ Get beam-target neutralization/cx rates
    type(LocalProfiles), intent(in)              :: plasma
        !+ Plasma parameters (for neutral temperature and vrot)
    real(Float64), dimension(nlevs), intent(in)  :: denn
        !+ Neutral density [cm^-3]
    real(Float64), intent(in)                    :: an
        !+ Neutral mass [amu]
    real(Float64), dimension(3),     intent(in)  :: vi
        !+ Ion velocity [cm/s]
    real(Float64), dimension(nlevs), intent(out) :: rates
        !+ Reaction rates [1/s]

    real(Float64) :: logEmin, dlogE, logeb_amu, eb_amu
    real(Float64) :: logTmin, dlogT, logti_amu, vrel
    integer :: neb, nt
    type(InterpolCoeffs2D) :: c
    real(Float64) :: b11, b12, b21, b22
    real(Float64), dimension(nlevs,nlevs) :: H_H_rate
    integer :: ebi, tii, n, err_status

    H_H_rate = 0.d0

    vrel=norm2(vi-plasma%vrot)
    eb_amu=v2_to_E_per_amu*vrel**2  ! [kev/amu]

    logeb_amu = log10(eb_amu)
    logti_amu = log10(plasma%ti/an)

    !!H_H
    err_status = 1
    logEmin = tables%H_H_cx_rate%logemin
    logTmin = tables%H_H_cx_rate%logtmin
    dlogE = tables%H_H_cx_rate%dlogE
    dlogT = tables%H_H_cx_rate%dlogT
    neb = tables%H_H_cx_rate%nenergy
    nt = tables%H_H_cx_rate%ntemp
    call interpol_coeff(logEmin, dlogE, neb, logTmin, dlogT, nt, &
                        logeb_amu, logti_amu, c, err_status)
    ebi = c%i
    tii = c%j
    b11 = c%b11
    b12 = c%b12
    b21 = c%b21
    b22 = c%b22
    if(err_status.eq.1) then
        if(inputs%verbose.ge.0) then
            write(*,'(a)') "BT_CX_RATES: Eb or Ti out of range of H_H_CX table. Setting H_H_CX rates to zero"
            write(*,'("eb/amu = ",ES10.3," [keV/amu]")') eb_amu
            write(*,'("ti/amu = ",ES10.3," [keV/amu]")') plasma%ti/an
        endif
        rates = 0.0
        return
    endif

    H_H_rate = (b11*tables%H_H_cx_rate%log_rate(:,:,ebi,tii)   + &
                b12*tables%H_H_cx_rate%log_rate(:,:,ebi,tii+1) + &
                b21*tables%H_H_cx_rate%log_rate(:,:,ebi+1,tii) + &
                b22*tables%H_H_cx_rate%log_rate(:,:,ebi+1,tii+1))

    where (H_H_rate.lt.tables%H_H_cx_rate%minlog_rate)
        H_H_rate = 0.d0
    elsewhere
        H_H_rate = exp(H_H_rate*log_10) !cm^3/s
    end where

    rates=matmul(H_H_rate,denn) !1/s

end subroutine bt_cx_rates