src/gray_core.f90: fast absorption computation
This implements an explicit formula for the absorption coefficient in the low temperature limit, available with iwarm=4.
This commit is contained in:
Michele Guerini Rocco
parent
899f524782
commit
ffed0dc1c5
GPG Key ID:
BFBAF4C975F76450
400
src/absorption.f90
Normal file
400
src/absorption.f90
Normal file
@ -0,0 +1,400 @@
|
|||||||
|
! This modules implements a simplified formula for the plasma absorption
|
||||||
|
! coefficient close to the harmonics valid for low temperature (< 10keV)
|
||||||
|
! and low harmonic numbers (<6).
|
||||||
|
!
|
||||||
|
! Reference: 1983 Nuclear Fusion 23 1153 (https://doi.org/bm6zcr)
|
||||||
|
!
|
||||||
|
module absorption
|
||||||
|
|
||||||
|
use const_and_precisions, only : wp_
|
||||||
|
|
||||||
|
implicit none
|
||||||
|
|
||||||
|
private
|
||||||
|
public :: alpha_fast
|
||||||
|
|
||||||
|
contains
|
||||||
|
|
||||||
|
pure function line_shape(Y, Npl2, theta) result(phi)
|
||||||
|
! The line shape function of the first harmonic,
|
||||||
|
! with exception of the X mode in the tenuous plasma regime:
|
||||||
|
!
|
||||||
|
! φ(ζ) = Im(-1/Z(ζ))
|
||||||
|
!
|
||||||
|
! From the same considerations as in `line_shape_harmonics`, we have:
|
||||||
|
!
|
||||||
|
! φ(ζ) = 1/√π exp(ζ²)/[1 + erfi(ζ)²]
|
||||||
|
!
|
||||||
|
! Note: Since we already have an implementation of Z(ζ) we
|
||||||
|
! use that and write Φ(ζ) = Im(Z)/[Re(Z)² + Im(Z)²].
|
||||||
|
!
|
||||||
|
use const_and_precisions, only: zero
|
||||||
|
use dispersion, only: zetac
|
||||||
|
|
||||||
|
implicit none
|
||||||
|
|
||||||
|
! function arguments
|
||||||
|
|
||||||
|
! CMA Y variable: Y=ω_c/ω
|
||||||
|
real(wp_), intent(in) :: Y
|
||||||
|
! squared parallel refractive index: N∥²=N²cos²θ
|
||||||
|
real(wp_), intent(in) :: Npl2
|
||||||
|
! adimensional temperature: Θ = kT/mc²
|
||||||
|
real(wp_), intent(in) :: theta
|
||||||
|
! line shape Φ(ζ)
|
||||||
|
real(wp_) :: phi
|
||||||
|
|
||||||
|
! local variables
|
||||||
|
real(wp_) :: a, b, z
|
||||||
|
integer :: err
|
||||||
|
|
||||||
|
z = (1 - Y)/sqrt(2 * theta * Npl2) ! ζ
|
||||||
|
|
||||||
|
call zetac(z, zero, a, b, err) ! a=Re(Z(ζ)), b=Im(Z(ζ))
|
||||||
|
phi = b/(a**2 + b**2) ! Im(-1/Z)
|
||||||
|
end function line_shape
|
||||||
|
|
||||||
|
|
||||||
|
pure function line_shape_harmonics(Y, Npl2, theta, n) result(phi)
|
||||||
|
! The line shape function of all modes at higher harmonics (n>1)
|
||||||
|
! and the X mode at the first harmonic:
|
||||||
|
!
|
||||||
|
! φ(ζ) = Im(Z(ζ))
|
||||||
|
!
|
||||||
|
! where Z(ζ) = 1/√π ∫dz exp(-z²)/(z-ζ);
|
||||||
|
! ζ = (ω - nω_c)/(√2 k∥ v_t)
|
||||||
|
! = (1 - nY)/(√(2Θ)⋅N∥).
|
||||||
|
!
|
||||||
|
! Since the plasma is weakly dissipative, the wavenumber has a small
|
||||||
|
! positive imaginary part, so Im(ζ)<0. By the Sokhotski–Plemelj
|
||||||
|
! theorem we have that:
|
||||||
|
!
|
||||||
|
! lim_Im(ζ)→0 Z(ζ) = 1/√π P ∫dz exp(-z²)/(z-ζ) + i√π exp(-ζ²)
|
||||||
|
! = -√π exp(-ζ²) erfi(ζ) + i√π exp(-ζ²)
|
||||||
|
!
|
||||||
|
! where erfi(ζ) = -i⋅erf(iζ). So, Φ(ζ) = √π exp(-ζ²).
|
||||||
|
!
|
||||||
|
use const_and_precisions, only : sqrt_pi
|
||||||
|
|
||||||
|
implicit none
|
||||||
|
|
||||||
|
! function arguments
|
||||||
|
|
||||||
|
! CMA Y variable: Y=ω_c/ω
|
||||||
|
real(wp_), intent(in) :: Y
|
||||||
|
! squared parallel refractive index: N∥²=N²cos²θ
|
||||||
|
real(wp_), intent(in) :: Npl2
|
||||||
|
! adimensional temperature: Θ = kT/mc²
|
||||||
|
real(wp_), intent(in) :: theta
|
||||||
|
! harmonic number
|
||||||
|
integer, intent(in) :: n
|
||||||
|
! line shape Φ(ζ)
|
||||||
|
real(wp_) :: phi
|
||||||
|
|
||||||
|
! local variables
|
||||||
|
real(wp_) :: z2
|
||||||
|
|
||||||
|
z2 = (1 - n*Y)**2/(2 * theta * Npl2)
|
||||||
|
phi = sqrt_pi * exp(-z2)
|
||||||
|
end function line_shape_harmonics
|
||||||
|
|
||||||
|
|
||||||
|
pure function warm_index(d, Npl2, Npr2, sox) result(M2)
|
||||||
|
! The real part of the warm plasma refractive index
|
||||||
|
! in the finite density regime near the first harmonic:
|
||||||
|
!
|
||||||
|
! M² = (1 - δ + ½ N⊥²/N² ∓ ½ √Δ) ⋅ N²/N⊥²
|
||||||
|
!
|
||||||
|
! where Δ = N⊥⁴/N⁴ + 4n²(1-δ)²N∥²/N²;
|
||||||
|
! δ = X/Y².
|
||||||
|
!
|
||||||
|
! Reference: 1983 Nuclear Fusion 23 1153 - table VII
|
||||||
|
!
|
||||||
|
implicit none
|
||||||
|
|
||||||
|
! function arguments
|
||||||
|
|
||||||
|
! parameter δ=(ω_pe/ω_ce)²
|
||||||
|
real(wp_), intent(in) :: d
|
||||||
|
! squared refractive index: N∥², N⊥² (cold dispersion)
|
||||||
|
real(wp_), intent(in) :: Npl2, Npr2
|
||||||
|
! sign of polarisation mode: -1 ⇒ O, +1 ⇒ X
|
||||||
|
integer, intent(in) :: sox
|
||||||
|
! warm refactive indexL M²
|
||||||
|
real(wp_) :: M2
|
||||||
|
|
||||||
|
! local variables
|
||||||
|
real(wp_) :: N2, delta
|
||||||
|
|
||||||
|
! total refractive index N²
|
||||||
|
N2 = Npl2 + Npr2
|
||||||
|
|
||||||
|
delta = (Npr2/N2)**2 + 4*(1-d)**2 * (Npl2/N2)
|
||||||
|
M2 = (1 - d + Npr2/N2/2 + sox*sqrt(delta)/2) * N2/Npr2
|
||||||
|
end function warm_index
|
||||||
|
|
||||||
|
|
||||||
|
pure function warm_index_harmonics(d, Npl2, Npr2, sox, n) result(M2)
|
||||||
|
! The real part of the warm plasma refractive index
|
||||||
|
! in the finite density regime near the harmonics (n>2)
|
||||||
|
!
|
||||||
|
! M² = 1 - δ⋅f
|
||||||
|
!
|
||||||
|
! where f = (1 - δ) / [(1 - δ) - (N⊥²/N² ∓ ρ) / (2n²)];
|
||||||
|
! ρ² = N⊥⁴/N⁴ + 4n²⋅(1 - δ)²⋅ N∥²/N²;
|
||||||
|
! δ = X/(nY)².
|
||||||
|
!
|
||||||
|
! Reference: 1983 Nuclear Fusion 23 1153 - table VIII
|
||||||
|
!
|
||||||
|
implicit none
|
||||||
|
|
||||||
|
! function arguments
|
||||||
|
|
||||||
|
! parameter δ=(ω_p/nω_c)²
|
||||||
|
real(wp_), intent(in) :: d
|
||||||
|
! squared refractive index: N∥², N⊥² (cold dispersion)
|
||||||
|
real(wp_), intent(in) :: Npl2, Npr2
|
||||||
|
! sign of polarisation mode: -1 ⇒ O, +1 ⇒ X
|
||||||
|
integer, intent(in) :: sox
|
||||||
|
! harmonic number
|
||||||
|
integer, intent(in) :: n
|
||||||
|
! warm refactive index
|
||||||
|
real(wp_) :: M2
|
||||||
|
|
||||||
|
! local variables
|
||||||
|
real(wp_) :: N2, f, r2
|
||||||
|
|
||||||
|
! total refractive index N²
|
||||||
|
N2 = Npl2 + Npr2
|
||||||
|
|
||||||
|
r2 = (Npr2/N2)**2 + 4*n**2 * (1 - d)**2 * Npl2/N2
|
||||||
|
f = (1 - d) / ((1 - d) - (Npr2/N2 + sox * sqrt(r2)) / (2*n**2))
|
||||||
|
M2 = 1 - d*f
|
||||||
|
end function warm_index_harmonics
|
||||||
|
|
||||||
|
|
||||||
|
pure function polarisation(X, Y, N2, Npl2, Npr2, sox) result(R)
|
||||||
|
! The polarisation function R(θ,O/X) in the finite density regime
|
||||||
|
! for the first harmonic:
|
||||||
|
!
|
||||||
|
! R(N∥, N⊥, O/X) = [(1 - δ)⋅(1 - δ/2 - M∥²) - M⊥²]²
|
||||||
|
! ⋅ 1/√(a² + b²) ⋅ M∥ ⋅ 1/d
|
||||||
|
!
|
||||||
|
! where δ = X/Y²
|
||||||
|
! a² = [(1 - δ - M⊥²)² + (1 - d)⋅M∥²]² ⋅ M⊥²
|
||||||
|
! b² = (2 - 2δ - M⊥²)² ⋅ (1 - d - M⊥²)² ⋅ M∥²
|
||||||
|
! M is the real part of the warm refractive index
|
||||||
|
!
|
||||||
|
! Reference: 1983 Nuclear Fusion 23 1153 - formula (3.1.68)
|
||||||
|
!
|
||||||
|
implicit none
|
||||||
|
|
||||||
|
! function arguments
|
||||||
|
|
||||||
|
! CMA diagram variables: X=(ω_pe/ω)², Y=ω_ce/ω
|
||||||
|
real(wp_), intent(in) :: X, Y
|
||||||
|
! squared refractive index: N², N∥², N⊥² (cold dispersion)
|
||||||
|
real(wp_), intent(in) :: N2, Npl2, Npr2
|
||||||
|
! sign of polarisation mode: -1 ⇒ O, +1 ⇒ X
|
||||||
|
integer, intent(in) :: sox
|
||||||
|
! polarisation function
|
||||||
|
real(wp_) :: R
|
||||||
|
|
||||||
|
! local variables
|
||||||
|
real(wp_) :: a2, b2, d
|
||||||
|
real(wp_) :: M2, Mpl2, Mpr2
|
||||||
|
|
||||||
|
d = X/Y**2 ! δ=(ω_pe/ω_ce)²
|
||||||
|
|
||||||
|
M2 = warm_index(d, Npl2, Npr2, sox)
|
||||||
|
Mpl2 = M2 * Npl2/N2 ! M∥² = M²cos²θ
|
||||||
|
Mpr2 = M2 * Npr2/N2 ! M⊥² = M²sin²θ
|
||||||
|
|
||||||
|
a2 = ((1 - d - Mpr2)**2 + (1 - d)*Mpl2)**2 * Mpr2
|
||||||
|
b2 = (2 - 2*d - Mpr2)**2 * (1 - d - Mpr2)**2 * Mpl2
|
||||||
|
R = ((1 - d)*(1 - d/2 - Mpl2) - Mpr2)**2 / sqrt(a2 + b2) &
|
||||||
|
* sqrt(Mpl2) / d
|
||||||
|
end function polarisation
|
||||||
|
|
||||||
|
|
||||||
|
pure function polarisation_harmonics(X, Y, N2, Npl2, Npr2, sox, n) result(mu)
|
||||||
|
! The polarisation function μ(n,θ,O/X) times 2(1+cos²θ)
|
||||||
|
! in the finite-density plasma regime for harmonics (n>1):
|
||||||
|
!
|
||||||
|
! μ⋅2(1+cos²θ) = M^(2n-3) (n-1)² (1 - (1+1/n)⋅f) / √(a²+b²)
|
||||||
|
!
|
||||||
|
! where f as in `warm_index_harmonics`;
|
||||||
|
! a = (N⊥/N) (1 + c/e² M∥²);
|
||||||
|
! b = (N∥/N) (1 + c/e);
|
||||||
|
! c = (1-δ) n²[1 - (1-1/n²)⋅f]²;
|
||||||
|
! e = 1-δ-M⊥²
|
||||||
|
! M∥ = M/N N∥ = M cosθ;
|
||||||
|
! M⊥ = M/N N⊥ = M sinθ;
|
||||||
|
! M is the real part of the warm refractive index.
|
||||||
|
!
|
||||||
|
! Reference: 1983 Nuclear Fusion 23 1153 - formula (3.1.77)
|
||||||
|
!
|
||||||
|
implicit none
|
||||||
|
|
||||||
|
! function arguments
|
||||||
|
|
||||||
|
! CMA diagram variables: X=(ω_pe/ω)², Y=ω_ce/ω
|
||||||
|
real(wp_), intent(in) :: X, Y
|
||||||
|
! squared refractive index: N², N∥², N⊥² (cold dispersion)
|
||||||
|
real(wp_), intent(in) :: N2, Npl2, Npr2
|
||||||
|
! sign of polarisation mode: -1 ⇒ O, +1 ⇒ X
|
||||||
|
integer, intent(in) :: sox
|
||||||
|
! harmonic number
|
||||||
|
integer, intent(in) :: n
|
||||||
|
! polarisation function
|
||||||
|
real(wp_) :: mu
|
||||||
|
|
||||||
|
! local variables
|
||||||
|
real(wp_) :: d, f, a2, b2, c, e
|
||||||
|
real(wp_) :: M2, Mpl2, Mpr2
|
||||||
|
|
||||||
|
! δ=(ω_pe/nω_ce)²
|
||||||
|
d = X/(n*Y)**2
|
||||||
|
|
||||||
|
M2 = warm_index_harmonics(d, Npl2, Npr2, sox, n)
|
||||||
|
Mpl2 = M2 * Npl2/N2 ! M∥ = Mcosθ
|
||||||
|
Mpr2 = M2 * Npr2/N2 ! M⊥ = Msinθ
|
||||||
|
|
||||||
|
f = (1 - M2)/d
|
||||||
|
c = (1 - d) * n**2 * (1 - (1 - 1/n**2) * f)**2
|
||||||
|
e = 1 - d - Mpr2
|
||||||
|
|
||||||
|
a2 = Npr2/N2 * (1 + c/e**2 * Mpl2)**2 ! a²
|
||||||
|
b2 = Npl2/N2 * (1 + c/e)**2 ! b²
|
||||||
|
|
||||||
|
mu = M2**(n-3.0_wp_/2) * (n-1)**2 &
|
||||||
|
* (1 - (1 + 1/n) * f)**2 &
|
||||||
|
/ sqrt(a2 + b2)
|
||||||
|
end function polarisation_harmonics
|
||||||
|
|
||||||
|
|
||||||
|
function alpha_fast(X, Y, theta, k0, Npl, Npr, sox, &
|
||||||
|
nhmin, nhmax) result(alpha)
|
||||||
|
! The absorption coefficient around ω=nω_c in the low
|
||||||
|
! temperature (< 10keV), low harmonic number limit (< 6),
|
||||||
|
! and oblique propagation:
|
||||||
|
!
|
||||||
|
! α(n,ω,θ,O/X) = α(n,θ)⋅μ(n,θ,O/X)⋅Φ(n,ω)
|
||||||
|
!
|
||||||
|
! where α(n,θ) is the nth-harmonic absorption coeff.;
|
||||||
|
! μ(n,θ,O/X) is the polarisation part;
|
||||||
|
! Φ(n,ω) is the line shape function.
|
||||||
|
!
|
||||||
|
! Finally, for the X1 mode:
|
||||||
|
!
|
||||||
|
! Notes:
|
||||||
|
! 1. The approximation is decent for
|
||||||
|
! δ = (ω_pe/nω_ce)² / (2 √Θ N∥) << 1 (tenuous plasma)
|
||||||
|
! >> 1 (finite density);
|
||||||
|
! 2. In the implementation the angular dependencies
|
||||||
|
! have been replaced by cosθ=N∥/N, sinθ=N⊥/N;
|
||||||
|
! 3. The special cases for n=1 have been recast
|
||||||
|
! from the original form (tbl. XI, eq. 3.1.67)
|
||||||
|
! to match with eq. 3.1.73.
|
||||||
|
!
|
||||||
|
implicit none
|
||||||
|
|
||||||
|
! function arguments
|
||||||
|
|
||||||
|
! CMA diagram variables: X=(ω_pe/ω)², Y=ω_ce/ω
|
||||||
|
real(wp_), intent(in) :: X, Y
|
||||||
|
! adimensional temperature: Θ = kT/mc²
|
||||||
|
real(wp_), intent(in) :: theta
|
||||||
|
! vacuum wavenumber k₀=ω/c
|
||||||
|
real(wp_), intent(in) :: k0
|
||||||
|
! refractive index: N∥, N⊥ (cold dispersion)
|
||||||
|
real(wp_), intent(in) :: Npl, Npr
|
||||||
|
! sign of polarisation mode: -1 ⇒ O, +1 ⇒ X
|
||||||
|
integer, intent(in) :: sox
|
||||||
|
! min/max harmonic number
|
||||||
|
integer, intent(in) :: nhmin, nhmax
|
||||||
|
! absorption coefficient
|
||||||
|
real(wp_) :: alpha
|
||||||
|
|
||||||
|
! local variables
|
||||||
|
integer :: n ! harmonic number
|
||||||
|
real(wp_) :: N2, Npr2, Npl2
|
||||||
|
real(wp_) :: delta
|
||||||
|
real(wp_) :: shape, polar, absorp
|
||||||
|
|
||||||
|
Npr2 = Npr**2 ! N⊥² refractive index
|
||||||
|
Npl2 = Npl**2 ! N∥²
|
||||||
|
N2 = Npr2 + Npl2 ! N²
|
||||||
|
|
||||||
|
! tenuous/finite density parameter:
|
||||||
|
! δ = (ω_p/ω_c)²/(2√Θ N∥)
|
||||||
|
delta = X/Y**2 / (2*abs(Npl)*sqrt(theta))
|
||||||
|
|
||||||
|
! First harmonic
|
||||||
|
if (delta < 1) then
|
||||||
|
! Tenuous plasma
|
||||||
|
if (sox < 0) then
|
||||||
|
! O mode
|
||||||
|
! absorption coeff: k₀⋅X
|
||||||
|
absorp = k0 * X
|
||||||
|
|
||||||
|
! polarisation part: N⊥⁴/N² ⋅ (N²+2N∥²)²/(N²+ N∥²)³
|
||||||
|
polar = Npr2**2/N2 * (N2 + 2*Npl2)**2/(N2 + Npl2)**3
|
||||||
|
|
||||||
|
! line shape part: Φ(ζ) = Im(-1/Z(ζ)) / (√(2/Θ)⋅Y⋅|N∥|)
|
||||||
|
shape = line_shape(Y, Npl2, theta) &
|
||||||
|
/ (sqrt(2/theta) * Y * abs(Npl))
|
||||||
|
else
|
||||||
|
! X mode
|
||||||
|
! polarisation part: (N² + N∥²)/N
|
||||||
|
polar = (N2 + Npl2) / sqrt(N2)
|
||||||
|
|
||||||
|
! line shape part: Φ(ζ) = Im(Z(ζ)) / (2√(2θ)⋅Y⋅|N∥|)
|
||||||
|
shape = line_shape_harmonics(Y, Npl2, theta, n) &
|
||||||
|
/ (2 * sqrt(2*theta) * Y * abs(Npl))
|
||||||
|
end if
|
||||||
|
else
|
||||||
|
! Finite density (δ > 1)
|
||||||
|
|
||||||
|
! absorption coeff: k₀⋅Y
|
||||||
|
absorp = k0 * Y
|
||||||
|
|
||||||
|
! polarisation part: R(N∥, N⊥, O/X)
|
||||||
|
polar = polarisation(X, Y, N2, Npl2, Npr2, sox)
|
||||||
|
|
||||||
|
! line shape part: Φ(ζ) = Im(-1/Z(ζ)) ⋅ 2√(2Θ)
|
||||||
|
shape = line_shape(Y, Npl2, theta) &
|
||||||
|
* 2 * sqrt(2*theta)
|
||||||
|
end if
|
||||||
|
alpha = absorp * polar * shape
|
||||||
|
|
||||||
|
if (nhmax == 1) return
|
||||||
|
|
||||||
|
! n-independent part of α: k₀⋅X
|
||||||
|
absorp = k0 * X
|
||||||
|
|
||||||
|
! Higher harmonics
|
||||||
|
! note: the loop starts from n=1 to compute
|
||||||
|
! the powers/factorials of n, even when nhmin>1.
|
||||||
|
do n=1,nhmax
|
||||||
|
! absorption part: α(n,θ)/[(1+cos²θ)⋅n^(2n)⋅(πω)]
|
||||||
|
absorp = absorp / (2*n) ! computes (2n)!!
|
||||||
|
|
||||||
|
if (n >= nhmin) then
|
||||||
|
! polarisation part: (1+cos²θ)⋅μ(n,θ)
|
||||||
|
polar = polarisation_harmonics(X, Y, N2, Npl2, Npr2, sox, n)
|
||||||
|
|
||||||
|
! line shape part: Φ(n,ω)⋅(πω)
|
||||||
|
shape = line_shape_harmonics(Y, Npl2, theta, n) &
|
||||||
|
/ (sqrt(2*theta) * n*Y * abs(Npl))
|
||||||
|
|
||||||
|
! full absorption coefficient: α(n,ω,θ)
|
||||||
|
alpha = alpha + absorp*n**(2*n) * polar * shape
|
||||||
|
end if
|
||||||
|
|
||||||
|
! computes * (Θ N⊥²/N²)^(n-1)
|
||||||
|
absorp = absorp * (theta * Npr2/N2)
|
||||||
|
end do
|
||||||
|
end function alpha_fast
|
||||||
|
|
||||||
|
end module absorption
|
||||||
@ -52,6 +52,7 @@ contains
|
|||||||
integer, intent(out) :: error
|
integer, intent(out) :: error
|
||||||
|
|
||||||
! local variables
|
! local variables
|
||||||
|
! taucr: max τ before considering a ray fully absorbed
|
||||||
real(wp_), parameter :: taucr = 12._wp_, etaucr = exp(-taucr)
|
real(wp_), parameter :: taucr = 12._wp_, etaucr = exp(-taucr)
|
||||||
character, dimension(2), parameter :: mode=['O','X']
|
character, dimension(2), parameter :: mode=['O','X']
|
||||||
|
|
||||||
@ -1985,7 +1986,7 @@ contains
|
|||||||
real(wp_), intent(in) :: psinv
|
real(wp_), intent(in) :: psinv
|
||||||
! CMA diagram variables: X=(ω_pe/ω)², Y=ω_ce/ω
|
! CMA diagram variables: X=(ω_pe/ω)², Y=ω_ce/ω
|
||||||
real(wp_), intent(in) :: X, Y
|
real(wp_), intent(in) :: X, Y
|
||||||
! densityity [10¹⁹ m⁻³], temperature [keV]
|
! density [10¹⁹ m⁻³], temperature [keV]
|
||||||
real(wp_), intent(in) :: density, temperature
|
real(wp_), intent(in) :: density, temperature
|
||||||
! vacuum wavenumber k₀=ω/c, resonant B field
|
! vacuum wavenumber k₀=ω/c, resonant B field
|
||||||
real(wp_), intent(in) :: k0, Bres
|
real(wp_), intent(in) :: k0, Bres
|
||||||
@ -2049,28 +2050,48 @@ contains
|
|||||||
end if
|
end if
|
||||||
end if
|
end if
|
||||||
|
|
||||||
! Compute α from the solution of the dispersion relation
|
if (params%iwarm < 4) then
|
||||||
! The absoption coefficient is defined as
|
! Compute α from the solution of the dispersion relation
|
||||||
!
|
accurate: block
|
||||||
! α = 2 Im(k̅)⋅s̅
|
! Compute α from the solution of the dispersion relation
|
||||||
!
|
! The absoption coefficient is defined as
|
||||||
! where s̅ = v̅_g/|v_g|, the direction of the energy flow.
|
!
|
||||||
! Since v̅_g = - ∂Λ/∂N̅ / ∂Λ/∂ω, using the cold dispersion
|
! α = 2 Im(k̅)⋅s̅
|
||||||
! relation, we have that
|
!
|
||||||
!
|
! where s̅ = v̅_g/|v_g|, the direction of the energy flow.
|
||||||
! s̅ = ∂Λ/∂N̅ / |∂Λ/∂N̅|
|
! Since v̅_g = - ∂Λ/∂N̅ / ∂Λ/∂ω, using the cold dispersion
|
||||||
! = [2N̅ - ∂(N²s)/∂N∥ b̅
|
! relation, we have that
|
||||||
! + ½(b̅⋅∇S_I)² ∂³(N²s)/∂N∥³ b̅] / |∂Λ/∂N̅|
|
!
|
||||||
!
|
! s̅ = ∂Λ/∂N̅ / |∂Λ/∂N̅|
|
||||||
! Assuming Im(k∥)=0:
|
! = [2N̅ - ∂(N²s)/∂N∥ b̅
|
||||||
!
|
! + ½(b̅⋅∇S_I)² ∂³(N²s)/∂N∥³ b̅] / |∂Λ/∂N̅|
|
||||||
! α = 4 Im(k⊥)⋅N⊥ / |∂Λ/∂N̅|
|
!
|
||||||
!
|
! Assuming Im(k∥)=0:
|
||||||
block
|
!
|
||||||
real(wp_) :: k_im
|
! α = 4 Im(k⊥)⋅N⊥ / |∂Λ/∂N̅|
|
||||||
k_im = k0 * Npr%im ! imaginary part of k⊥
|
!
|
||||||
alpha = 4 * k_im*Npr_cold / derdnm
|
real(wp_) :: k_im
|
||||||
end block
|
k_im = k0 * Npr%im ! imaginary part of k⊥
|
||||||
|
alpha = 4 * k_im*Npr_cold / derdnm
|
||||||
|
end block accurate
|
||||||
|
else
|
||||||
|
! Compute a fast and approximate α
|
||||||
|
fast: block
|
||||||
|
use logger, only: log_warning
|
||||||
|
use absorption, only: alpha_fast
|
||||||
|
|
||||||
|
character(256) :: msg
|
||||||
|
|
||||||
|
! the absorption coefficient in the tenuous plasma limit
|
||||||
|
alpha = alpha_fast(X, Y, 1/mu, k0, Npl, Npr_cold, sox, nhmin, nhmax)
|
||||||
|
|
||||||
|
if (temperature > 10 .or. nhmax > 5) then
|
||||||
|
write (msg, '(a,g0.3,a,g0)') &
|
||||||
|
'iwarm=4 is inaccurate: Te=', temperature, ' nhmax=', nhmax
|
||||||
|
call log_warning(msg, mod='gray_core', proc='alpha_effj')
|
||||||
|
end if
|
||||||
|
end block fast
|
||||||
|
end if
|
||||||
|
|
||||||
if (alpha < 0) then
|
if (alpha < 0) then
|
||||||
error = raise_error(error, negative_absorption)
|
error = raise_error(error, negative_absorption)
|
||||||
|
|||||||
Loading…
Reference in New Issue
Block a user