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

fix incorrect coupling for X mode beams

Browse Source
This commit is contained in:
Michele Guerini Rocco 2024-03-02 17:32:03 +01:00
parent 7eeb34c7dc
commit 5966ee8e9b
Signed by: rnhmjoj
GPG Key ID: BFBAF4C975F76450
3 changed files with 116 additions and 199 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

145
src/gray_core.f90
View File

@ -64,7 +64,7 @@ contains
real(wp_) :: rhotp,drhotp,rhotj,drhotj,dpdvmx,jphimx,ratjamx,ratjbmx real(wp_) :: rhotp,drhotp,rhotj,drhotj,dpdvmx,jphimx,ratjamx,ratjbmx
real(wp_) :: pabs_beam,icd_beam,cpl_beam1,cpl_beam2,cpl_cbeam1,cpl_cbeam2 real(wp_) :: pabs_beam,icd_beam,cpl_beam1,cpl_beam2,cpl_cbeam1,cpl_cbeam2
real(wp_), dimension(2) :: pabs_pass,icd_pass,cpl,cpl0 real(wp_), dimension(2) :: pabs_pass,icd_pass,cpl
real(wp_), dimension(3) :: xv,anv0,anv,bv,derxg real(wp_), dimension(3) :: xv,anv0,anv,bv,derxg
! Ray variables ! Ray variables
@ -74,8 +74,8 @@ contains
! i: integration step, jk: global ray index ! i: integration step, jk: global ray index
integer :: i, jk integer :: i, jk
integer :: iox,nharm,nhf,nnd,iokhawa,ierrn,ierrhcd,index_rt integer :: iox,nharm,nhf,nnd,iokhawa,ierrn,ierrhcd, index_rt, parent_index_rt
integer :: ip,ib,iopmin,ipar,iO integer :: ip,ib,iopmin,ipar,child_index_rt
integer :: igrad_b,istop_pass,nbeam_pass,nlim integer :: igrad_b,istop_pass,nbeam_pass,nlim
logical :: ins_pl,ins_wl,ent_pl,ext_pl,ent_wl,ext_wl,iboff logical :: ins_pl,ins_wl,ent_pl,ext_pl,ent_wl,ext_wl,iboff
@ -180,21 +180,15 @@ contains
iroff,yynext,yypnext,yw0,ypw0, & iroff,yynext,yypnext,yw0,ypw0, &
stnext,p0ray,taus,tau1,etau1,cpls,cpl1,lgcpl1,psipv,chipv) stnext,p0ray,taus,tau1,etau1,cpls,cpl1,lgcpl1,psipv,chipv)
if(params%raytracing%ipol == 0) then
if(params%antenna%iox == 2) then ! only X mode on 1st pass
cpl0 = [zero, one]
else ! only O mode on 1st pass
cpl0 = [one, zero]
end if
end if
sox=1 ! mode inverted for each beam sox=1 ! mode inverted for each beam
iox=2 ! start with O: sox=-1, iox=1 iox=2 ! start with O: sox=-1, iox=1
psipol = params%antenna%psi
chipol = params%antenna%chi
call pweight(params%antenna%power, p0jk) call pweight(params%antenna%power, p0jk)
! Set original polarisation
psipv(0) = params%antenna%psi
chipv(0) = params%antenna%chi
nbeam_pass=1 ! max n of beam per pass nbeam_pass=1 ! max n of beam per pass
index_rt=0 ! global beam index: 1,O 2,X 1st pass index_rt=0 ! global beam index: 1,O 2,X 1st pass
! | | | | ! | | | |
@ -222,9 +216,9 @@ contains
sox = -sox ! invert mode sox = -sox ! invert mode
iox = 3-iox ! O-mode at ip=1,ib=1 iox = 3-iox ! O-mode at ip=1,ib=1
index_rt = index_rt +1 index_rt = index_rt +1
iO = 2*index_rt +1 ! * index_rt of O-mode derived ray (iX=iO+1) child_index_rt = 2*index_rt + 1 ! * index_rt of O-mode child beam
parent_index_rt = ceiling(index_rt / 2.0_wp_) - 1 ! * index_rt of parent beam
call initbeam(params, index_rt, iroff, iboff, iwait, stv, dst, & call initbeam(params, index_rt, iroff, iboff, iwait, stv, dst, &
jphi_beam, pins_beam, currins_beam, dpdv_beam, jcd_beam) jphi_beam, pins_beam, currins_beam, dpdv_beam, jcd_beam)
@ -249,10 +243,7 @@ contains
lgcpl1 = zero lgcpl1 = zero
p0ray = p0jk ! * initial beam power p0ray = p0jk ! * initial beam power
call ic_gb(params, anv0, ak0, yw, ypw, & ! * initial conditions call ic_gb(params, anv0, ak0, yw, ypw, xc, du1, gri, ggri, index_rt) ! * initial conditions
xc, du1, gri, ggri, index_rt)
call set_pol(yw, bres, sox, params%raytracing%ipol, & ! * initial polarization
psipol, chipol, ext, eyt)
do jk=1,params%raytracing%nray do jk=1,params%raytracing%nray
zzm = yw(3,jk)*0.01_wp_ zzm = yw(3,jk)*0.01_wp_
@ -341,14 +332,13 @@ 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 plasma_in(jk,xv,anv,bres,sox,cpl,psipol,chipol,iop,ext,eyt) call set_pol(psipv(parent_index_rt), chipv(parent_index_rt), ext, eyt) ! compute polarisation and couplings
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)
if(iop(jk) == 1 .and. ip==1) then ! * 1st entrance on 1st pass (ray hasn't entered in plasma yet) => continue current pass if(iop(jk) == 1 .and. ip==1) then ! * 1st entrance on 1st pass (ray hasn't entered in plasma yet) => continue current pass
if(params%raytracing%ipol == 0) then ! + IF single mode propagation if(cpl(iox) < etaucr) then ! + IF low coupled power for current mode => de-activate derived rays
cpl = cpl0
p0ray(jk) = p0ray(jk)*cpl(iox)
else if(cpl(iox) < etaucr) then ! + ELSE IF low coupled power for current mode => de-activate derived rays
call turnoffray(jk,ip+1,2*ib+2-iox,iroff) call turnoffray(jk,ip+1,2*ib+2-iox,iroff)
iwait(jk) = .true. ! . stop advancement and H&CD computation for current ray iwait(jk) = .true. ! . stop advancement and H&CD computation for current ray
if(cpl(iox).le.comp_tiny) cpl(iox)=etaucr if(cpl(iox).le.comp_tiny) cpl(iox)=etaucr
@ -361,6 +351,9 @@ contains
write (msg,'(" 1st pass - central ray (",a1,"-mode) c=",g0.4)') & write (msg,'(" 1st pass - central ray (",a1,"-mode) c=",g0.4)') &
mode(iox), cpl(iox) mode(iox), cpl(iox)
call log_info(msg, mod='gray_core', proc='gray_main') call log_info(msg, mod='gray_core', proc='gray_main')
write (msg,'(" polarisation: ψ=",g0.5,"°, χ=",g0.5,"°")') &
psipol, chipol
call log_debug(msg, mod='gray_core', proc='gray_main')
psipv(index_rt) = psipol ! + polarization angles at plasma boundary for central ray psipv(index_rt) = psipol ! + polarization angles at plasma boundary for central ray
chipv(index_rt) = chipol chipv(index_rt) = chipol
end if end if
@ -383,13 +376,13 @@ contains
if(cpl(2).le.comp_tiny) cpl(2)=etaucr if(cpl(2).le.comp_tiny) cpl(2)=etaucr
end if end if
taus(jk,iO:iO+1) = tau1(jk) + tau0(jk) ! . starting tau for next O-mode pass taus(jk,child_index_rt:child_index_rt+1) = tau1(jk) + tau0(jk) ! . starting tau for next O-mode pass
cpls(jk,iO) = cpl1(jk) * cpl(1) ! . cumulative coupling for next O-mode pass cpls(jk,child_index_rt) = cpl1(jk) * cpl(1) ! . cumulative coupling for next O-mode pass
cpls(jk,iO+1) = cpl1(jk) * cpl(2) ! . cumulative coupling for next X-mode pass cpls(jk,child_index_rt+1) = cpl1(jk) * cpl(2) ! . cumulative coupling for next X-mode pass
if(jk == 1) then ! . polarization angles at plasma boundary for central ray if(jk == 1) then ! . polarization angles at plasma boundary for central ray
psipv(iO:iO+1) = psipol psipv(child_index_rt:child_index_rt+1) = psipol
chipv(iO:iO+1) = chipol chipv(child_index_rt:child_index_rt+1) = chipol
end if end if
else ! * 1st entrance on 2nd+ pass (ray hasn't entered in plasma since end of previous pass) => continue current pass else ! * 1st entrance on 2nd+ pass (ray hasn't entered in plasma since end of previous pass) => continue current pass
cpl = [zero, zero] cpl = [zero, zero]
@ -431,7 +424,8 @@ contains
yypnext(:,jk,index_rt) = ypw(:,jk) ! for next pass = reflection point yypnext(:,jk,index_rt) = ypw(:,jk) ! for next pass = reflection point
stnext(jk,index_rt) = stv(jk) ! . starting step for next pass = step after reflection stnext(jk,index_rt) = stv(jk) ! . starting step for next pass = step after reflection
call plasma_in(jk,xv,anv,bres,sox,cpl,psipol,chipol,iop,ext,eyt) ! . ray re-enters plasma after reflection call plasma_in(jk, xv, anv, bres, sox, cpl, psipol, chipol, & ! . ray re-enters plasma after reflection
iop, ext, eyt, perfect=.false.)
if(cpl(1) < etaucr) then ! . low coupled power for O-mode? => de-activate derived rays if(cpl(1) < etaucr) then ! . low coupled power for O-mode? => de-activate derived rays
call turnoffray(jk,ip+1,2*ib-1,iroff) call turnoffray(jk,ip+1,2*ib-1,iroff)
@ -442,13 +436,13 @@ contains
if(cpl(2).le.comp_tiny) cpl(2)=etaucr if(cpl(2).le.comp_tiny) cpl(2)=etaucr
end if end if
taus(jk,iO:iO+1) = tau1(jk) + tau0(jk) ! . starting tau for next O-mode pass taus(jk,child_index_rt:child_index_rt+1) = tau1(jk) + tau0(jk) ! . starting tau for next O-mode pass
cpls(jk,iO) = cpl1(jk) * cpl(1) ! . cumulative coupling for next O-mode pass cpls(jk,child_index_rt) = cpl1(jk) * cpl(1) ! . cumulative coupling for next O-mode pass
cpls(jk,iO+1) = cpl1(jk) * cpl(2) ! . cumulative coupling for next X-mode pass cpls(jk,child_index_rt+1) = cpl1(jk) * cpl(2) ! . cumulative coupling for next X-mode pass
if(jk == 1) then ! + polarization angles at plasma boundary for central ray if(jk == 1) then ! + polarization angles at plasma boundary for central ray
psipv(iO:iO+1) = psipol psipv(child_index_rt:child_index_rt+1) = psipol
chipv(iO:iO+1) = chipol chipv(child_index_rt:child_index_rt+1) = chipol
end if end if
end if end if
end if end if
@ -601,12 +595,13 @@ contains
icd_pass(iox) = icd_pass(iox) + icd_beam icd_pass(iox) = icd_pass(iox) + icd_beam
if(ip < params%raytracing%ipass .and. iopmin > 2) then ! not last pass AND at least one ray re-entered plasma if(ip < params%raytracing%ipass .and. iopmin > 2) then ! not last pass AND at least one ray re-entered plasma
cpl_beam1 = sum(p0ray * exp(-tau0) * cpls(:,iO)/cpl1, MASK=iop > 2) / & cpl_beam1 = sum(&
sum(p0ray * exp(-tau0), MASK=iop > 2) ! * average O-mode coupling for next beam (on active rays) p0ray * exp(-tau0) * cpls(:,child_index_rt)/cpl1, MASK=iop > 2) / &
sum(p0ray * exp(-tau0), MASK=iop > 2) ! * average O-mode coupling for next beam (on active rays)
cpl_beam2 = one - cpl_beam1 ! * average X-mode coupling for next beam cpl_beam2 = one - cpl_beam1 ! * average X-mode coupling for next beam
if(iop(1) > 2) then ! * central ray O/X-mode coupling for next beam if(iop(1) > 2) then ! * central ray O/X-mode coupling for next beam
cpl_cbeam1 = cpls(1,iO)/cpl1(1) cpl_cbeam1 = cpls(1,child_index_rt)/cpl1(1)
cpl_cbeam2 = one - cpl_cbeam1 cpl_cbeam2 = one - cpl_cbeam1
end if end if
else ! last pass OR no ray re-entered plasma else ! last pass OR no ray re-entered plasma
@ -635,6 +630,7 @@ contains
if(iop(1) > 2) then if(iop(1) > 2) then
write(msg, '(3x,a,(g0.4,", ",g0.4))') & write(msg, '(3x,a,(g0.4,", ",g0.4))') &
'coupling [ctr ray, O/X]:', cpl_cbeam1, cpl_cbeam2 ! central ray coupling for next O/X beams 'coupling [ctr ray, O/X]:', cpl_cbeam1, cpl_cbeam2 ! central ray coupling for next O/X beams
call log_info(msg, mod='gray_core', proc='gray_main')
end if end if
end if end if
@ -2149,68 +2145,27 @@ contains
end subroutine alpha_effj end subroutine alpha_effj
subroutine set_pol(psi, chi, e_x, e_y)
subroutine set_pol(ywrk0, bres, sox, ipol, psipol0, chipol0, ext0, eyt0) ! Compute the polarisation vector for each ray from the ellipse angles ψ, χ
use const_and_precisions, only : degree, zero, one, half, im use const_and_precisions, only : degree, im
use beamdata, only : nray,nrayth use polarization, only : stokes_ell
use equilibrium, only : bfield
use polarization, only : pol_limit, polellipse, &
stokes_ce, stokes_ell
implicit none
! subroutine arguments ! subroutine arguments
real(wp_), dimension(6, nray), intent(in) :: ywrk0 real(wp_), intent(in) :: psi, chi
real(wp_), intent(in) :: bres complex(wp_), intent(out) :: e_x(:), e_y(:)
integer, intent(in) :: sox, ipol
real(wp_), intent(inout) :: psipol0, chipol0
complex(wp_), dimension(nray), intent(out) :: ext0, eyt0
! local variables ! local variables
integer :: j,k,jk real(wp_) :: qq, uu, vv, deltapol
real(wp_), dimension(3) :: xmv, anv, bv
real(wp_) :: rm, csphi, snphi, bphi, br, bz, qq, uu, vv, deltapol
j=1 call stokes_ell(chi*degree, psi*degree, qq, uu, vv)
k=0 if(qq**2 < 1) then
do jk=1,nray deltapol = asin(vv/sqrt(1 - qq**2))
k=k+1 e_x = sqrt((1 + qq)/2)
if(jk == 2 .or. k > nrayth) then e_y = sqrt((1 - qq)/2) * exp(-im * deltapol)
j=j+1 else
k=1 e_x = 1
end if e_y = 0
end if
if(ipol == 0) then
xmv=ywrk0(1:3,jk)*0.01_wp_ ! convert from cm to m
anv=ywrk0(4:6,jk)
rm=sqrt(xmv(1)**2+xmv(2)**2)
csphi=xmv(1)/rm
snphi=xmv(2)/rm
call bfield(rm,xmv(3),bphi,br,bz)
! bv(i) = B_i in cartesian coordinates
bv(1)=br*csphi-bphi*snphi
bv(2)=br*snphi+bphi*csphi
bv(3)=bz
call pol_limit(anv,bv,bres,sox,ext0(jk),eyt0(jk))
if (jk == 1) then
call stokes_ce(ext0(jk),eyt0(jk),qq,uu,vv)
call polellipse(qq,uu,vv,psipol0,chipol0)
psipol0=psipol0/degree ! convert from rad to degree
chipol0=chipol0/degree
end if
else
call stokes_ell(chipol0*degree,psipol0*degree,qq,uu,vv)
if(qq**2 < one) then
deltapol=asin(vv/sqrt(one - qq**2))
ext0(jk)= sqrt(half*(one + qq))
eyt0(jk)= sqrt(half*(one - qq))*exp(-im*deltapol)
else
ext0(jk)= one
eyt0(jk)= zero
end if
endif
end do
end subroutine set_pol end subroutine set_pol

116
src/multipass.f90
View File

@ -13,56 +13,73 @@ module multipass
contains contains
! ------------------------------ subroutine plasma_in(i, x, N, Bres, sox, cpl, psi, chi, iop, ext, eyt, perfect)
subroutine plasma_in(i,xv,anv,bres,sox,cpl,psipol1,chipol1,iop,ext,eyt) ! ray enters plasma ! Computes the ray polarisation and power couplings when it enteres the plasma
implicit none use const_and_precisions, only: cm
! arguments
integer, intent(in) :: i ! ray index ! subroutine arguments
real(wp_), dimension(3), intent(in) :: xv,anv integer, intent(in) :: i ! ray index
real(wp_), intent(in) :: bres real(wp_), intent(in) :: x(3), N(3) ! position, refactive index
integer, intent(in) :: sox real(wp_), intent(in) :: Bres ! resonant B field
real(wp_), dimension(2), intent(out) :: cpl ! coupling (O/X) integer, intent(in) :: sox ! sign of polarisation mode: -1 ⇒ O, +1 ⇒ X
real(wp_), intent(out) :: psipol1,chipol1 real(wp_), intent(out) :: cpl(2) ! power coupling vector (O, X)
integer, dimension(:), intent(inout), pointer :: iop ! in/out plasma flag real(wp_), intent(out) :: psi, chi ! polarisation ellipse angles
complex(wp_), dimension(:), intent(inout), pointer :: ext,eyt integer, intent(inout), pointer :: iop(:) ! inside/outside plasma flag
! local variables complex(wp_), intent(inout), pointer :: ext(:), eyt(:) ! ray polarisation vector (e_x, e_y)
real(wp_) :: rm,csphi,snphi,bphi,br,bz logical, intent(in) :: perfect ! whether to assume perfect coupling
real(wp_) :: qq1,uu1,vv1,qq,uu,vv,powcpl
real(wp_), dimension(3) :: bv,xmv ! local variables
complex(wp_) :: ext1,eyt1 real(wp_) :: R, z, cosphi, sinphi, B_phi, B_R, B_z
! real(wp_) :: B(3)
iop(i)=iop(i)+1 ! out->in real(wp_) :: c
complex(wp_) :: e_mode(2), e_ray(2)
xmv=xv*0.01_wp_ ! convert from cm to m
rm=sqrt(xmv(1)**2+xmv(2)**2) ! Update the inside/outside flag
csphi=xmv(1)/rm iop(i) = iop(i) + 1
snphi=xmv(2)/rm
call bfield(rm,xmv(3),bphi,br,bz) ! Compute B in cartesian coordinates
bv(1)=br*csphi-bphi*snphi R = norm2(x(1:2)) * cm
bv(2)=br*snphi+bphi*csphi z = x(3) * cm
bv(3)=bz cosphi = x(1)/R * cm
sinphi = x(2)/R * cm
call pol_limit(anv,bv,bres,sox,ext1,eyt1) call bfield(R, z, B_phi, B_R, B_z)
call stokes_ce(ext1,eyt1,qq1,uu1,vv1) ! stokes parameter at plasma entrance B(1) = B_R*cosphi - B_phi*sinphi
call stokes_ce(ext(i),eyt(i),qq,uu,vv) ! stokes parameter at plasma exit/wall reflection B(2) = B_R*sinphi + B_phi*cosphi
powcpl = half*(one + vv*vv1+uu*uu1+qq*qq1) ! coupling for incoming mode B(3) = B_z
if(sox.eq.-one) then ! incoming mode = O ! Get the polarisation vector of the given mode
cpl=(/powcpl,one-powcpl/) call pol_limit(N, B, Bres, sox, e_mode(1), e_mode(2))
else ! incoming mode = X
cpl=(/one-powcpl,powcpl/) if(i == 1) then
end if ! For the central ray, compute the polarization ellipse
chi = asin(2*imag(e_mode(1) * conjg(e_mode(2)))) / 2
if(i.eq.1) then ! polarization angles at plasma entrace for central ray psi = atan2(real(2 * e_mode(1) * conjg(e_mode(2))), &
call polellipse(qq1,uu1,vv1,psipol1,chipol1) real(dot_product(e_mode, conjg(e_mode)))) / 2
psipol1=psipol1/degree ! convert from rad to degree ! and convert from radians to degree
chipol1=chipol1/degree psi = psi / degree
chi = chi / degree
else else
psipol1 = zero psi = 0
chipol1 = zero chi = 0
end if end if
if (perfect) then
! Ignore the given vector and use the expected one
! Note: this will give 100% coupling to the current mode
ext(i) = e_mode(1)
eyt(i) = e_mode(2)
end if
! Compute the power coupling with the current mode
e_ray = [ext(i), eyt(i)]
c = abs(dot_product(e_mode, e_ray))**2
! Store both O and X couplings, in this order
c = merge(c, 1-c, sox == -1)
cpl = [c, 1-c]
end subroutine plasma_in end subroutine plasma_in
! ------------------------------
subroutine plasma_out(i,xv,anv,bres,sox,iop,ext,eyt) ! ray exits plasma subroutine plasma_out(i,xv,anv,bres,sox,iop,ext,eyt) ! ray exits plasma
implicit none implicit none
! arguments ! arguments
@ -110,7 +127,6 @@ contains
complex(wp_), dimension(:), intent(inout), pointer :: ext,eyt complex(wp_), dimension(:), intent(inout), pointer :: ext,eyt
! local variables ! local variables
integer :: irfl integer :: irfl
real(wp_) :: qq1,uu1,vv1
real(wp_), dimension(3) :: xvrfl,anvrfl,walln real(wp_), dimension(3) :: xvrfl,anvrfl,walln
complex(wp_) :: ext1,eyt1 complex(wp_) :: ext1,eyt1
! !
@ -293,8 +309,8 @@ contains
allocate(dpdv_beam(dim)) allocate(dpdv_beam(dim))
allocate(jcd_beam(dim)) allocate(jcd_beam(dim))
allocate(psipv(nbeam_tot)) allocate(psipv(0:nbeam_tot))
allocate(chipv(nbeam_tot)) allocate(chipv(0:nbeam_tot))
end subroutine alloc_multipass end subroutine alloc_multipass

54
src/polarization.f90
View File

@ -98,59 +98,5 @@ contains
! gam = atan2(sngam,csgam)/degree ! gam = atan2(sngam,csgam)/degree
end subroutine pol_limit end subroutine pol_limit
subroutine polarcold(anpl,anpr,xg,yg,sox,exf,eyif,ezf,elf,etf)
use const_and_precisions, only : wp_,zero,one
implicit none
! arguments
real(wp_), intent(in) :: anpl,anpr,xg,yg,sox
real(wp_), intent(out) :: exf,eyif,ezf,elf,etf
! local variables
real(wp_) :: anpl2,anpr2,an2,yg2,dy2,aa,e3,qq,p
if(xg <= zero) then
exf = zero
if(sox < zero) then
ezf = one
eyif = zero
else
ezf = zero
eyif = one
end if
elf = zero
etf = one
else
anpl2 = anpl**2
anpr2 = anpr**2
an2 = anpl2 + anpr2
yg2=yg**2
aa=1.0_wp_-xg-yg2
dy2 = one - yg2
qq = xg*yg/(an2*dy2 - aa)
if (anpl == zero) then
if(sox < zero) then
exf = zero
eyif = zero
ezf = one
else
qq = -aa/(xg*yg)
exf = one/sqrt(one + qq**2)
eyif = qq*exf
ezf = zero
end if
else
e3 = one - xg
p = (anpr2 - e3)/(anpl*anpr) ! undef for anpr==0
exf = p*ezf
eyif = qq*exf
ezf = one/sqrt(one + p**2*(one + qq**2))
end if
elf = (anpl*ezf + anpr*exf)/sqrt(an2)
etf = sqrt(one - elf**2)
end if
end subroutine polarcold
end module polarization end module polarization