src/polarization.f90: rewrite
This commit is contained in:
Michele Guerini Rocco
parent
5966ee8e9b
commit
39b019db1d
GPG Key ID:
BFBAF4C975F76450
@ -22,6 +22,7 @@ contains
|
|||||||
use const_and_precisions, only : zero, one, degree, comp_tiny
|
use const_and_precisions, only : zero, one, degree, comp_tiny
|
||||||
use coreprofiles, only : temp, fzeff
|
use coreprofiles, only : temp, fzeff
|
||||||
use dispersion, only : expinit
|
use dispersion, only : expinit
|
||||||
|
use polarization, only : ellipse_to_field
|
||||||
use gray_params, only : gray_parameters, gray_data, gray_results, &
|
use gray_params, only : gray_parameters, gray_data, gray_results, &
|
||||||
print_parameters
|
print_parameters
|
||||||
use beams, only : xgygcoeff, launchangles2n
|
use beams, only : xgygcoeff, launchangles2n
|
||||||
@ -332,7 +333,8 @@ contains
|
|||||||
if(ent_pl) then ! ray enters plasma
|
if(ent_pl) then ! ray enters plasma
|
||||||
write (msg, '(" ray ",g0," entered plasma (",g0," steps)")') jk, i
|
write (msg, '(" ray ",g0," entered plasma (",g0," steps)")') jk, i
|
||||||
call log_debug(msg, mod='gray_core', proc='gray_main')
|
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, &
|
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)
|
perfect=params%raytracing%ipol == 0 .and. params%antenna%iox == iox .and. ip == 1)
|
||||||
|
|
||||||
@ -2145,30 +2147,6 @@ contains
|
|||||||
end subroutine alpha_effj
|
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)
|
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.
|
! v2.01 12/07/95 -- written by d v bartlett, jet joint undertaking.
|
||||||
! (based on an older code)
|
! (based on an older code)
|
||||||
|
|||||||
@ -2,7 +2,7 @@ module multipass
|
|||||||
use const_and_precisions, only : wp_, zero, half, one, degree, czero
|
use const_and_precisions, only : wp_, zero, half, one, degree, czero
|
||||||
use beamdata, only : dst, nray
|
use beamdata, only : dst, nray
|
||||||
use gray_params, only : ipass
|
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 reflections, only : wall_refl
|
||||||
use equilibrium, only : bfield
|
use equilibrium, only : bfield
|
||||||
|
|
||||||
@ -52,12 +52,7 @@ contains
|
|||||||
|
|
||||||
if(i == 1) then
|
if(i == 1) then
|
||||||
! For the central ray, compute the polarization ellipse
|
! For the central ray, compute the polarization ellipse
|
||||||
chi = asin(2*imag(e_mode(1) * conjg(e_mode(2)))) / 2
|
call field_to_ellipse(e_mode(1), e_mode(2), psi, chi)
|
||||||
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
|
|
||||||
else
|
else
|
||||||
psi = 0
|
psi = 0
|
||||||
chi = 0
|
chi = 0
|
||||||
@ -138,11 +133,8 @@ contains
|
|||||||
xv = xvrfl
|
xv = xvrfl
|
||||||
anv = anvrfl
|
anv = anvrfl
|
||||||
|
|
||||||
if(i.eq.1) then ! polarization angles at wall reflection for central ray
|
if(i == 1) then ! polarization angles at wall reflection for central ray
|
||||||
call stokes_ce(ext1,eyt1,qq1,uu1,vv1)
|
call field_to_ellipse(ext1, eyt1, psipol1, chipol1)
|
||||||
call polellipse(qq1,uu1,vv1,psipol1,chipol1)
|
|
||||||
psipol1=psipol1/degree ! convert from rad to degree
|
|
||||||
chipol1=chipol1/degree
|
|
||||||
else
|
else
|
||||||
psipol1 = zero
|
psipol1 = zero
|
||||||
chipol1 = zero
|
chipol1 = zero
|
||||||
|
|||||||
@ -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
|
module polarization
|
||||||
interface stokes
|
|
||||||
module procedure stokes_ce,stokes_ell
|
use const_and_precisions, only : wp_, pi, im
|
||||||
end interface
|
|
||||||
|
implicit none
|
||||||
|
|
||||||
|
private
|
||||||
|
public ellipse_to_field, field_to_ellipse ! Converting between descriptions
|
||||||
|
public pol_limit ! Plasma modes polarisations
|
||||||
|
|
||||||
contains
|
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
|
pure subroutine ellipse_to_field(psi, chi, e_x, e_y)
|
||||||
uu = 2*real(ext*conjg(eyt))
|
! Computes the normalised Jones vector from the
|
||||||
vv = 2*imag(ext*conjg(eyt))
|
! polarisation ellipse angles ψ, χ
|
||||||
end subroutine stokes_ce
|
!
|
||||||
|
! 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)
|
pure subroutine field_to_ellipse(e_x, e_y, psi, chi)
|
||||||
use const_and_precisions, only : wp_,two
|
! Computes the polarisation ellipse angles ψ, χ
|
||||||
implicit none
|
! from the normalised Jones vector
|
||||||
! arguments
|
|
||||||
real(wp_), intent(in) :: chi,psi
|
|
||||||
real(wp_), intent(out) :: qq,uu,vv
|
|
||||||
|
|
||||||
qq=cos(two*chi)*cos(two*psi)
|
! subroutine arguments
|
||||||
uu=cos(two*chi)*sin(two*psi)
|
complex(wp_), intent(in) :: e_x, e_y
|
||||||
vv=sin(two*chi)
|
real(wp_), intent(out) :: psi, chi
|
||||||
end subroutine stokes_ell
|
|
||||||
|
|
||||||
|
! 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)
|
pure subroutine pol_limit(N, B, Bres, sox, e_x, e_y)
|
||||||
use const_and_precisions, only : wp_,half
|
! Computes the Jones vectors of the cold plasma dispersion
|
||||||
implicit none
|
! relation in the limit of vanishing electron density
|
||||||
! arguments
|
!
|
||||||
real(wp_), intent(in) :: qq,uu,vv
|
! Note: the Jones vectors are given in the local beam frame,
|
||||||
real(wp_), intent(out) :: psi,chi
|
! that is, the z axis is aligned with the wave vector and x axis
|
||||||
! real(wp_) :: ll,aa,bb,ell
|
! 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)
|
! subroutine arguments
|
||||||
! aa = sqrt(half*(1 + ll))
|
real(wp_), intent(in) :: N(3) ! N̅ refractive index
|
||||||
! bb = sqrt(half*(1 - ll))
|
real(wp_), intent(in) :: B(3) ! B̅ magnetic field
|
||||||
! ell = bb/aa
|
real(wp_), intent(in) :: Bres ! resonant magnetic field
|
||||||
psi = half*atan2(uu,qq)
|
! sign of polarisation mode: -1 ⇒ O, +1 ⇒ X
|
||||||
chi = half*asin(vv)
|
integer, intent(in) :: sox
|
||||||
end subroutine polellipse
|
! Components of the Jones vector
|
||||||
|
complex(wp_), intent(out) :: e_x, e_y
|
||||||
|
|
||||||
subroutine pol_limit(anv,bv,bres,sox,ext,eyt) !,gam)
|
! local variables
|
||||||
use const_and_precisions, only : wp_,ui=>im,zero,one
|
real(wp_) :: Y, Npl, delta
|
||||||
implicit none
|
real(wp_) :: z(3), R(2,2), gamma
|
||||||
! arguments
|
complex(wp_) :: f, e(2)
|
||||||
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
|
|
||||||
|
|
||||||
anx = anv(1)
|
Y = norm2(B) / Bres ! Y = ω_ce/ω
|
||||||
any = anv(2)
|
Npl = dot_product(N, B) / norm2(B) ! N∥ = N̅⋅B̅/B
|
||||||
anz = anv(3)
|
|
||||||
an2 = anx**2 + any**2 + anz**2
|
|
||||||
an = sqrt(an2)
|
|
||||||
anxy = sqrt(anx**2 + any**2)
|
|
||||||
|
|
||||||
anpl = (anv(1)*bnv(1) + anv(2)*bnv(2) + anv(3)*bnv(3))
|
! Coordinate systems in use
|
||||||
anpl2= anpl**2
|
! 1. plasma basis, definition of dielectric tensor
|
||||||
anpr = sqrt(an2 - anpl2)
|
! 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
|
! The cold plasma dispersion relation tensor in the plasma
|
||||||
del0 = sqrt(dnl**2 + 4.0_wp_*anpl2/yg**2)
|
! basis (the usual one) is:
|
||||||
|
!
|
||||||
sngam = (anz*anpl - an2*bnv(3))/(an*anxy*anpr)
|
! [ (R+L)/2 - N∥² -i(R-L)/2 N∥N⊥ ]
|
||||||
csgam = -(any*bnv(1) - anx*bnv(2))/ (anxy*anpr)
|
! Λ = [i(R-L)/2 (R+L)/2 - N² 0 ]
|
||||||
|
! [N∥N⊥ 0 P - N⊥²]
|
||||||
ff = 0.5_wp_*yg*(dnl - sox*del0)
|
!
|
||||||
ff2 = ff**2
|
! where P = 1 - X
|
||||||
den = ff2 + anpl2
|
! L = 1 - X/(1+Y)
|
||||||
if (den>zero) then
|
! R = 1 - X/(1-Y)
|
||||||
ext = (ff*csgam - ui*anpl*sngam)/sqrt(den)
|
! X = (ω_p/ω)²
|
||||||
eyt = (-ff*sngam - ui*anpl*csgam)/sqrt(den)
|
! Y = ω_c/ω
|
||||||
den = sqrt(abs(ext)**2+abs(eyt)**2)
|
!
|
||||||
ext = ext/den
|
! To compute the polarisation (the ratio Ẽ₁/Ẽ₂ in the wave
|
||||||
eyt = eyt/den
|
! transverse plane) we switch to the polarisation basis.
|
||||||
else ! only for XM (sox=+1) when N//=0
|
! The relations between the new and old unit vectors are:
|
||||||
ext = -ui*sngam
|
!
|
||||||
eyt = -ui*csgam
|
! 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
|
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 subroutine pol_limit
|
||||||
|
|
||||||
|
|
||||||
end module polarization
|
end module polarization
|
||||||
|
|||||||
Loading…
Reference in New Issue
Block a user