From 0a63a20e734a945e3dada46b69c4925a0a2b8bc8 Mon Sep 17 00:00:00 2001 From: Michele Guerini Rocco Date: Fri, 20 May 2022 13:10:49 +0200 Subject: [PATCH] src/dispersion.f90: cleanup - merge branch with a method to control the speed of iteration and improve the convergence of `warmdisp` (thanks Thomas) - unify `diel_tens_fr` and `diel_tens_wr` into a single subroutine, `dielectric_tensor` - stay as close as possible to the notation of Daniela Farina's paper - make `sox` an integer - mark more subroutines as pure - add more comments --- src/dispersion.f90 | 1055 ++++++++++++++++++++++++------------------ src/eccd.f90 | 16 +- src/gray_core.f90 | 144 +++--- src/multipass.f90 | 9 +- src/polarization.f90 | 3 +- 5 files changed, 708 insertions(+), 519 deletions(-) diff --git a/src/dispersion.f90 b/src/dispersion.f90 index 1403a97..30a1c8e 100644 --- a/src/dispersion.f90 +++ b/src/dispersion.f90 @@ -1,416 +1,668 @@ +! This module contains subroutines related to the dispersion relation +! for EC waves in a plasma. +! +! In particular, +! `colddisp` gives the analytical solution of the cold plasma +! dispersion; +! `warmdisp` solves numerically the warm plasma dispersion, either in +! full or in the weakly relativistic approximation; +! `harmnumber` computes the possible harmonic numbers and can be used +! to limit the number of calls to `warmdisp`. +! +! Note: the module must be initialised by calling `expinit` when using +! the fully relativistic model (argument `model` > 1). module dispersion -! - use const_and_precisions, only : wp_,zero,one,half,two,im,czero,cunit,pi, & - sqrt_pi + + use const_and_precisions, only : wp_, zero, one, half, two, & + im, czero, cunit, pi, sqrt_pi implicit none -! local constants - integer, parameter :: npts=500 - real(wp_), parameter :: tmax=5.0_wp_,dtex=2.0_wp_*tmax/dble(npts) -! variables - real(wp_), dimension(npts+1), save :: ttv,extv -! + + ! global constants + integer, parameter :: npts = 500 + real(wp_), parameter :: tmax = 5.0_wp_ + real(wp_), parameter :: dtex = 2*tmax/dble(npts) + + ! global variables + real(wp_), dimension(npts+1), save :: ttv, extv + + private + public expinit, colddisp, warmdisp + public zetac, harmnumber + contains -pure function colddisp(xg, yg, npl, sox) result(npr) +pure function colddisp(X, Y, Npl, sox) result(npr) ! Cold plasma dispersion relation - - ! Given the parallel component of refractive index returns - ! the orthogonal one. + ! + ! Given the parallel component of refractive index N∥ + ! returns the orthogonal one N⊥. ! ! Reference: IFP-CNR Internal Report FP 05/1 - App. A - ! implicit none ! subroutine arguments ! CMA diagram variables: X=(ω_pe/ω)², Y=ω_ce/ω - real(wp_), intent(in) :: xg, yg + real(wp_), intent(in) :: X, Y ! cold refractive index: N∥, N⊥ - real(wp_), intent(in) :: npl - real(wp_) :: npr + real(wp_), intent(in) :: Npl + real(wp_) :: Npr ! sign of polarisation mode: -1 ⇒ O, +1 ⇒ X - real(wp_), intent(in) :: sox + integer, intent(in) :: sox ! local variables - real(wp_) :: yg2, npl2, dnl, del, dxg + real(wp_) :: Y2, Npl2, delta, dX - npl2 = npl**2 - yg2 = yg**2 - dnl = one - npl2 - dxg = one - xg - del = sqrt(dnl**2 + 4.0_wp_*npl2*dxg/yg2) - npr = sqrt(dxg - npl2 - xg*yg2 * (one + npl2 + sox*del)/(dxg - yg2)/2.0_wp_) + ! The formula is: + ! + ! N⊥ = √(N² - N∥²) + ! + ! where + ! N² = 1 - X - ½ XY² (1 + N∥² ± Δ)/(1 - X - Y²) + ! Δ² = (1 - N∥²)² + 4N∥²(1 - X)/Y² + ! + Npl2 = Npl**2 + Y2 = Y**2 + dX = 1 - X + delta = sqrt((1 - Npl2)**2 + 4*Npl2*dX/Y2) + Npr = sqrt(dx - Npl2 - X*Y2 * (1 + Npl2 + sox*delta)/(dX - Y2)/2) end function colddisp -! -subroutine harmnumber(yg,mu,npl,nhmin,nhmax,iwr) -! computation of minimum and maximum harmonic -! - implicit none -! local constants -! expcr = maximum value for mu*(gamma-1) above which the distribution function -! is considered to be 0 - real(wp_), parameter :: expcr=16.0_wp_ -! arguments -! yg = omegac/omega -! npl = parallel N -! mu = mc^2/Te -! nh = number of the armonic (min/max) -! iwr = weakly (iwr=1) or fully relativistic approximation - integer :: nhmin,nhmax,iwr - real(wp_) :: yg,npl,mu -! local variables - integer :: nh,nhc - real(wp_) :: ygc,ygn,npl2,gg,dnl,rdu2,argexp,uu2 -! - nhmin=0 - nhmax=0 - npl2=npl**2 - dnl=one-npl2 -! - if(iwr.eq.1) then - ygc=max(one-0.5_wp_*npl2,zero) - else - ygc=sqrt(max(dnl,zero)) - end if - nhc=int(ygc/yg) - if (nhc*yg zero) then - s = sox - else - s = -sox - end if - - npr2 = (s*rdiscr - half*cc2)/cc4 - - errnpr=abs(one-abs(npr2)/abs(npra2)) - if(i.gt.1.and.errnpr.lt.1.0e-5_wp_) exit - npra2=npr2 - end do iterations_loop - - if(i.gt.imxx.and.imxx.gt.1) then - if (imx.lt.0) then - ! first try failed - err=1 - imxx=1 - npra2=nprf**2 - else - ! first try failed, no retry - err=2 - exit - end if - else - ! converged or second try - exit - end if - - end do tries_loop - - if(dble(sqrt(npr2)).lt.zero.or.npr2.ne.npr2.or.abs(npr2).ge.huge(one).or. & - abs(npr2).le.tiny(one)) then - npr2=czero - err=err+4 - end if -! - npr=sqrt(npr2) - nprr=dble(npr) - npri=dimag(npr) -! - ex=czero - ey=czero - ez=czero -! - if (abs(npl).gt.1.0e-6_wp_) then - den=e12*e23-(e13+npr*npl)*(e22-npr2-npl2) - ey=-(e12*(e13+npr*npl)+(e11-npl2)*e23)/den - ez=(e12*e12+(e22-npr2-npl2)*(e11-npl2))/den - ez2=abs(ez)**2 - ey2=abs(ey)**2 - enx2=one/(one+ey2+ez2) - ex=dcmplx(sqrt(enx2),zero) - ez2=ez2*enx2 - ey2=ey2*enx2 - ez=ez*ex - ey=ey*ex - else - if(sox.lt.zero) then - ez=cunit - ez2=abs(ez)**2 - else - ex2=one/(one+abs(-e11/e12)**2) - ex=sqrt(ex2) - ey=-ex*e11/e12 - ey2=abs(ey)**2 - ez2=zero - end if - end if -! -end subroutine warmdisp - - -subroutine diel_tens_fr(yg,mu,npl,a330,epsl,lrm,fast) - ! Fully relativistic case computation of dielectric tensor elements - ! up to third order in Larmor radius for hermitian part +pure subroutine harmnumber(Y, mu, Npl2, weakly, nhmin, nhmax) + ! Computes the min/max possible harmonic numbers ! + ! The condition for the wave-particle resonance to take place at + ! the nth harmonic (including relativity and Doppler effects) is: + ! + ! γ(ω - k∥⋅v∥) = n ω_ce + ! + ! where ω_ce = eB/m, γ = 1/√(1 - β²). + ! Defining the adimensional velocity u̅ = p̅/mc = γβ̅, the Lorentz + ! factor can be written as + ! + ! γ² = 1 + γ²β² = 1 + u∥² + u⊥² + ! + ! Substituting this expression for γ, the condition can be rewritten + ! as the equation for an ellipse in the (u∥, u⊥) plane: + ! + ! (u∥ - u∥₀)²/α∥² + u⊥²/α⊥² = 1 + ! + ! where u∥₀ = nYN∥/(1 - N∥²) + ! α∥ = √(N∥² + n²Y² - 1)/(1 - N∥²) + ! α⊥ = √(N∥² + n²Y² - 1)/√(1 - N∥²) + ! Y = ω_ce/ω + ! N̅ = ck̅/ω + ! Assuming thermal equilibrium, the electrons are distributed + ! according to a Boltzmann factor + ! + ! exp(-E/kT) = exp(-mc²(γ-1)/kT) + ! = exp(-μ⋅√(1+u²)), (μ=mc²/kT) + ! + ! which is constant along circles centered at the origin: + ! + ! u∥² + u⊥² = u² + ! + ! A harmonic number is possible when its ellipse intersects a region + ! containing a significant fraction of the electron population. + ! + implicit none + + ! subroutine arguments + + ! CMA Y variable: Y=ω_c/ω + real(wp_), intent(in) :: Y + ! squared parallel refractive index: N∥²=N²cos²θ + real(wp_), intent(in) :: Npl2 + ! reciprocal of adimensional temperature: μ=mc²/kT + real(wp_), intent(in) :: mu + ! whether to use the weakly-relativistic limit + logical, intent(in) :: weakly + ! min, max harmonic numbers + integer, intent(out) :: nhmin, nhmax + + ! local variables + integer :: nh, nhc + real(wp_) :: Yc, Yn, gamma, rdu2, argexp, uu2 + + ! local constants + + ! Maximum value for μ⋅(γ-1) above which the + ! distribution function is considered to be 0. + real(wp_), parameter :: expcr = 16.0_wp_ + + nhmin = 0 + nhmax = 0 + + ! Compute the critical Y value for which the first + ! harmonic ellipse collapses into a single point. + ! + ! The condition is α∥⋅α⊥=0 (with n=1), from which: + ! + ! Yc = √(1 - N∥²) (fully relativistic) + ! Yc = 1 -½ N∥² (weakly relativistic) + ! + if (weakly) then + Yc = max(one - npl2/2, zero) + else + Yc = sqrt(max(1 - npl2, zero)) + end if + + ! At the given Y and Yc, the "critical" harmonic + ! numeber, i.e. the first non-collapsed, is ⌈Yc/Y⌉. + ! Note: to see why consider the axes can be written + ! + ! α∥ = √(n²Y² - Yc²)/Yc² + ! α⊥ = √(n²Y² - Yc²)/Yc + ! + nhc = ceiling(Yc/Y) + + ! Check a few numbers starting with nhc + do nh = nhc, nhc + 10 + Yn = Y*nh + + ! Find the closest point to the origin + ! + ! The solutions are simply: + ! + ! u⊥ = 0 + ! u∥ = u∥₀ - sgn(N∥)⋅α∥ + ! + ! Note that N∥ can be negative, in which case u∥₀ is + ! in the second quadrant and we must swap the sign. + ! The value of γ in this point is + ! + ! γ² = 1 + u∥² + ! = 1 + u∥₀² + α∥² - 2⋅sgn(N∥)⋅u∥₀α∥ + ! + ! The sign can be simplified by noting that + ! + ! sgn(N∥)⋅u∥₀ = sgn(N∥)⋅N∥ nY/(1 - N∥²) + ! = |N∥|⋅Yn/Yc² + ! + ! Putting everything back together: + ! + ! γ = √[1 + u∥₀² + α∥² - 2⋅sgn(N∥)⋅u∥₀α∥] + ! = √[Yc⁴ + N∥²Yn² + (Yn²-Yc²) - 2⋅Yn √(N∥²(Yn²-Yc²))]/Yc² + ! = √[Yn² + N∥²(Yn² - Yc²) - 2⋅Yn √(N∥²(Yn²-Yc²))]/Yc² + ! = |Yn - √(N∥²(Yn²-Yc²))|/Yc² + ! + if (weakly) then + rdu2 = 2*(Yn - Yc) + uu2 = npl2 + rdu2 - 2*sqrt(Npl2*rdu2) + argexp = mu*uu2/2 + else + rdu2 = Yn**2 - Yc**2 + gamma = (Yn - sqrt(Npl2*rdu2))/(1 - Npl2) + argexp = mu*(gamma - one) + end if + + if (argexp <= expcr) then + ! The are enough electrons with this energy, + ! update the max harmonic number + nhmax = nh + ! Also update the min if this is the first check + if (nhmin == 0) nhmin = nh + else + ! Too few electrons, exit now as higher harmonics + ! will have even fewer particles than the current + if (nhmin > 0) exit + end if + + end do + +end subroutine harmnumber + + +subroutine warmdisp(X, Y, mu, Npl, Npr_cold, sox, & + error, Npr, e, & + model, nlarmor, max_iters, fallback) + ! Solves the warm plasma dispersion relation using an iterative + ! method starting from the cold plasma solution + ! + ! Reference: https://doi.org/10.13182/FST08-A1660 + use, intrinsic :: ieee_arithmetic, only: ieee_is_finite implicit none ! subroutine arguments - integer :: lrm,fast - real(wp_) :: yg,mu,npl - complex(wp_) :: a330 - complex(wp_), dimension(3,3,lrm) :: epsl + + ! Inputs + + ! CMA diagram variables: X=(ω_pe/ω)², Y=ω_ce/ω + real(wp_), intent(in) :: X, Y + ! reciprocal of adimensional temperature: μ=mc²/kT + real(wp_), intent(in) :: mu + ! solution of the cold plasma dispersion: N⊥ + real(wp_), intent(in) :: Npr_cold + ! parallel refractive index: N∥ + real(wp_), intent(in) :: Npl + ! sign of polarisation mode: -1 ⇒ O, +1 ⇒ X + integer, intent(in) :: sox + + ! Outputs + + ! error code: + ! 0: converged within `max_iters` + ! 1: failed to converge, returned fallback value + ! 2: failed to converge, returned last value + ! 4: converged to invalid value + ! 5: failed to converge, fallback is also invalid + ! 6: failed to converge, last value is invalid + ! Note: invalid means either IEEE 754 NaN or ±Infinity. + integer, intent(out) :: error + ! orthogonal refractive index: N⊥ + complex(wp_), intent(out) :: Npr + ! polarization vector: ε (local cartesian frame, z∥B̅) + complex(wp_), intent(out), optional :: e(3) + + ! Parameters + + ! switch for the dielectic tensor model: + ! 1: weakly relativistic + ! 2: fully relativistic (faster variant) + ! 3: fully relativistic (slower variant) + integer, intent(in) :: model + ! order of the electron Larmor radius expansion + integer, intent(in) :: nlarmor + ! max number of iterations + integer, intent(in) :: max_iters + ! whether to fall back to the result of the first + ! iteration in case the method fails to converge + logical, intent(in) :: fallback ! local variables - integer :: i,j,l,is,k,lm - real(wp_) :: cr,ci - real(wp_) :: npl2,dnl,w,asl,bsl,cmxw,fal - real(wp_), dimension(-lrm:lrm,0:2,0:lrm) :: rr - real(wp_), dimension(lrm,0:2,lrm) :: ri - complex(wp_) :: ca11,ca12,ca22,ca13,ca23,ca33,cq0p,cq0m,cq1p,cq1m,cq2p - npl2=npl**2 - dnl=one-npl2 + ! iteration counter + integer :: i - cmxw=one+15.0_wp_/(8.0_wp_*mu)+105.0_wp_/(128.0_wp_*mu**2) - cr=-mu*mu/(sqrt_pi*cmxw) - ci=sqrt(2.0_wp_*pi*mu)*mu**2/cmxw + ! squares of the refractive index + complex(wp_) :: Npr2 + real(wp_) :: Npl2 - do l=1,lrm - do j=1,3 - do i=1,3 - epsl(i,j,l)=czero + ! full dielectric tensor: ε̅ + complex(wp_) :: eps(3, 3) + ! N⊥-independent parts: ε₃₃⁰, ε̅^(l) + complex(wp_) :: e330, epsl(3, 3, nlarmor) + ! N⊥, N⊥² from the previous iteration + complex(wp_) :: Npr_prev, Npr2_prev + ! N⊥² from the first iteration, needed when `fallback` is true + complex(wp_) :: Npr2_fallback + ! "error" estimated by the change in successive iterations + real(wp_) :: Npr2_err + ! change in N⊥ between successive iterations + complex(wp_) :: delta_Npr + ! change in complex argument of N⊥ between successive iterations + real(wp_) :: angle, angle_norm + + ! factor for controlling the speed of the iterations + ! Note: the factor is used as follow: + ! + ! N⊥²(i) ← N⊥²(i-1)⋅(1 - speed) + N⊥²(i)⋅speed + ! + real(wp_) :: speed + + ! Use the cold N⊥ as the initial "previous" N⊥, + ! set the ΔN⊥ to a small value + Npr2_prev = Npr_cold**2 + delta_Npr = 1e-10_wp_ + + Npl2 = Npl**2 + Npr2_err = 1 + + ! Compute the N⊥-independent part of the dielectric tensor + call dielectric_tensor(X, Y, mu, Npl, model, nlarmor, e330, epsl) + + iterate_to_solve_disperion: do i = 1, max_iters + Npr_prev = sqrt(Npr2_prev) + + ! Compute the "new" dielectric tensor + ! + ! ε̅ = Σ_l N⊥^(2l - 2) ε̅^(l) (eq. 10) + compute_eps: block + integer :: l + eps = 0 + do l = 1, nlarmor + eps = eps + Npr2_prev**(l-1) * epsl(:,:,l) end do - end do - end do + end block compute_eps - if (fast<4) then - call hermitian(rr,yg,mu,npl,cr,fast,lrm) - else - call hermitian_2(rr,yg,mu,npl,cr,fast,lrm) + ! Compute the speed factor + if (i < 10) then + speed = 1 + elseif (i < 50) then + ! Linear function that is 1 at i=10 and 0.1 at i=50 + speed = 1 - 0.9_wp_*real(i - 10, wp_)/40_wp_ + else + speed = 0.1_wp_ + end if + + ! Check if convergence is reached + ! + ! Make sure we do at least two steps in the iteration + ! and the relative error is within tolerance + if (i > 2 .and. Npr2_err < 1.0e-3_wp_*speed) exit + + solve_quadratic: block + ! Solve the quadratic equation for N⊥² + ! + ! AN⊥⁴ + Bn⊥² + C = 0 (eq. 29) + ! + ! with complex coefficients A,B,C. + integer :: s + complex(wp_) :: A, B, C, delta + + A = (eps(1,1) - Npl2)*(1 - eps(3,3)) + (eps(1,3) + Npl)**2 + B = -(1 - eps(3,3)) * ( (eps(1,1) - Npl2)*(eps(2,2) - Npl2) + eps(1,2)**2 ) & + -e330*(eps(1,1) - Npl2) + 2*eps(1,2)*eps(2,3) * (eps(1,3) + Npl) & + +eps(2,3)**2 * (eps(1,1) - Npl2) - (eps(1,3) + Npl)**2 * (eps(2,2) - Npl2) + C = e330 * ((eps(1,1) - Npl2)*(eps(2,2) - Npl2) + eps(1,2)**2) + + ! polynomial discriminant + delta = B**2 - 4*A*C + + ! Choose the sign of the solution + if (Y > 1) then + s = sox + if (delta%im <= 0) s = -s + else + s = -sox + if (delta%re <= 0 .and. delta%im >= 0) s = -s + end if + + ! Solution + Npr2 = (-B + s*sqrt(delta)) / (2*A) + end block solve_quadratic + + ! A new value for Npr2 has been established: estimate the error + ! and store the new value as the latest prediction, Npr2_prev + Npr2_err = abs(1 - abs(Npr2)/abs(Npr2_prev)) + + ! Apply a sequence of modification to N⊥ to improve the stability + modify_fixed_point: block + ! fixed-point prediction of N⊥ + complex(wp_) :: Npr_fix, change + + ! Save the original fixed-point iteration + Npr_fix = sqrt(Npr2) + + ! If real and imaginary parts of the refractive index have opposite + ! sign, then the wave is a backward wave, typically a Bernstein + ! branch (or we have an instability). We then have to find our way + ! back to the EM branch, or at least make sure this wave number is + ! as emission. + if (Npr_fix%re * Npr_fix%im < 0) Npr2 = conjg(Npr_fix)**2 + + ! Change the speed when the current/previous step go out of phase + if (i > 1) then + change = (Npr_prev - sqrt(Npr2)) / delta_Npr + angle = atan2(change%im, change%re) + angle_norm = modulo(angle/pi + 1, two) - 1 + if (abs(angle_norm) > 0.4_wp_) then + ! Change in speed is linear from 1.0 at abs(angle)=0.4*angle + ! to 0.3 at angle=π + speed = speed * (1 - 0.7_wp_ * (abs(angle_norm) - 0.4_wp_) / (1 - 0.4_wp_)) + end if + end if + + ! In the iteration process we do not want huge steps; do not allow steps > 0.05 + Npr2 = Npr2_prev*(1 - speed) + Npr2*speed + if (abs(Npr2 - Npr2_prev) > 0.05_wp_) then + Npr2 = Npr2_prev + 0.05_wp_ * (Npr2 - Npr2_prev) / abs(Npr2 - Npr2_prev) + ! Again, make sure that we have a damped EM wave and not a + ! Bernstein-like wave (see above) + if (real(sqrt(Npr2)) * aimag(sqrt(Npr2)) < zero) & + Npr2 = conjg(sqrt(Npr2))**2 + end if + end block modify_fixed_point + + ! Store the value of N⊥² and its change between iterations + Npr2_prev = Npr2 + delta_Npr = Npr_prev - sqrt(Npr2_prev) + + ! Store the result of the first iteration as a fallback + ! in case the method fails to converge + if (i == 1 .and. fallback) Npr2_fallback = Npr2 + end do iterate_to_solve_disperion + + ! Handle convergence failures + if (i > max_iters) then + if (fallback) then + ! first try failed, use fallback + error = 1 + Npr2 = Npr2_fallback + else + ! first try failed, no retry + error = 2 + end if + end if + + ! Catch NaN and ±Infinity values + if (.not. ieee_is_finite(Npr2%re) .or. .not. ieee_is_finite(Npr2%im)) then + Npr2 = 0 + error = error + 4 + end if + + if (Npr2%re < 0 .and. Npr2%im < 0) then + Npr2 = sqrt(0.5_wp_) ! Lorenzo, what should we have here?? + error = 99 + end if + + ! Final result + Npr = sqrt(Npr2) + + if (.not. present(e)) return + + ! Compute the polarization vector + ! TODO: find a reference for the formula + polarization: block + complex(wp_) :: den + e = 0 + if (abs(Npl) > 1.0e-6_wp_) then + den = Npr*eps(1,2)*eps(2,3) - Npr*(eps(1,3) + Npl)*(eps(2,2) - Npr2 - Npl2) + e(2) = -(eps(1,2)*Npr*(eps(1,3) + Npl) + (eps(1,1) - Npl2)* Npr*eps(2,3))/den + e(3) = (eps(1,2)**2 + (eps(2,2) - Npr2 - Npl2)*(eps(1,1) - Npl2))/den + e(1) = sqrt(1/(1 + abs(e(3))**2 + abs(e(2))**2)) + e(2) = e(2)*e(1) + e(3) = e(3)*e(1) + else + if (sox < 0) then + e(3) = 1 + else + e(1) = sqrt(one/(1 + abs(-eps(1,1)/eps(1,2))**2)) + e(2) = -e(1)*eps(1,1)/eps(1,2) + end if + end if + end block polarization + +end subroutine warmdisp + + +subroutine dielectric_tensor(X, Y, mu, Npl, model, nlarmor, e330, epsl) + ! Returns the N⊥-independent parts of the dielectric tensor, namely + ! ε₃₃⁰ and the 3x3 tensors ε̅^(l) for l=1,nlarmor. + ! + ! Reference: https://doi.org/10.13182/FST08-A1660 + + implicit none + + ! subroutine arguments + + ! Inputs + + ! CMA diagram variables: X=(ω_pe/ω)², Y=ω_ce/ω + real(wp_), intent(in) :: X, Y + ! reciprocal of adimensional temperature: μ=mc²/kT + real(wp_), intent(in) :: mu + ! parallel refractive index: N∥ + real(wp_), intent(in) :: Npl + + ! Parameters + + ! switch for the dielectic tensor model: + ! 1: weakly relativistic + ! 2: fully relativistic (faster variant) + ! 3: fully relativistic (slower variant) + integer, intent(in) :: model + + ! order of the electron Larmor radius expansion + integer, intent(in) :: nlarmor + + ! Outputs + + ! element ε₃₃⁰, defined in eq. 9 + complex(wp_), intent(out) :: e330 + ! tensors ε̅_ij^(l), defined in eq. 10 + complex(wp_), intent(out) :: epsl(3,3,nlarmor) + + ! local variables + integer :: l, n + real(wp_) :: Npl2, alpha_nl, beta_nl, A_l + complex(wp_) :: Qp_0l, Qm_0l, Qp_1l, Qm_1l, Qp_2l + + ! weakly relativistic functions + complex(wp_), dimension(0:nlarmor,0:2) :: Fp, Fm + + ! fully relativistic functions + real(wp_) :: rr(-nlarmor:nlarmor, 0:2, 0:nlarmor) + real(wp_) :: ri(nlarmor, 0:2, nlarmor) + + Npl2 = Npl**2 + + if (model == 1) then + ! Compute the (combinations of) Shkarofsky functions F_q(n) + ! in which the functions Q⁺_hl(n), Q⁻_hl(n) can be expressed + call fsup(nlarmor, Y, Npl, mu, Fp, Fm) + else + ! Compute the real and imaginary parts of the integrals + ! defining the functions Q_hl(n), Q_hl(n) + block + real(wp_) :: cr, ci, cmxw + cmxw = 1 + 15/(8*mu) + 105/(128 * mu**2) + cr = -mu**2/(sqrt_pi * cmxw) + ci = sqrt(2*pi*mu) * mu**2/cmxw + + ! real part + if (model == 2) then + ! fast variant + call hermitian(rr, Y, mu, Npl, cr, model, nlarmor) + else + ! slow variant + call hermitian_2(rr, Y, mu, Npl, cr, model, nlarmor) + end if + ! imaginary part + call antihermitian(ri, Y, mu, Npl, ci, nlarmor) + end block end if - call antihermitian(ri,yg,mu,npl,ci,lrm) block - ! compute all factorials incrementally - integer :: fact_l, fact_2l, fact_l2 - integer :: fact_lpis, fact_lmis + ! Compute all factorials incrementally... + integer :: fact_l ! l! + integer :: fact_2l ! (2l)! + integer :: fact_2l_l ! (2l)! / l! + integer :: fact_l2 ! (l!)² + integer :: fact_lpn ! (l+n)! + integer :: fact_lmn ! (l-n)! - fact_l = 1 ! l! - fact_2l = 1 ! (2l)! - fact_l2 = 1 ! (l!)² + ! ...starting from 1 + fact_l = 1 + fact_2l = 1 + fact_2l_l = 1 + fact_l2 = 1 - do l=1,lrm - fact_l = fact_l * l - fact_2l = fact_2l * 2*l * (2*l - 1) - fact_l2 = fact_l2 * l**2 + epsl = 0 + do l = 1, nlarmor + fact_l = fact_l * l - lm=l-1 - fal=-0.25_wp_**l * fact_2l/(fact_l2 * yg**(2*lm)) - ca11=czero - ca12=czero - ca13=czero - ca22=czero - ca23=czero - ca33=czero + ! Coefficient A_l, eq. 11 + if (model == 1) then + fact_2l_l = fact_2l_l * (4*l - 2) ! see http://oeis.org/A001813 + A_l = (1/(2 * Y**2 * mu))**(l - 1)/2 * fact_2l_l + else + fact_2l = fact_2l * 2*l * (2*l - 1) + fact_l2 = fact_l2 * l**2 + A_l = - fact_2l/(4**l * Y**(2*(l - 1)) * fact_l2) + end if - fact_lpis = fact_l ! (l+is)! - fact_lmis = fact_l ! (l-is)! + ! When n=0, (l+n)! = (l-n)! = l! + fact_lpn = fact_l + fact_lmn = fact_l - do is=0,l - k=l-is - w=(-one)**k + do n = 0, l + ! Coefficients α_nl and β_nl, eq. 12 + alpha_nl = one * (-1)**(l - n)/(fact_lpn * fact_lmn) + beta_nl = alpha_nl*(n**2 + one * (2*(l - n)*(l + n)*(l - 1)) / (2*l - 1)) - asl=w/(fact_lpis * fact_lmis * one) - bsl=asl*(is*is+dble(2*k*lm*(l+is))/(2.0_wp_*l-one)) + ! Compute the next factorial values + fact_lpn = fact_lpn * (l + n+1) + fact_lmn = merge(fact_lmn / (l - n), 1, n < l) - fact_lpis = fact_lpis * (l + is+1) - fact_lmis = merge(fact_lmis / (l - is), 1, is < l) - - if(is.gt.0) then - cq0p=rr(is,0,l)+rr(-is,0,l)+im*ri(is,0,l) - cq0m=rr(is,0,l)-rr(-is,0,l)+im*ri(is,0,l) - cq1p=rr(is,1,l)+rr(-is,1,l)+im*ri(is,1,l) - cq1m=rr(is,1,l)-rr(-is,1,l)+im*ri(is,1,l) - cq2p=rr(is,2,l)+rr(-is,2,l)+im*ri(is,2,l) + ! Functions Q⁺_hl(n) Q⁻_hl(n), eq. 7 + if (model == 1) then + ! Weakly relativistic approximation using F_q(n) + Qp_0l = mu * Fp(n,0) + Qm_0l = mu * Fm(n,0) + Qp_1l = mu * Npl * (Fp(n,0) - Fp(n,1)) + Qm_1l = mu * Npl * (Fm(n,0) - Fm(n,1)) + Qp_2l = Fp(n,1) + mu * Npl2 * (Fp(n,2) + Fp(n,0) - 2*Fp(n,1)) else - cq0p=rr(is,0,l) - cq0m=rr(is,0,l) - cq1p=rr(is,1,l) - cq1m=rr(is,1,l) - cq2p=rr(is,2,l) + ! Fully relativistic computation + Qp_0l = rr(n,0,l) + merge(im*ri(n,0,l) + rr(-n,0,l), czero, n > 0) + Qm_0l = rr(n,0,l) + merge(im*ri(n,0,l) - rr(-n,0,l), czero, n > 0) + Qp_1l = rr(n,1,l) + merge(im*ri(n,1,l) + rr(-n,1,l), czero, n > 0) + Qm_1l = rr(n,1,l) + merge(im*ri(n,1,l) - rr(-n,1,l), czero, n > 0) + Qp_2l = rr(n,2,l) + merge(im*ri(n,2,l) + rr(-n,2,l), czero, n > 0) end if - ca11=ca11+is**2*asl*cq0p - ca12=ca12+is*l*asl*cq0m - ca22=ca22+bsl*cq0p - ca13=ca13+is*asl*cq1m/yg - ca23=ca23+l*asl*cq1p/yg - ca33=ca33+asl*cq2p/yg**2 + ! Components of the ε̅^(l) tensors, eq. 11 + epsl(1,1,l) = epsl(1,1,l) - n**2 * alpha_nl * Qp_0l + epsl(1,2,l) = epsl(1,2,l) + im * l * n*alpha_nl * Qm_0l + epsl(2,2,l) = epsl(2,2,l) - beta_nl * Qp_0l + epsl(1,3,l) = epsl(1,3,l) - n*alpha_nl * Qm_1l / Y + epsl(2,3,l) = epsl(2,3,l) - im * l*alpha_nl * Qp_1l / Y + epsl(3,3,l) = epsl(3,3,l) - alpha_nl * Qp_2l / Y**2 end do - epsl(1,1,l) = - ca11*fal - epsl(1,2,l) = + im*ca12*fal - epsl(2,2,l) = - ca22*fal - epsl(1,3,l) = - ca13*fal - epsl(2,3,l) = - im*ca23*fal - epsl(3,3,l) = - ca33*fal + + ! Add the common factor A_l + epsl(:,:,l) = epsl(:,:,l) * A_l end do end block - cq2p=rr(0,2,0) - a330=cq2p + ! Multiply by X and add the Kroneker δ_l,1 to diagonal elements + epsl = epsl * X + epsl(1,1,1) = epsl(1,1,1) + 1 + epsl(2,2,1) = epsl(2,2,1) + 1 - do l=1,lrm - epsl(2,1,l) = - epsl(1,2,l) - epsl(3,1,l) = epsl(1,3,l) - epsl(3,2,l) = - epsl(2,3,l) - end do + ! The special ε₃₃⁰ element + if (model == 1) then + Qp_2l = mu*Fp(0,1) + mu**2 * Npl2*(Fp(0,2) + Fp(0,0) - 2*Fp(0,1)) + else + Qp_2l = -rr(0,2,0) + end if + e330 = 1 - X * Qp_2l -end subroutine diel_tens_fr + ! Other components given by symmetry + epsl(2,1,:) = - epsl(1,2,:) + epsl(3,1,:) = + epsl(1,3,:) + epsl(3,2,:) = - epsl(2,3,:) + +end subroutine dielectric_tensor subroutine hermitian(rr,yg,mu,npl,cr,fast,lrm) use eierf, only : calcei3 implicit none -! arguments + + ! subroutine arguments integer :: lrm,fast real(wp_) :: yg,mu,npl,cr real(wp_), dimension(-lrm:lrm,0:2,0:lrm) :: rr -! local variables + + ! local variables integer :: i,k,n,n1,nn,m,llm real(wp_) :: mu2,mu4,mu6,npl2,npl4,bth,bth2,bth4,bth6,t,rxt,upl2, & upl,gx,exdx,x,gr,s,zm,zm2,zm3,fe0m,ffe,sy1,sy2,sy3 @@ -964,105 +1216,18 @@ subroutine expinit end subroutine expinit -subroutine diel_tens_wr(yg,mu,npl,a330,epsl,lrm) - ! Weakly relativistic dielectric tensor computation of dielectric - ! tensor elements (Krivenki and Orefice, JPP 30,125 - 1983) - ! - +pure subroutine fsup(lrm, yg, npl, mu, cefp, cefm) implicit none ! subroutine arguments - integer :: lrm - real(wp_) :: yg,npl,mu - complex(wp_) :: a330,epsl(3,3,lrm) + integer, intent(in) :: lrm + real(wp_), intent(in) :: yg, npl, mu + complex(wp_), intent(out), dimension(0:lrm,0:2) :: cefp, cefm + + ! local constants + real(wp_), parameter :: apsicr=0.7_wp_ ! local variables - integer :: l,lm,is,k - real(wp_) :: npl2,fcl,w,asl,bsl - complex(wp_) :: ca11,ca12,ca13,ca22,ca23,ca33,cq0p,cq0m,cq1p,cq1m,cq2p - complex(wp_), dimension(0:lrm,0:2) :: cefp,cefm - - npl2=npl*npl - - call fsup(cefp,cefm,lrm,yg,npl,mu) - - block - ! compute all factorials incrementally - integer :: fact_l, fact_2l_l - integer :: fact_lpis, fact_lmis - - fact_l = 1 ! l! - fact_2l_l = 1 ! (2l)!/l! - - do l=1,lrm - fact_l = fact_l * l - fact_2l_l = fact_2l_l * (4*l - 2) ! see http://oeis.org/A001813 - - lm=l-1 - fcl=0.5_wp_**l*((one/yg)**2/mu)**lm*fact_2l_l - ca11=czero - ca12=czero - ca13=czero - ca22=czero - ca23=czero - ca33=czero - - fact_lpis = fact_l ! (l+is)! - fact_lmis = fact_l ! (l-is)! - - do is=0,l - k=l-is - w=(-one)**k - - asl=w/(fact_lpis * fact_lmis * one) - bsl=asl*(is*is+dble(2*k*lm*(l+is))/(2.0_wp_*l-one)) - - fact_lpis = fact_lpis * (l + is+1) - fact_lmis = merge(fact_lmis / (l - is), 1, is < l) - - cq0p=mu*cefp(is,0) - cq0m=mu*cefm(is,0) - cq1p=mu*npl*(cefp(is,0)-cefp(is,1)) - cq1m=mu*npl*(cefm(is,0)-cefm(is,1)) - cq2p=cefp(is,1)+mu*npl2*(cefp(is,2)+cefp(is,0)-2.0_wp_*cefp(is,1)) - - ca11=ca11+is**2*asl*cq0p - ca12=ca12+is*l*asl*cq0m - ca22=ca22+bsl*cq0p - ca13=ca13+is*asl*cq1m/yg - ca23=ca23+l*asl*cq1p/yg - ca33=ca33+asl*cq2p/yg**2 - end do - epsl(1,1,l) = -ca11*fcl - epsl(1,2,l) = +im*ca12*fcl - epsl(2,2,l) = -ca22*fcl - epsl(1,3,l) = -ca13*fcl - epsl(2,3,l) = -im*ca23*fcl - epsl(3,3,l) = -ca33*fcl - end do - end block - - cq2p=cefp(0,1)+mu*npl2*(cefp(0,2)+cefp(0,0)-2.0_wp_*cefp(0,1)) - a330=-mu*cq2p - - do l=1,lrm - epsl(2,1,l) = - epsl(1,2,l) - epsl(3,1,l) = epsl(1,3,l) - epsl(3,2,l) = - epsl(2,3,l) - end do - -end subroutine diel_tens_wr - - -subroutine fsup(cefp,cefm,lrm,yg,npl,mu) - implicit none -! local constants - real(wp_), parameter :: apsicr=0.7_wp_ -! arguments - integer :: lrm - real(wp_) :: yg,npl,mu - complex(wp_), dimension(0:lrm,0:2) :: cefp,cefm -! local variables integer :: is,l,iq,ir,iflag real(wp_) :: psi,apsi,alpha,phi2,phim,xp,yp,xm,ym,x0,y0, & zrp,zip,zrm,zim,zr0,zi0 @@ -1208,7 +1373,7 @@ end subroutine fsup ! this work published in transactions on mathematical software, ! vol. 16, no. 1, pp. 47. ! - pure subroutine zetac (xi, yi, zr, zi, iflag) +pure subroutine zetac (xi, yi, zr, zi, iflag) ! ! given a complex number z = (xi,yi), this subroutine computes ! the value of the faddeeva-function w(z) = exp(-z**2)*erfc(-i*z), diff --git a/src/eccd.f90 b/src/eccd.f90 index 520eb7a..880a959 100644 --- a/src/eccd.f90 +++ b/src/eccd.f90 @@ -164,7 +164,7 @@ contains eccdpar(4+nlm:npar) = chlm end subroutine setcdcoeff_ncl - subroutine eccdeff(yg,anpl,anprre,dens,amu,ex,ey,ez,nhmn,nhmx, & + subroutine eccdeff(yg,anpl,anprre,dens,amu,e,nhmn,nhmx, & ithn,cst2,fcur,eccdpar,effjcd,iokhawa,ierr) use const_and_precisions, only : pi,qesi=>e_,mesi=>me_, & vcsi=>c_,qe=>ecgs_,me=>mecgs_,vc=>ccgs_,mc2=>mc2_ @@ -181,7 +181,7 @@ contains integer :: i,nhmn,nhmx,ithn,iokhawa,ierr real(wp_) :: yg,anpl,anprre,dens,amu,cst2,effjcd real(wp_), dimension(:) :: eccdpar - complex(wp_) :: ex,ey,ez + complex(wp_) :: e(3) ! local variables integer :: nhn,neval,ier,last,npar integer, dimension(liw) :: iw @@ -202,12 +202,12 @@ contains apar(2) = anpl apar(3) = amu apar(4) = anprre - apar(5) = dble(ex) - apar(6) = dimag(ex) - apar(7) = dble(ey) - apar(8) = dimag(ey) - apar(9) = dble(ez) - apar(10) = dimag(ez) + apar(5) = dble(e(1)) + apar(6) = dimag(e(1)) + apar(7) = dble(e(2)) + apar(8) = dimag(e(2)) + apar(9) = dble(e(3)) + apar(10) = dimag(e(3)) apar(11) = dble(ithn) npar=size(apar) diff --git a/src/gray_core.f90 b/src/gray_core.f90 index 320bd14..50da22b 100644 --- a/src/gray_core.f90 +++ b/src/gray_core.f90 @@ -44,7 +44,8 @@ contains real(wp_), parameter :: taucr = 12._wp_, etaucr = exp(-taucr) character, dimension(2), parameter :: mode=['O','X'] - real(wp_) :: sox,ak0,bres,xgcn,xg,yg,rrm,zzm,alpha,didp,anpl,anpr,anprim,anprre + integer :: sox + real(wp_) :: ak0,bres,xgcn,xg,yg,rrm,zzm,alpha,didp,anpl,anpr,anprim,anprre real(wp_) :: chipol,psipol,btot,psinv,dens,tekev,dersdst,derdnm real(wp_) :: tau,pow,dids,ddr,ddi,taumn,taumx real(wp_) :: rhotpav,drhotpav,rhotjava,drhotjava,dpdvp,jphip @@ -160,7 +161,7 @@ contains end if end if - sox=one ! mode inverted for each beam + sox=1 ! mode inverted for each beam iox=2 ! start with O: sox=-1, iox=1 psipol = params%antenna%psi @@ -425,9 +426,14 @@ contains if(ierrn==0 .and. params%ecrh_cd%iwarm>0 .and. ins_pl .and. & (tau1(jk)+tau0(jk)+lgcpl1(jk))<=taucr .and. .not.iwait(jk)) then ! H&CD computation check tekev=temp(psinv) - call alpha_effj(params%ecrh_cd, psinv, xg, yg, dens, tekev, & - ak0, bres, derdnm, anpl, anpr, sox, anprre, & - anprim, alpha, didp, nharm, nhf, iokhawa, ierrhcd) + block + complex(wp_) :: Npr_warm + call alpha_effj(params%ecrh_cd, psinv, xg, yg, dens, tekev, & + ak0, bres, derdnm, anpl, anpr, sox, Npr_warm, & + alpha, didp, nharm, nhf, iokhawa, ierrhcd) + anprre = Npr_warm%re + anprim = Npr_warm%im + end block if(ierrhcd/=0) then error = ior(error,ierrhcd) call print_errhcd(ierrhcd,i,anprre,anprim,alpha) @@ -965,12 +971,12 @@ contains ! use gray_params, only : igrad use beamdata, only : h,hh,h6 implicit none - real(wp_), intent(in) :: sox,bres,xgcn + real(wp_), intent(in) :: bres,xgcn real(wp_), dimension(6), intent(inout) :: y real(wp_), dimension(6), intent(in) :: yp real(wp_), dimension(3), intent(in) :: dgr real(wp_), dimension(3,3), intent(in) :: ddgr - integer, intent(in) :: igrad + integer, intent(in) :: igrad,sox real(wp_), dimension(6) :: yy,fk1,fk2,fk3,fk4 real(wp_) :: gr2 @@ -998,11 +1004,11 @@ contains implicit none ! arguments real(wp_), dimension(6), intent(in) :: y - real(wp_), intent(in) :: sox,bres,xgcn,gr2 + real(wp_), intent(in) :: bres,xgcn,gr2 real(wp_), dimension(3), intent(in) :: dgr2,dgr real(wp_), dimension(3,3), intent(in) :: ddgr real(wp_), dimension(6), intent(out) :: dery - integer, intent(in) :: igrad + integer, intent(in) :: igrad,sox ! local variables real(wp_) :: psinv,dens,btot,xg,yg real(wp_), dimension(3) :: xv,anv,bv,derxg,deryg @@ -1026,7 +1032,8 @@ contains real(wp_), dimension(3), intent(in) :: xv,anv real(wp_), dimension(3), intent(in) :: dgr real(wp_), dimension(3,3), intent(in) :: ddgr - real(wp_), intent(in) :: sox,bres,xgcn + integer, intent(in) :: sox + real(wp_), intent(in) :: bres,xgcn real(wp_), dimension(6), intent(out) :: dery real(wp_), intent(out) :: psinv,dens,btot,xg,yg,anpl,anpr real(wp_), intent(out) :: ddr,ddi,dersdst,derdnm @@ -1434,7 +1441,7 @@ contains ! refractive index N̅ vector, b̅ = B̅/B magnetic field direction real(wp_), dimension(3), intent(in) :: anv, bv ! sign of polarisation mode: -1 ⇒ O, +1 ⇒ X - real(wp_), intent(in) :: sox + integer, intent(in) :: sox ! CMA diagram variables: X=(ω_pe/ω)², Y=ω_ce/ω real(wp_), intent(in) :: xg, yg ! gradients of X, Y @@ -1690,12 +1697,12 @@ contains - subroutine alpha_effj(params, psinv, xg, yg, dens, tekev, ak0, bres, & - derdnm, anpl, anpr, sox, anprre, anprim, alpha, & - didp, nhmin, nhmax, iokhawa, error) + subroutine alpha_effj(params, psinv, X, Y, density, temperature, & + k0, Bres, derdnm, Npl, Npr_cold, sox, & + Npr, alpha, dIdp, nhmin, nhmax, iokhawa, error) ! Computes the absorption coefficient and effective current - use const_and_precisions, only : zero, pi, mc2=>mc2_ + use const_and_precisions, only : pi, mc2=>mc2_ use gray_params, only : ecrh_cd_parameters use coreprofiles, only : fzeff use equilibrium, only : sgnbphi @@ -1715,24 +1722,24 @@ contains ! poloidal flux ψ real(wp_), intent(in) :: psinv ! CMA diagram variables: X=(ω_pe/ω)², Y=ω_ce/ω - real(wp_), intent(in) :: xg, yg - ! density [10¹⁹ m⁻³], temperature [keV] - real(wp_), intent(in) :: dens, tekev + real(wp_), intent(in) :: X, Y + ! densityity [10¹⁹ m⁻³], temperature [keV] + real(wp_), intent(in) :: density, temperature ! vacuum wavenumber k₀=ω/c, resonant B field - real(wp_), intent(in) :: ak0, bres + real(wp_), intent(in) :: k0, Bres ! group velocity: |∂Λ/∂N̅| where Λ=c²/ω²⋅D_R real(wp_), intent(in) :: derdnm ! refractive index: N∥, N⊥ (cold dispersion) - real(wp_), intent(in) :: anpl, anpr + real(wp_), intent(in) :: Npl, Npr_cold ! sign of polarisation mode: -1 ⇒ O, +1 ⇒ X - real(wp_), intent(in) :: sox + integer, intent(in) :: sox ! Outputs - ! real/imaginary parts of N⊥ (warm dispersion) - real(wp_), intent(out) :: anprre, anprim + ! orthogonal refractive index N⊥ (solution of the warm dispersion) + complex(wp_), intent(out) :: Npr ! absorption coefficient, current density - real(wp_), intent(out) :: alpha, didp + real(wp_), intent(out) :: alpha, dIdp ! lowest/highest resonant harmonic numbers integer, intent(out) :: nhmin, nhmax ! ECCD/overall error codes @@ -1740,18 +1747,16 @@ contains ! local variables real(wp_) :: rbavi, rrii, rhop - integer :: lrm, ithn, ierrcd - real(wp_) :: amu, rbn, rbx + integer :: nlarmor, ithn, ierrcd + real(wp_) :: mu, rbn, rbx real(wp_) :: zeff, cst2, bmxi, bmni, fci real(wp_), dimension(:), pointer :: eccdpar=>null() - real(wp_) :: effjcd, effjcdav, btot - real(wp_) :: akim - complex(wp_) :: ex, ey, ez + real(wp_) :: effjcd, effjcdav, Btot + complex(wp_) :: e(3) - alpha = zero - anprim = zero - anprre = zero - didp = zero + alpha = 0 + Npr = 0 + dIdp = 0 nhmin = 0 nhmax = 0 iokhawa = 0 @@ -1760,37 +1765,52 @@ contains ! Absorption computation ! Skip when the temperature is zero (no plasma) - if (tekev <= zero) return + if (temperature <= 0) return ! Skip when the lowest resonant harmonic number is zero - amu = mc2/tekev ! μ=(mc²/kT) - call harmnumber(yg, amu, anpl, nhmin, nhmax, params%iwarm) + mu = mc2/temperature ! μ=(mc²/kT) + call harmnumber(Y, mu, Npl**2, params%iwarm == 1, nhmin, nhmax) if (nhmin <= 0) return - ! Solve the dispersion relation - lrm = max(params%ilarm, nhmax) - call warmdisp(xg, yg, amu, anpl, anpr, sox, lrm, error, & - anprre, anprim, params%iwarm, params%imx, & - ex, ey, ez) + ! Solve the full dispersion only when needed + if (params%iwarm /= 4 .or. params%ieccd /= 0) then + nlarmor = max(params%ilarm, nhmax) + if (params%ieccd /= 0) then + ! Compute the polarisation vector only when current drive is on + call warmdisp(X, Y, mu, Npl, Npr_cold, sox, error, Npr, e, & + model=params%iwarm, nlarmor=nlarmor, & + max_iters=abs(params%imx), fallback=params%imx < 0) + else + call warmdisp(X, Y, mu, Npl, Npr_cold, sox, error, Npr, & + model=params%iwarm, nlarmor=nlarmor, & + max_iters=abs(params%imx), fallback=params%imx < 0) + end if + end if + ! Compute α from the solution of the dispersion relation ! The absoption coefficient is defined as ! ! α = 2 Im(k̅)⋅s̅ ! ! where s̅ = v̅_g/|v_g|, the direction of the energy flow. - ! Since v̅_g = - ∂Λ/∂N̅ / ∂Λ/∂ω we have: + ! Since v̅_g = - ∂Λ/∂N̅ / ∂Λ/∂ω, using the cold dispersion + ! relation, we have that ! ! s̅ = ∂Λ/∂N̅ / |∂Λ/∂N̅| - ! = -2N̅ / |∂Λ/∂N̅| (using the cold dispersion) + ! = [2N̅ - ∂(N²s)/∂N∥ b̅ + ! + ½(b̅⋅∇S_I)² ∂³(N²s)/∂N∥³ b̅] / |∂Λ/∂N̅| ! ! Assuming Im(k∥)=0: ! ! α = 4 Im(k⊥)⋅N⊥ / |∂Λ/∂N̅| ! - akim = ak0 * anprim ! imaginary part of k⊥ = k₀N⊥ - alpha = 4 * akim * anpr/derdnm + block + real(wp_) :: k_im + k_im = k0 * Npr%im ! imaginary part of k⊥ + alpha = 4 * k_im*Npr_cold / derdnm + end block - if (alpha < zero) then + if (alpha < 0) then error = ibset(error, palph) return end if @@ -1800,36 +1820,36 @@ contains zeff = fzeff(psinv) ithn = 1 - if (lrm > nhmin) ithn = 2 + if (nlarmor > nhmin) ithn = 2 rhop = sqrt(psinv) call fluxval(rhop, rri=rrii, rbav=rbavi, bmn=bmni, bmx=bmxi, fc=fci) - btot = yg*bres - rbn = btot/bmni - rbx = btot/bmxi + Btot = Y*Bres + rbn = Btot/bmni + rbx = Btot/bmxi select case(params%ieccd) case(1) ! Cohen model - call setcdcoeff(zeff,rbn,rbx,cst2,eccdpar) - call eccdeff(yg,anpl,anprre,dens,amu,ex,ey,ez,nhmin,nhmax, & - ithn,cst2,fjch,eccdpar,effjcd,iokhawa,ierrcd) + call setcdcoeff(zeff, rbn, rbx, cst2, eccdpar) + call eccdeff(Y, Npl, Npr%re, density, mu, e, nhmin, nhmax, & + ithn, cst2, fjch, eccdpar, effjcd, iokhawa, ierrcd) case(2) ! No trapping - call setcdcoeff(zeff,cst2,eccdpar) - call eccdeff(yg,anpl,anprre,dens,amu,ex,ey,ez,nhmin,nhmax, & - ithn,cst2,fjch0,eccdpar,effjcd,iokhawa,ierrcd) + call setcdcoeff(zeff, cst2, eccdpar) + call eccdeff(Y, Npl, Npr%re, density, mu, e, nhmin, nhmax, & + ithn, cst2, fjch0, eccdpar, effjcd, iokhawa, ierrcd) case default ! Neoclassical model - call setcdcoeff(zeff,rbx,fci,amu,rhop,cst2,eccdpar) - call eccdeff(yg,anpl,anprre,dens,amu,ex,ey,ez,nhmin,nhmax, & - ithn,cst2,fjncl,eccdpar,effjcd,iokhawa,ierrcd) + call setcdcoeff(zeff, rbx, fci, mu, rhop, cst2, eccdpar) + call eccdeff(Y, Npl, Npr%re, density, mu, e, nhmin, nhmax, & + ithn, cst2, fjncl, eccdpar, effjcd, iokhawa, ierrcd) end select error = error + ierrcd if (associated(eccdpar)) deallocate(eccdpar) ! current drive efficiency [A⋅m/W] effjcdav = rbavi*effjcd - didp = sgnbphi * effjcdav / (2.0_wp_*pi*rrii) + dIdp = sgnbphi * effjcdav / (2*pi*rrii) end subroutine alpha_effj @@ -1846,8 +1866,8 @@ contains ! subroutine arguments real(wp_), dimension(6, nray), intent(in) :: ywrk0 - real(wp_), intent(in) :: bres, sox - integer, intent(in) :: ipol + real(wp_), intent(in) :: bres + integer, intent(in) :: sox, ipol real(wp_), intent(inout) :: psipol0, chipol0 complex(wp_), dimension(nray), intent(out) :: ext0, eyt0 diff --git a/src/multipass.f90 b/src/multipass.f90 index 45c069a..7df3349 100644 --- a/src/multipass.f90 +++ b/src/multipass.f90 @@ -19,7 +19,8 @@ contains ! arguments integer, intent(in) :: i ! ray index real(wp_), dimension(3), intent(in) :: xv,anv - real(wp_), intent(in) :: bres,sox + real(wp_), intent(in) :: bres + integer, intent(in) :: sox real(wp_), dimension(2), intent(out) :: cpl ! coupling (O/X) real(wp_), intent(out) :: psipol1,chipol1 integer, dimension(:), intent(inout), pointer :: iop ! in/out plasma flag @@ -67,7 +68,8 @@ contains ! arguments integer, intent(in) :: i ! ray index real(wp_), dimension(3), intent(in) :: xv,anv - real(wp_), intent(in) :: bres,sox + real(wp_), intent(in) :: bres + integer, intent(in) :: sox integer, dimension(:), intent(inout), pointer :: iop ! in/out plasma flag complex(wp_), dimension(:), intent(out), pointer :: ext,eyt ! local variables @@ -101,7 +103,8 @@ contains integer, intent(in) :: i ! ray index logical, intent(in) :: ins ! inside plasma? (ins=1 plasma/wall overlap) real(wp_), dimension(3), intent(inout) :: xv,anv - real(wp_), intent(in) :: bres,sox + real(wp_), intent(in) :: bres + integer, intent(in) :: sox real(wp_), intent(out) :: psipol1,chipol1 integer, dimension(:), intent(inout), pointer :: iow,iop ! in/out vessel and plasma flags complex(wp_), dimension(:), intent(inout), pointer :: ext,eyt diff --git a/src/polarization.f90 b/src/polarization.f90 index 3ceb059..5c28598 100644 --- a/src/polarization.f90 +++ b/src/polarization.f90 @@ -51,7 +51,8 @@ contains implicit none ! arguments real(wp_), dimension(3), intent(in) :: anv,bv - real(wp_), intent(in) :: bres,sox + real(wp_), intent(in) :: bres + integer, intent(in) :: sox complex(wp_), intent(out) :: ext,eyt ! real(wp_), optional, intent(out) :: gam ! local variables