istp/gray
2
0
Fork 0
Code Issues Pull Requests Packages Projects Releases 1 Wiki Activity

src/polarization.f90: rewrite

Browse Source 
- Remove the Stokes parameters as an intermediate step in the
  conversion between Jones vectors and polarisation ellipses.

- Document every single step performed when converting between
  different parametrisations and how the polarisation at the
  plasma boundary is computed. This includes how everything
  was derived from first principles.

- Mark the subroutines as pure.

- Remove `set_pol` entirely.
This commit is contained in:
Michele Guerini Rocco 2024-03-12 16:42:19 +01:00
parent f82f91bc8d
commit 0a87a3ef76
Signed by: rnhmjoj
GPG Key ID: BFBAF4C975F76450
3 changed files with 277 additions and 116 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

28
src/gray_core.f90
View File

@ -11,6 +11,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, EQ_VACUUM
use beams, only : xgygcoeff, launchangles2n
@ -296,7 +297,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=.not. params%raytracing%ipol &
.and. params%antenna%iox == iox &
@ -1838,30 +1840,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)

17
src/multipass.f90
View File

@ -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
@ -124,7 +119,6 @@ contains
complex(wp_), dimension(:), intent(inout), pointer :: ext,eyt
! local variables
integer :: irfl
real(wp_) :: qq1,uu1,vv1
real(wp_), dimension(3) :: xvrfl,anvrfl,walln
complex(wp_) :: ext1,eyt1
!
@ -136,11 +130,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

348
src/polarization.f90
View File

@ -1,100 +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
use const_and_precisions, only : wp_, pi, im
implicit none
interface stokes
module procedure stokes_ce,stokes_ell
end interface
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_
! 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
! 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
! 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
! 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