Gets rate matrix for use in colrad
| Type | Intent | Optional | Attributes | Name | ||
|---|---|---|---|---|---|---|
| type(LocalProfiles), | intent(in) | :: | plasma | Plasma parameters |
||
| real(kind=Float64), | intent(in) | :: | ab | "Beam" ion mass [amu] |
||
| real(kind=Float64), | intent(in) | :: | eb | "Beam" ion energy [keV] |
||
| real(kind=Float64), | intent(out), | dimension(nlevs,nlevs) | :: | rmat | Rate matrix |
subroutine get_rate_matrix(plasma, ab, eb, rmat)
!+ Gets rate matrix for use in [[libfida:colrad]]
type(LocalProfiles), intent(in) :: plasma
!+ Plasma parameters
real(Float64), intent(in) :: ab
!+ "Beam" ion mass [amu]
real(Float64), intent(in) :: eb
!+ "Beam" ion energy [keV]
real(Float64), dimension(nlevs,nlevs), intent(out) :: rmat
!+ Rate matrix
real(Float64) :: ai
real(Float64) :: logEmin, dlogE, logeb, logeb_amu
real(Float64) :: logTmin, dlogT, logti, logti_amu, logte
integer :: neb, nt, i
type(InterpolCoeffs2D) :: c
real(Float64) :: b11, b12, b21, b22, dene, deni(max_species), denimp
real(Float64), dimension(nlevs,nlevs) :: H_H_pop_i, H_H_pop, H_e_pop, H_Aq_pop
real(Float64), dimension(nlevs) :: H_H_depop_i, H_H_depop, H_e_depop, H_Aq_depop
integer :: ebi, tii, tei, n, err_status
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
deni = plasma%deni
dene = plasma%dene
denimp = plasma%denimp
logeb_amu = log10(eb/ab)
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
do i=1, n_thermal
H_H_pop_i = 0.d0
H_H_depop_i = 0.d0
logti_amu = log10(plasma%ti/thermal_mass(i))
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)') "GET_RATE_MATRIX: Eb or Ti out of range of H_H table. Setting H_H rates to zero"
write(*,'("eb/amu = ",ES10.3," [keV/amu]")') eb/ab
write(*,'("ti/amu = ",ES10.3," [keV/amu]")') plasma%ti/thermal_mass(i)
endif
else
H_H_pop_i = (b11*tables%H_H%log_pop(:,:,ebi,tii) + &
b12*tables%H_H%log_pop(:,:,ebi,tii+1) + &
b21*tables%H_H%log_pop(:,:,ebi+1,tii) + &
b22*tables%H_H%log_pop(:,:,ebi+1,tii+1))
where (H_H_pop_i.lt.tables%H_H%minlog_pop)
H_H_pop_i = 0.d0
elsewhere
H_H_pop_i = deni(i) * exp(H_H_pop_i*log_10)
end where
H_H_depop_i = (b11*tables%H_H%log_depop(:,ebi,tii) + &
b12*tables%H_H%log_depop(:,ebi,tii+1) + &
b21*tables%H_H%log_depop(:,ebi+1,tii) + &
b22*tables%H_H%log_depop(:,ebi+1,tii+1))
where (H_H_depop_i.lt.tables%H_H%minlog_depop)
H_H_depop_i = 0.d0
elsewhere
H_H_depop_i = deni(i) * exp(H_H_depop_i*log_10)
end where
endif
H_H_pop = H_H_pop + H_H_pop_i
H_H_depop = H_H_depop + H_H_depop_i
enddo
!!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_amu, 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/amu = ",ES10.3," [keV/amu]")') eb/ab
write(*,'("te = ",ES10.3," [keV]")') plasma%te
endif
else
H_e_pop = (b11*tables%H_e%log_pop(:,:,ebi,tei) + &
b12*tables%H_e%log_pop(:,:,ebi,tei+1) + &
b21*tables%H_e%log_pop(:,:,ebi+1,tei) + &
b22*tables%H_e%log_pop(:,:,ebi+1,tei+1))
where (H_e_pop.lt.tables%H_e%minlog_pop)
H_e_pop = 0.d0
elsewhere
H_e_pop = dene * exp(H_e_pop*log_10)
end where
H_e_depop = (b11*tables%H_e%log_depop(:,ebi,tei) + &
b12*tables%H_e%log_depop(:,ebi,tei+1) + &
b21*tables%H_e%log_depop(:,ebi+1,tei) + &
b22*tables%H_e%log_depop(:,ebi+1,tei+1))
where (H_e_depop.lt.tables%H_e%minlog_depop)
H_e_depop = 0.d0
elsewhere
H_e_depop = dene * exp(H_e_depop*log_10)
end where
endif
!!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_amu, 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/amu = ",ES10.3," [keV/amu]")') eb
write(*,'("ti = ",ES10.3," [keV]")') plasma%ti
endif
else
H_Aq_pop = (b11*tables%H_Aq%log_pop(:,:,ebi,tii) + &
b12*tables%H_Aq%log_pop(:,:,ebi,tii+1) + &
b21*tables%H_Aq%log_pop(:,:,ebi+1,tii) + &
b22*tables%H_Aq%log_pop(:,:,ebi+1,tii+1))
where (H_Aq_pop.lt.tables%H_Aq%minlog_pop)
H_Aq_pop = 0.d0
elsewhere
H_Aq_pop = denimp * exp(H_Aq_pop*log_10)
end where
H_Aq_depop = (b11*tables%H_Aq%log_depop(:,ebi,tii) + &
b12*tables%H_Aq%log_depop(:,ebi,tii+1) + &
b21*tables%H_Aq%log_depop(:,ebi+1,tii) + &
b22*tables%H_Aq%log_depop(:,ebi+1,tii+1))
where (H_Aq_depop.lt.tables%H_Aq%minlog_depop)
H_Aq_depop = 0.d0
elsewhere
H_Aq_depop = denimp * exp(H_Aq_depop*log_10)
end where
endif
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