bt_cx_rates Subroutine

  1. fidasim.f90
  2. libfida
  3. bt_cx_rates

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

Get beam-target neutralization/cx rates

Arguments

Type IntentOptional AttributesName
type(LocalProfiles), intent(in) :: plasma

Plasma parameters
real(kind=Float64), intent(in), dimension(nlevs):: denn

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

Ion velocity [cm/s]
integer, intent(in) :: i_type

Ion type
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 proc~interpol2d_coeff interpol2D_coeff interface~interpol_coeff->proc~interpol2d_coeff proc~interpol1d_coeff interpol1D_coeff interface~interpol_coeff->proc~interpol1d_coeff proc~interpol2d_coeff_arr interpol2D_coeff_arr interface~interpol_coeff->proc~interpol2d_coeff_arr proc~interpol1d_coeff_arr interpol1D_coeff_arr interface~interpol_coeff->proc~interpol1d_coeff_arr proc~interpol2d_coeff_arr->proc~interpol2d_coeff proc~interpol1d_coeff_arr->proc~interpol1d_coeff

Called by

proc~~bt_cx_rates~~CalledByGraph proc~bt_cx_rates bt_cx_rates proc~fida_weights_los fida_weights_los proc~fida_weights_los->proc~bt_cx_rates proc~get_beam_cx_rate get_beam_cx_rate proc~get_beam_cx_rate->proc~bt_cx_rates proc~npa_weights npa_weights proc~npa_weights->proc~get_beam_cx_rate proc~npa_f npa_f proc~npa_f->proc~get_beam_cx_rate program~fidasim fidasim program~fidasim->proc~fida_weights_los program~fidasim->proc~npa_weights program~fidasim->proc~npa_f proc~npa_mc npa_mc program~fidasim->proc~npa_mc proc~npa_mc->proc~get_beam_cx_rate

Contents

Source Code

bt_cx_rates

Source Code

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

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

    real(Float64) :: mullog10 = log(10.0d0)

    H_H_rate = 0.d0

    if(i_type.eq.beam_ion) then
        b_amu = inputs%ab
    else
        b_amu = inputs%ai
    endif
    vrel=norm2(vi-plasma%vrot)
    eb=b_amu*v2_to_E_per_amu*vrel**2  ! [kev/amu]

    logeb = log10(eb)
    logti = log10(plasma%ti)

    !!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, logti, 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 = ",ES10.3," [keV]")') eb
            write(*,'("ti = ",ES10.3," [keV]")') plasma%ti
        endif
        rates = 0.0
        return
    endif

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

    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*mullog10) !cm^3/s
    end where

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

end subroutine bt_cx_rates