From 39b019db1d99f16b91fb02aad17d1c496353140e Mon Sep 17 00:00:00 2001 From: Michele Guerini Rocco Date: Tue, 12 Mar 2024 16:42:19 +0100 Subject: [PATCH] src/polarization.f90: rewrite --- src/gray_core.f90 | 28 +--- src/multipass.f90 | 16 +- src/polarization.f90 | 356 +++++++++++++++++++++++++++++++++---------- 3 files changed, 280 insertions(+), 120 deletions(-) diff --git a/src/gray_core.f90 b/src/gray_core.f90 index cb93646..f4fb4fb 100644 --- a/src/gray_core.f90 +++ b/src/gray_core.f90 @@ -22,6 +22,7 @@ contains use const_and_precisions, only : zero, one, degree, comp_tiny use coreprofiles, only : temp, fzeff use dispersion, only : expinit + use polarization, only : ellipse_to_field use gray_params, only : gray_parameters, gray_data, gray_results, & print_parameters use beams, only : xgygcoeff, launchangles2n @@ -332,7 +333,8 @@ contains if(ent_pl) then ! ray enters plasma write (msg, '(" ray ",g0," entered plasma (",g0," steps)")') jk, i call log_debug(msg, mod='gray_core', proc='gray_main') - call set_pol(psipv(parent_index_rt), chipv(parent_index_rt), ext, eyt) ! compute polarisation and couplings + call ellipse_to_field(psipv(parent_index_rt), chipv(parent_index_rt), & ! compute polarisation and couplings + ext, eyt) call plasma_in(jk, xv, anv, bres, sox, cpl, psipol, chipol, iop, ext, eyt, & perfect=params%raytracing%ipol == 0 .and. params%antenna%iox == iox .and. ip == 1) @@ -2145,30 +2147,6 @@ contains end subroutine alpha_effj - subroutine set_pol(psi, chi, e_x, e_y) - ! Compute the polarisation vector for each ray from the ellipse angles ψ, χ - use const_and_precisions, only : degree, im - use polarization, only : stokes_ell - - ! subroutine arguments - real(wp_), intent(in) :: psi, chi - complex(wp_), intent(out) :: e_x(:), e_y(:) - - ! local variables - real(wp_) :: qq, uu, vv, deltapol - - call stokes_ell(chi*degree, psi*degree, qq, uu, vv) - if(qq**2 < 1) then - deltapol = asin(vv/sqrt(1 - qq**2)) - e_x = sqrt((1 + qq)/2) - e_y = sqrt((1 - qq)/2) * exp(-im * deltapol) - else - e_x = 1 - e_y = 0 - end if - end subroutine set_pol - - subroutine cniteq(rqgrid,zqgrid,matr2dgrid,nr,nz,h,ncon,npts,icount,rcon,zcon) ! v2.01 12/07/95 -- written by d v bartlett, jet joint undertaking. ! (based on an older code) diff --git a/src/multipass.f90 b/src/multipass.f90 index b1c4bba..6efc7d2 100644 --- a/src/multipass.f90 +++ b/src/multipass.f90 @@ -2,7 +2,7 @@ module multipass use const_and_precisions, only : wp_, zero, half, one, degree, czero use beamdata, only : dst, nray use gray_params, only : ipass - use polarization, only : pol_limit, stokes_ce, polellipse + use polarization, only : pol_limit, field_to_ellipse use reflections, only : wall_refl use equilibrium, only : bfield @@ -52,12 +52,7 @@ contains if(i == 1) then ! For the central ray, compute the polarization ellipse - chi = asin(2*imag(e_mode(1) * conjg(e_mode(2)))) / 2 - psi = atan2(real(2 * e_mode(1) * conjg(e_mode(2))), & - real(dot_product(e_mode, conjg(e_mode)))) / 2 - ! and convert from radians to degree - psi = psi / degree - chi = chi / degree + call field_to_ellipse(e_mode(1), e_mode(2), psi, chi) else psi = 0 chi = 0 @@ -138,11 +133,8 @@ contains xv = xvrfl anv = anvrfl - if(i.eq.1) then ! polarization angles at wall reflection for central ray - call stokes_ce(ext1,eyt1,qq1,uu1,vv1) - call polellipse(qq1,uu1,vv1,psipol1,chipol1) - psipol1=psipol1/degree ! convert from rad to degree - chipol1=chipol1/degree + if(i == 1) then ! polarization angles at wall reflection for central ray + call field_to_ellipse(ext1, eyt1, psipol1, chipol1) else psipol1 = zero chipol1 = zero diff --git a/src/polarization.f90 b/src/polarization.f90 index db46042..30d0a28 100644 --- a/src/polarization.f90 +++ b/src/polarization.f90 @@ -1,102 +1,292 @@ +! This module contains subroutines to convert between different descriptions +! of the wave polarisation and to compute the polarisation of the plasma modes module polarization - interface stokes - module procedure stokes_ce,stokes_ell - end interface + + use const_and_precisions, only : wp_, pi, im + + implicit none + + private + public ellipse_to_field, field_to_ellipse ! Converting between descriptions + public pol_limit ! Plasma modes polarisations contains - subroutine stokes_ce(ext,eyt,qq,uu,vv) - use const_and_precisions, only : wp_ - implicit none -! arguments - complex(wp_), intent(in) :: ext,eyt - real(wp_), intent(out) :: qq,uu,vv - qq = abs(ext)**2 - abs(eyt)**2 - uu = 2*real(ext*conjg(eyt)) - vv = 2*imag(ext*conjg(eyt)) - end subroutine stokes_ce + pure subroutine ellipse_to_field(psi, chi, e_x, e_y) + ! Computes the normalised Jones vector from the + ! polarisation ellipse angles ψ, χ + ! + ! Notes: + ! - ψ∈[-π/2, π/2] is the angle between the x and the major axis + ! - χ∈[-π/4, π/4] is defined by tan(χ) = b/a, where a,b are the + ! major,minor semi-axis, respectively; χ>0 for positive helicity + ! (left-handed wave), χ<0 for negative helicity (right-handed wave). + + ! subroutine arguments + real(wp_), intent(in) :: psi, chi + complex(wp_), intent(out) :: e_x(:), e_y(:) + + ! The Eikonal ansatz is: + ! + ! E̅(r̅, t) = Re e̅(r̅) exp(-ik₀S(r̅) + iωt) + ! + ! where e̅(r̅) = [|e₁|exp(iφ₁), |e₂|exp(iφ₂), 0], since the wave + ! is transversal in vacuum. At a fixed position r̅=0, ignoring + ! the third component, we have: + ! + ! E̅(0, t) = [|e₁|cos(φ₁ + ωt), |e₂|cos(φ₂ + ωt)] + ! = [|e₁|cos(φ₁)cos(ωt) - |e₁|sin(φ₁)sin(ωt), + ! |e₂|cos(φ₂)cos(ωt) - |e₂|sin(φ₂)sin(ωt)] + ! + ! Then, we compare this to the parametric equation of + ! an ellipse rotated by ψ through the origin, + ! + ! P̅(t) = R(ψ) [acos(ωt), bsin(ωt)] + ! = [cos(ψ)a⋅cos(ωt), -sin(ψ)b⋅sin(ωt), + ! sin(ψ)a⋅cos(ωt), cos(ψ)b⋅sin(ωt)] + ! + ! at ωt=0 and ωt=π/2, so: + ! + ! 1. |e₁|cos(φ₁) = a⋅cos(ψ) + ! 2. |e₂|cos(φ₂) = a⋅sin(ψ) + ! 3. |e₁|sin(φ₁) = b⋅sin(ψ) + ! 4. |e₂|sin(φ₂) = -b⋅cos(ψ) + ! + ! From 1²+3² and 2²+4² we have + ! + ! |e₁|² = a²cos(ψ)²+b²sin(ψ)² + ! |e₂|² = a²sin(ψ)²+b²cos(ψ)² + ! + ! ⇒ |e₁|² + |e₂|² = a² + b² + ! + ! Assuming e̅ is normalised, that is e̅⋅e̅*=|e₁|²+|e₂|²=1, + ! we have a²+b²=1, so we can define a=cos(χ), b=sin(χ). + ! + ! We can then rewrite: + ! + ! 1. Re e₁ = cos(χ)cos(ψ) + ! 2. Re e₂ = cos(χ)sin(ψ) + ! 3. Im e₁ = +sin(χ)sin(ψ) + ! 4. Im e₂ = -sin(χ)cos(ψ) + ! + ! from which: + ! + ! e₁ = cos(χ)cos(ψ) + i⋅sin(χ)sin(ψ) + ! e₂ = cos(χ)sin(ψ) - i⋅sin(χ)cos(ψ) + ! + e_x = cosd(chi)*cosd(psi) + im * sind(chi)*sind(psi) + e_y = cosd(chi)*sind(psi) - im * sind(chi)*cosd(psi) + end subroutine ellipse_to_field - subroutine stokes_ell(chi,psi,qq,uu,vv) - use const_and_precisions, only : wp_,two - implicit none -! arguments - real(wp_), intent(in) :: chi,psi - real(wp_), intent(out) :: qq,uu,vv + pure subroutine field_to_ellipse(e_x, e_y, psi, chi) + ! Computes the polarisation ellipse angles ψ, χ + ! from the normalised Jones vector - qq=cos(two*chi)*cos(two*psi) - uu=cos(two*chi)*sin(two*psi) - vv=sin(two*chi) - end subroutine stokes_ell + ! subroutine arguments + complex(wp_), intent(in) :: e_x, e_y + real(wp_), intent(out) :: psi, chi + ! First, to make the map (ψ, χ) → unit ellipse + ! unique we restricted its domain to + ! [-π/2, π/2]×[-π/4, π/4]. Then, starting from + ! + ! e₁ = cos(χ)cos(ψ) + i⋅sin(χ)sin(ψ) + ! e₂ = cos(χ)sin(ψ) - i⋅sin(χ)cos(ψ) + ! + ! (see `ellipse_to_field`), we obtain χ doing + ! + ! 2e₁e₂* = cos(2χ)sin(2ψ) + i sin(2χ) + ! ⇒ χ = ½ asin(Im(2e₁e₂*)) + ! + ! This gives the angle χ in the correct interval since + ! the range of asin is [-π/2, π/2] by definition. + ! + ! For ψ notice that: + ! + ! e⋅e = e₁²-e₂² = cos(χ)²cos(2ψ) - sin(χ)²cos(2ψ) + ! = cos(2χ)cos(2ψ) + ! + ! Combining this to the previous expression: + ! + ! tan(2ψ) = Re(2e₁e₂*)/(e⋅e) + ! + ! Finally, since atan2 gives the principal branch + ! (-π, π], ψ is given within the correct interval as: + ! + ! ψ = ½ atan2(Re(2e₁e₂*), (e⋅e)) + ! + chi = asind(imag(2 * e_x * conjg(e_y))) / 2 + psi = atan2d(real(2 * e_x * conjg(e_y)), abs(e_x)**2 - abs(e_y)**2) / 2 + end subroutine field_to_ellipse - subroutine polellipse(qq,uu,vv,psi,chi) - use const_and_precisions, only : wp_,half - implicit none -! arguments - real(wp_), intent(in) :: qq,uu,vv - real(wp_), intent(out) :: psi,chi -! real(wp_) :: ll,aa,bb,ell + pure subroutine pol_limit(N, B, Bres, sox, e_x, e_y) + ! Computes the Jones vectors of the cold plasma dispersion + ! relation in the limit of vanishing electron density + ! + ! Note: the Jones vectors are given in the local beam frame, + ! that is, the z axis is aligned with the wave vector and x axis + ! lies in the tokamak equatorial plane. + ! This allows to directly compare the beam polarisation with + ! the plasma modes Jones vectors to obtain the power couplings. -! ll = sqrt(qq**2 + uu**2) -! aa = sqrt(half*(1 + ll)) -! bb = sqrt(half*(1 - ll)) -! ell = bb/aa - psi = half*atan2(uu,qq) - chi = half*asin(vv) - end subroutine polellipse + ! subroutine arguments + real(wp_), intent(in) :: N(3) ! N̅ refractive index + real(wp_), intent(in) :: B(3) ! B̅ magnetic field + real(wp_), intent(in) :: Bres ! resonant magnetic field + ! sign of polarisation mode: -1 ⇒ O, +1 ⇒ X + integer, intent(in) :: sox + ! Components of the Jones vector + complex(wp_), intent(out) :: e_x, e_y - subroutine pol_limit(anv,bv,bres,sox,ext,eyt) !,gam) - use const_and_precisions, only : wp_,ui=>im,zero,one - implicit none -! arguments - real(wp_), dimension(3), intent(in) :: anv,bv - real(wp_), intent(in) :: bres - integer, intent(in) :: sox - complex(wp_), intent(out) :: ext,eyt -! real(wp_), optional, intent(out) :: gam -! local variables - real(wp_), dimension(3) :: bnv - real(wp_) :: anx,any,anz,an2,an,anpl2,anpl,anpr,anxy, & - btot,yg,den,dnl,del0,ff,ff2,sngam,csgam -! - btot = sqrt(bv(1)**2+bv(2)**2+bv(3)**2) - bnv = bv/btot - yg = btot/bres + ! local variables + real(wp_) :: Y, Npl, delta + real(wp_) :: z(3), R(2,2), gamma + complex(wp_) :: f, e(2) - anx = anv(1) - any = anv(2) - anz = anv(3) - an2 = anx**2 + any**2 + anz**2 - an = sqrt(an2) - anxy = sqrt(anx**2 + any**2) + Y = norm2(B) / Bres ! Y = ω_ce/ω + Npl = dot_product(N, B) / norm2(B) ! N∥ = N̅⋅B̅/B - anpl = (anv(1)*bnv(1) + anv(2)*bnv(2) + anv(3)*bnv(3)) - anpl2= anpl**2 - anpr = sqrt(an2 - anpl2) + ! Coordinate systems in use + ! 1. plasma basis, definition of dielectric tensor + ! x̅ = N̅⊥/N⊥ + ! y̅ = B̅/B×N̅/N + ! z̅ = B̅/B + ! 2. polarisation basis, definition of Ẽ₁/Ẽ₂ + ! x̅ = -N̅/N×(B̅/B×N̅/N)= -B̅⊥/B⊥ + ! y̅ = B̅/B×N̅/N + ! z̅ = N̅/N + ! 3. beam basis, definition of rays initial conditions + ! x̅ = N̅/N×e̅₃ + ! y̅ = N̅/N×(N̅/N×e̅₃) + ! z̅ = N̅/N - dnl = one - anpl2 - del0 = sqrt(dnl**2 + 4.0_wp_*anpl2/yg**2) - - sngam = (anz*anpl - an2*bnv(3))/(an*anxy*anpr) - csgam = -(any*bnv(1) - anx*bnv(2))/ (anxy*anpr) - - ff = 0.5_wp_*yg*(dnl - sox*del0) - ff2 = ff**2 - den = ff2 + anpl2 - if (den>zero) then - ext = (ff*csgam - ui*anpl*sngam)/sqrt(den) - eyt = (-ff*sngam - ui*anpl*csgam)/sqrt(den) - den = sqrt(abs(ext)**2+abs(eyt)**2) - ext = ext/den - eyt = eyt/den - else ! only for XM (sox=+1) when N//=0 - ext = -ui*sngam - eyt = -ui*csgam + ! The cold plasma dispersion relation tensor in the plasma + ! basis (the usual one) is: + ! + ! [ (R+L)/2 - N∥² -i(R-L)/2 N∥N⊥ ] + ! Λ = [i(R-L)/2 (R+L)/2 - N² 0 ] + ! [N∥N⊥ 0 P - N⊥²] + ! + ! where P = 1 - X + ! L = 1 - X/(1+Y) + ! R = 1 - X/(1-Y) + ! X = (ω_p/ω)² + ! Y = ω_c/ω + ! + ! To compute the polarisation (the ratio Ẽ₁/Ẽ₂ in the wave + ! transverse plane) we switch to the polarisation basis. + ! The relations between the new and old unit vectors are: + ! + ! y̅' = y̅ + ! z̅' = N̅/N = (N̅∥ + N̅⊥)/N = (N∥/N)z̅ + (N⊥/N)x̅ + ! x̅' = y̅'×z̅' = y̅×z̅' = (N∥/N)x̅ - (N⊥/N)z̅ + ! + ! so the change of basis matrix is: + ! + ! [ N∥/N 0 N⊥/N] + ! M = [ 0 1 0 ] + ! [-N⊥/N 0 N∥/N] + ! + ! which is equivalent to a clockwise rotation of the + ! x̅∧z̅ plane around y̅ by the pitch angle θ: N̅∥ = Ncosθ. + ! + ! Consequently the tensor Λ in this basis is given by: + ! + ! Λ' = M⁻¹ Λ M + ! + ! Substituting the solution of the dispersion relation + ! N⊥=N⊥(X, Y, N∥, ±) into Λ' and taking lim X→0 yields: + ! + ! [ ½XY²(N∥²-1∓Δ) -iXYN∥ XY²N∥√(1-N∥²)] + ! Λ' = [ iXYN∥ ½XY²(1-N∥²∓Δ) iXY√(1-N∥²)] + o(X²) + ! [XY²N∥√(1-N∥²) -iXY√(1-N∥²) -XY²N∥²+X+Y²-1] + ! + ! where Δ = √[4N∥²/Y² + (1 - N∥²)²] + ! + ! At exactly X=0, Λ' reduces to diag(0,0,1), so for Λ'E=0 we must + ! have E₃=0, meaning we can effectively ignore the third column. + ! + ! The polarisation is determined by the ratio + ! + ! Ẽ₁/Ẽ₂ = -Λ₁₂/Λ₁₁ = -i Y[N∥² -1 ∓ Δ(X=0)]/(2N∥) + ! = i [Y(N∥²-1) ± Δ']/(2N∥) + ! ≡ i f± + ! + ! where Δ' = YΔ(X=0) = √[4N∥² + Y²(1 - N∥²)²] + ! + if (sox == -1 .and. Npl == 0) then + ! If Ẽ₂=0, the direction of the electric field is e̅ = (1, 0). + ! This happens for O mode at N∥=0, the polarisation here is + ! linear with Ẽ ∥ x̅ ∥ B̅. + ! This case is handled separately to avoid a division by zero. + e = [1, 0] + else + ! In general, if Ẽ₂≠0, the direction of the electric field is + ! + ! (Ẽ₁, Ẽ₂) ~ (Ẽ₁/Ẽ₂, 1) ~ (iẼ₁/Ẽ₂, i) = (f±, i) + ! + ! so, the polarisation is elliptical and the Jones vector is: + ! + ! e̅ = (f±, i) / √(1 + f±²) + ! + delta = sqrt(4*Npl**2 + Y**2*(1 - Npl**2)**2) + f = (Y*(Npl**2 - 1) + sox * delta) / (2*Npl) + e = [f, im] / sqrt(1 + f**2) end if -! gam = atan2(sngam,csgam)/degree + ! To obtain the polarisation in the beam basis we must compute the + ! angle γ such that a rotation of γ around z̅=N̅/N aligns the x̅ axis + ! of the wave transverse plane to the equatorial plane of the + ! tokamak, that is: + ! + ! R(γ,z̅)(x̅) ⋅ e̅₃ = 0 + ! + ! where R(γ,n̅) is given by Rodrigues' rotation formula; + ! x̅ = -B̅⊥/B⊥ (x axis in polarisation basis); + ! z̅ = N̅/N is (z axis in polarization basis) + ! e̅₃ is the vertical direction (standard basis). + ! + ! Since the rotation is a linear operation, + ! + ! R(γ,z̅)(-B̅⊥/B⊥) ⋅ e̅₃ = 0 ⇒ R(γ,z̅)B̅⊥ ⋅ e̅₃ = 0 + ! + ! and the condition simplifies to: + ! + ! R(γ,z̅)B̅⊥ ⋅ e̅₃ + ! = [B̅⊥cosγ + z̅×B̅⊥sinγ + z̅(z̅⋅B̅⊥)(1-cosγ)] ⋅ e̅₃ + ! = (B̅⊥cosγ + z̅×B̅⊥sinγ) ⋅ e̅₃ + ! = B⊥₃cosγ + (z̅×B̅⊥)₃sinγ + ! = 0 + ! + ! By definition B̅⊥ = z̅×(B̅×z̅)=B̅-(z̅⋅B̅)z̅, so: + ! + ! B⊥₃ = B̅₃ - (z̅⋅B̅)z̅₃ + ! + ! (z̅×B̅⊥) = [z̅×(z̅×(B̅×z̅))]₃ + ! = [-(z̅×(B̅×z̅))×z̅]₃ + ! = -(B̅×z̅)₃ + ! = (z₁B₂-z₂B₁) + ! + ! Substituting, we obtain: + ! + ! [B̅₃ - (z̅⋅B̅)z̅₃]cosγ + (z₁B₂-z₂B₁)sinγ = 0 + ! + ! ⇒ γ = atan2((z̅⋅B̅)z̅₃-B̅₃, z₁B₂-z₂B₁) + ! + z = N / norm2(N) + gamma = atan2(dot_product(B, z) * z(3) - B(3), z(1)*B(2) - z(2)*B(1)) + + ! Apply the rotation R(γ,z̅) to the Jones vector + ! + ! Note: Jones vectors transform as normal vectors under rotations + ! in physical space (SO(3) group); but as spinors on the Poincaré + ! sphere (SU(2) group). For example, a rotation by γ around z̅ + ! changes the longitude φ=2ψ by 2γ. + R = reshape([cos(gamma), -sin(gamma), sin(gamma), cos(gamma)], [2, 2]) + e = matmul(R, e) + e_x = e(1) + e_y = e(2) end subroutine pol_limit - end module polarization