colrad Subroutine

  1. fidasim.f90
  2. libfida
  3. colrad

public subroutine colrad(plasma, i_type, vn, dt, states, dens, photons)

Evolve density of states in time dt via collisional radiative model

Arguments

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

Plasma parameters
integer, intent(in) :: i_type

Ion/Neutral type (beam,thermal)
real(kind=Float64), intent(in), dimension(:):: vn

Neutral velocitiy [cm/s]
real(kind=Float64), intent(in) :: dt

Time interval [s]
real(kind=Float64), intent(inout), dimension(:):: states

Density of states
real(kind=Float64), intent(out), dimension(nlevs):: dens

Density of neutrals
real(kind=Float64), intent(out) :: photons

Emitted photons(3->2)

Calls

proc~~colrad~~CallsGraph proc~colrad colrad proc~eigen eigen proc~colrad->proc~eigen proc~get_rate_matrix get_rate_matrix proc~colrad->proc~get_rate_matrix proc~matinv matinv proc~colrad->proc~matinv proc~hqr2 hqr2 proc~eigen->proc~hqr2 proc~balback balback proc~eigen->proc~balback proc~elmtrans elmtrans proc~eigen->proc~elmtrans proc~elmhes elmhes proc~eigen->proc~elmhes proc~balance balance proc~eigen->proc~balance interface~interpol_coeff interpol_coeff proc~get_rate_matrix->interface~interpol_coeff proc~ludcmp ludcmp proc~matinv->proc~ludcmp 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~interpol2d_coeff interpol2D_coeff interface~interpol_coeff->proc~interpol2d_coeff proc~interpol1d_coeff_arr interpol1D_coeff_arr interface~interpol_coeff->proc~interpol1d_coeff_arr proc~hqrvec hqrvec proc~hqr2->proc~hqrvec proc~outerprod outerprod proc~ludcmp->proc~outerprod proc~interpol2d_coeff_arr->proc~interpol2d_coeff proc~interpol1d_coeff_arr->proc~interpol1d_coeff

Called by

proc~~colrad~~CalledByGraph proc~colrad colrad proc~fida_weights_los fida_weights_los proc~fida_weights_los->proc~colrad proc~attenuate attenuate proc~attenuate->proc~colrad proc~npa_weights npa_weights proc~npa_weights->proc~attenuate proc~npa_f npa_f proc~npa_f->proc~attenuate 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~attenuate

Contents

Source Code

colrad

Source Code

subroutine colrad(plasma,i_type,vn,dt,states,dens,photons)
    !+ Evolve density of states in time `dt` via collisional radiative model
    type(LocalProfiles), intent(in)              :: plasma
        !+ Plasma parameters
    integer, intent(in)                          :: i_type
        !+ Ion/Neutral type (beam,thermal)
    real(Float64), dimension(:), intent(in)      :: vn
        !+ Neutral velocitiy [cm/s]
    real(Float64), intent(in)                    :: dt
        !+ Time interval [s]
    real(Float64), dimension(:), intent(inout)   :: states
        !+ Density of states
    real(Float64), dimension(nlevs), intent(out) :: dens
        !+ Density of neutrals
    real(Float64), intent(out)                   :: photons
        !+ Emitted photons(3->2)

    real(Float64), dimension(nlevs,nlevs) :: matrix  !! Matrix
    real(Float64) :: b_amu
    real(Float64) :: vnet_square    !! net velocity of neutrals squared
    real(Float64) :: eb             !! Energy of the fast neutral

    real(Float64), dimension(nlevs,nlevs) :: eigvec, eigvec_inv
    real(Float64), dimension(nlevs) :: eigval, coef
    real(Float64), dimension(nlevs) :: exp_eigval_dt
    real(Float64) :: iflux !!Initial total flux
    integer :: n

    photons=0.d0
    dens=0.d0

    iflux=sum(states)
    if(iflux.lt.colrad_threshold .and. inputs%calc_npa.eq.0)then
        return
    endif

    if(.not.plasma%in_plasma) then
        dens = states*dt
        return
    endif

    if(i_type.eq.beam_ion) then
        b_amu = inputs%ab
    else
        b_amu = inputs%ai
    endif
    vnet_square=dot_product(vn-plasma%vrot,vn-plasma%vrot)  ![cm/s]
    eb = v2_to_E_per_amu*b_amu*vnet_square ![kev]
    call get_rate_matrix(plasma, i_type, eb, matrix)

    call eigen(nlevs,matrix, eigvec, eigval)
    call matinv(eigvec, eigvec_inv)
    coef = matmul(eigvec_inv, states)!coeffs determined from states at t=0
    exp_eigval_dt = exp(eigval*dt)   ! to improve speed (used twice)
    do n=1,nlevs
        if(eigval(n).eq.0.0) eigval(n)=eigval(n)+1 !protect against dividing by zero
    enddo

    states = matmul(eigvec, coef * exp_eigval_dt)  ![neutrals/cm^3/s]!
    dens   = matmul(eigvec,coef*(exp_eigval_dt-1.d0)/eigval)

    if ((minval(states).lt.0).or.(minval(dens).lt.0)) then
        do n=1,nlevs
            if(states(n).lt.0) states(n)=0.d0
            if(dens(n).lt.0) dens(n)=0.d0
        enddo
    endif

    photons=dens(3)*tables%einstein(2,3) !! - [Ph/(s*cm^3)] - !!

end subroutine colrad