Gets rate matrix for use in colrad
| Type | Intent | Optional | Attributes | Name | ||
|---|---|---|---|---|---|---|
| type(LocalProfiles), | intent(in) | :: | plasma | Plasma parameters |
||
| integer, | intent(in) | :: | i_type | Ion type |
||
| real(kind=Float64), | intent(in) | :: | eb | Ion energy [keV] |
||
| real(kind=Float64), | intent(out), | dimension(nlevs,nlevs) | :: | rmat | Rate matrix |
subroutine get_rate_matrix(plasma, i_type, eb, rmat)
!+ Gets rate matrix for use in [[libfida:colrad]]
type(LocalProfiles), intent(in) :: plasma
!+ Plasma parameters
integer, intent(in) :: i_type
!+ Ion type
real(Float64), intent(in) :: eb
!+ Ion energy [keV]
real(Float64), dimension(nlevs,nlevs), intent(out) :: rmat
!+ Rate matrix
real(Float64) :: logEmin, dlogE, logeb
real(Float64) :: logTmin, dlogT, logti, logte
integer :: neb, nt
type(InterpolCoeffs2D) :: c
real(Float64) :: b11, b12, b21, b22, dene, denp, denimp
real(Float64), dimension(nlevs,nlevs) :: H_H_pop, H_e_pop, H_Aq_pop
real(Float64), dimension(nlevs) :: H_H_depop, H_e_depop, H_Aq_depop
integer :: ebi, tii, tei, n, err_status
real(Float64) :: mullog10 = log(10.0d0)
H_H_pop = 0.d0
H_e_pop = 0.d0
H_Aq_pop = 0.d0
H_H_depop = 0.d0
H_e_depop = 0.d0
H_Aq_depop = 0.d0
denp = plasma%denp
dene = plasma%dene
denimp = plasma%denimp
logeb = log10(eb)
logti = log10(plasma%ti)
logte = log10(plasma%te)
!!H_H
err_status = 1
logEmin = tables%H_H%logemin
logTmin = tables%H_H%logtmin
dlogE = tables%H_H%dlogE
dlogT = tables%H_H%dlogT
neb = tables%H_H%nenergy
nt = tables%H_H%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)') "GET_RATE_MATRIX: Eb or Ti out of range of H_H table. Setting H_H rates to zero"
write(*,'("eb = ",ES10.3," [keV]")') eb
write(*,'("ti = ",ES10.3," [keV]")') plasma%ti
endif
denp = 0.d0
endif
H_H_pop = (b11*tables%H_H%log_pop(:,:,ebi,tii,i_type) + &
b12*tables%H_H%log_pop(:,:,ebi,tii+1,i_type) + &
b21*tables%H_H%log_pop(:,:,ebi+1,tii,i_type) + &
b22*tables%H_H%log_pop(:,:,ebi+1,tii+1,i_type))
where (H_H_pop.lt.tables%H_H%minlog_pop)
H_H_pop = 0.d0
elsewhere
H_H_pop = denp * exp(H_H_pop*mullog10)
end where
H_H_depop = (b11*tables%H_H%log_depop(:,ebi,tii,i_type) + &
b12*tables%H_H%log_depop(:,ebi,tii+1,i_type) + &
b21*tables%H_H%log_depop(:,ebi+1,tii,i_type) + &
b22*tables%H_H%log_depop(:,ebi+1,tii+1,i_type))
where (H_H_depop.lt.tables%H_H%minlog_depop)
H_H_depop = 0.d0
elsewhere
H_H_depop = denp * exp(H_H_depop*mullog10)
end where
!!H_e
err_status = 1
logEmin = tables%H_e%logemin
logTmin = tables%H_e%logtmin
dlogE = tables%H_e%dlogE
dlogT = tables%H_e%dlogT
neb = tables%H_e%nenergy
nt = tables%H_e%ntemp
call interpol_coeff(logEmin, dlogE, neb, logTmin, dlogT, nt, &
logeb, logte, c, err_status)
ebi = c%i
tei = 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)') "GET_RATE_MATRIX: Eb or Te out of range of H_e table. Setting H_e rates to zero"
write(*,'("eb = ",ES10.3," [keV]")') eb
write(*,'("te = ",ES10.3," [keV]")') plasma%te
endif
dene = 0.d0
endif
H_e_pop = (b11*tables%H_e%log_pop(:,:,ebi,tei,i_type) + &
b12*tables%H_e%log_pop(:,:,ebi,tei+1,i_type) + &
b21*tables%H_e%log_pop(:,:,ebi+1,tei,i_type) + &
b22*tables%H_e%log_pop(:,:,ebi+1,tei+1,i_type))
where (H_e_pop.lt.tables%H_e%minlog_pop)
H_e_pop = 0.d0
elsewhere
H_e_pop = dene * exp(H_e_pop*mullog10)
end where
H_e_depop = (b11*tables%H_e%log_depop(:,ebi,tei,i_type) + &
b12*tables%H_e%log_depop(:,ebi,tei+1,i_type) + &
b21*tables%H_e%log_depop(:,ebi+1,tei,i_type) + &
b22*tables%H_e%log_depop(:,ebi+1,tei+1,i_type))
where (H_e_depop.lt.tables%H_e%minlog_depop)
H_e_depop = 0.d0
elsewhere
H_e_depop = dene * exp(H_e_depop*mullog10)
end where
!!H_Aq
err_status = 1
logEmin = tables%H_Aq%logemin
logTmin = tables%H_Aq%logtmin
dlogE = tables%H_Aq%dlogE
dlogT = tables%H_Aq%dlogT
neb = tables%H_Aq%nenergy
nt = tables%H_Aq%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)') "GET_RATE_MATRIX: Eb or Ti out of range of H_Aq table. Setting H_Aq rates to zero"
write(*,'("eb = ",ES10.3," [keV]")') eb
write(*,'("ti = ",ES10.3," [keV]")') plasma%ti
endif
denimp = 0.d0
endif
H_Aq_pop = (b11*tables%H_Aq%log_pop(:,:,ebi,tii,i_type) + &
b12*tables%H_Aq%log_pop(:,:,ebi,tii+1,i_type) + &
b21*tables%H_Aq%log_pop(:,:,ebi+1,tii,i_type) + &
b22*tables%H_Aq%log_pop(:,:,ebi+1,tii+1,i_type))
where (H_Aq_pop.lt.tables%H_Aq%minlog_pop)
H_Aq_pop = 0.d0
elsewhere
H_Aq_pop = denimp * exp(H_Aq_pop*mullog10)
end where
H_Aq_depop = (b11*tables%H_Aq%log_depop(:,ebi,tii,i_type) + &
b12*tables%H_Aq%log_depop(:,ebi,tii+1,i_type) + &
b21*tables%H_Aq%log_depop(:,ebi+1,tii,i_type) + &
b22*tables%H_Aq%log_depop(:,ebi+1,tii+1,i_type))
where (H_Aq_depop.lt.tables%H_Aq%minlog_depop)
H_Aq_depop = 0.d0
elsewhere
H_Aq_depop = denimp * exp(H_Aq_depop*mullog10)
end where
rmat = tables%einstein + H_H_pop + H_e_pop + H_Aq_pop
do n=1,nlevs
rmat(n,n) = -sum(tables%einstein(:,n)) - H_H_depop(n) - H_e_depop(n) - H_Aq_depop(n)
enddo
end subroutine get_rate_matrix