From f82f91bc8dbf141c4310861af2f742813a9fd0a7 Mon Sep 17 00:00:00 2001 From: Michele Guerini Rocco Date: Sat, 2 Mar 2024 17:32:03 +0100 Subject: [PATCH] fix coupling for subsequent beams In situations when multiple beams are traced, either when allowing multiple plasma crossings (raytracing.ipass > 0) or the initial polarisation is mixed (raytracing.ipol == .true.), the couplings of all but the first beam (with least index_rt) were invalid. The bug is due to the re-use of the psipol,chipol variables as the beams are traced sequentially over the beam_loop. For the first beam being traced the psipol,chipol are correctly initialised to the user-defined value and the resulting coupling is correct. However, in each subsequent beam the values were not set to those of the parent beam (or to the user-defined value in the case of the first X mode beam), but to those of the previous beams (current index_rt - 1). This change repurposes the psipv,chipv arrays to store the polarisation of the parent beams, including the initial user-defined value and makes plasma_in always use these to compute the coupling. In addition, in the case the polarisation is not immediately known (i.e. if raytracing.ipol == .false.), this change postpones the computation of the Jones vector (ext, eyt) from the launch point, if the magnetic equilibrium is available, to when the ray actually crosses the plasma boundary. The original code, besides being strictly incorrect, can lead to non-negligible alterations to the coupling. This change also mean: 1. most of the functionality of `set_pol` has been merged with `plasma_in` 2. the polarisation is undefined and the Jones vector is set to the placeholder value [1, 0] till `plasma_im` is called Finally, `polarcold` is removed because it's unused. --- src/gray_core.f90 | 153 +++++++++++++++---------------------------- src/multipass.f90 | 112 ++++++++++++++++++------------- src/polarization.f90 | 54 --------------- 3 files changed, 117 insertions(+), 202 deletions(-) diff --git a/src/gray_core.f90 b/src/gray_core.f90 index e0dbdaa..7b56a53 100644 --- a/src/gray_core.f90 +++ b/src/gray_core.f90 @@ -49,7 +49,7 @@ contains 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_), 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 ! Ray variables @@ -59,8 +59,8 @@ contains ! i: integration step, jk: global ray index integer :: i, jk - integer :: iox,nharm,nhf,nnd,iokhawa,ierrn,ierrhcd,index_rt - integer :: ip,ib,iopmin,ipar,iO + integer :: iox,nharm,nhf,nnd,iokhawa,ierrn,ierrhcd, index_rt, parent_index_rt + integer :: ip,ib,iopmin,ipar,child_index_rt integer :: igrad_b,istop_pass,nbeam_pass,nlim logical :: ins_pl,ins_wl,ent_pl,ext_pl,ent_wl,ext_wl,iboff @@ -154,21 +154,15 @@ contains iroff,yynext,yypnext,yw0,ypw0, & stnext,p0ray,taus,tau1,etau1,cpls,cpl1,lgcpl1,psipv,chipv) - if(.not. params%raytracing%ipol) 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 iox=2 ! start with O: sox=-1, iox=1 - psipol = params%antenna%psi - chipol = params%antenna%chi 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 index_rt=0 ! global beam index: 1,O 2,X 1st pass ! | | | | @@ -196,9 +190,9 @@ contains sox = -sox ! invert mode iox = 3-iox ! O-mode at ip=1,ib=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(index_rt,iroff,iboff,iwait,stv,jphi_beam, & pins_beam,currins_beam,dpdv_beam,jcd_beam) @@ -223,10 +217,7 @@ contains lgcpl1 = zero p0ray = p0jk ! * initial beam power - call ic_gb(params, anv0, ak0, yw, ypw, & ! * initial conditions - stv, xc, du1, gri, ggri, index_rt) - call set_pol(yw, bres, sox, params%raytracing%ipol, & ! * initial polarization - psipol, chipol, ext, eyt) + call ic_gb(params, anv0, ak0, yw, ypw, stv, xc, du1, gri, ggri, index_rt) ! * initial conditions do jk=1,params%raytracing%nray zzm = yw(3,jk)*0.01_wp_ @@ -305,14 +296,15 @@ 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 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=.not. params%raytracing%ipol & + .and. params%antenna%iox == iox & + .and. iop(jk) == 0 .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(.not. params%raytracing%ipol) then ! + IF single mode propagation - 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 + if(cpl(iox) < etaucr) then ! + IF low coupled power for current mode => de-activate derived rays call turnoffray(jk,ip+1,2*ib+2-iox,iroff) iwait(jk) = .true. ! . stop advancement and H&CD computation for current ray if(cpl(iox).le.comp_tiny) cpl(iox)=etaucr @@ -325,6 +317,9 @@ contains write (msg,'(" 1st pass - central ray (",a1,"-mode) c=",g0.4)') & mode(iox), cpl(iox) 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 chipv(index_rt) = chipol end if @@ -347,13 +342,13 @@ contains if(cpl(2).le.comp_tiny) cpl(2)=etaucr end if - taus(jk,iO:iO+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,iO+1) = cpl1(jk) * cpl(2) ! . cumulative coupling for next X-mode pass + taus(jk,child_index_rt:child_index_rt+1) = tau1(jk) + tau0(jk) ! . starting tau for next O-mode pass + cpls(jk,child_index_rt) = cpl1(jk) * cpl(1) ! . cumulative coupling for next O-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 - psipv(iO:iO+1) = psipol - chipv(iO:iO+1) = chipol + psipv(child_index_rt:child_index_rt+1) = psipol + chipv(child_index_rt:child_index_rt+1) = chipol end if else ! * 1st entrance on 2nd+ pass (ray hasn't entered in plasma since end of previous pass) => continue current pass cpl = [zero, zero] @@ -395,7 +390,8 @@ contains 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 - 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 call turnoffray(jk,ip+1,2*ib-1,iroff) @@ -406,13 +402,13 @@ contains if(cpl(2).le.comp_tiny) cpl(2)=etaucr end if - taus(jk,iO:iO+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,iO+1) = cpl1(jk) * cpl(2) ! . cumulative coupling for next X-mode pass + taus(jk,child_index_rt:child_index_rt+1) = tau1(jk) + tau0(jk) ! . starting tau for next O-mode pass + cpls(jk,child_index_rt) = cpl1(jk) * cpl(1) ! . cumulative coupling for next O-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 - psipv(iO:iO+1) = psipol - chipv(iO:iO+1) = chipol + psipv(child_index_rt:child_index_rt+1) = psipol + chipv(child_index_rt:child_index_rt+1) = chipol end if end if end if @@ -532,12 +528,13 @@ contains 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 - cpl_beam1 = sum(p0ray * exp(-tau0) * cpls(:,iO)/cpl1, MASK=iop > 2) / & - sum(p0ray * exp(-tau0), MASK=iop > 2) ! * average O-mode coupling for next beam (on active rays) + cpl_beam1 = sum(& + 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 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 end if else ! last pass OR no ray re-entered plasma @@ -566,6 +563,7 @@ contains if(iop(1) > 2) then 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 + call log_info(msg, mod='gray_core', proc='gray_main') end if end if @@ -1840,74 +1838,27 @@ contains end subroutine alpha_effj - - subroutine set_pol(ywrk0, bres, sox, ipol, psipol0, chipol0, ext0, eyt0) - use const_and_precisions, only : degree, zero, one, half, im - use beamdata, only : nray,nrayth - use equilibrium, only : bfield - use polarization, only : pol_limit, polellipse, & - stokes_ce, stokes_ell + 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_), dimension(6, nray), intent(in) :: ywrk0 - real(wp_), intent(in) :: bres - integer, intent(in) :: sox - logical, intent(in) :: ipol - real(wp_), intent(inout) :: psipol0, chipol0 - complex(wp_), dimension(nray), intent(out) :: ext0, eyt0 + real(wp_), intent(in) :: psi, chi + complex(wp_), intent(out) :: e_x(:), e_y(:) ! local variables - integer :: j,k,jk - real(wp_), dimension(3) :: xmv, anv, bv - real(wp_) :: rm, csphi, snphi, bphi, br, bz, qq, uu, vv, deltapol + real(wp_) :: qq, uu, vv, deltapol - j=1 - k=0 - do jk=1,nray - k=k+1 - if(jk == 2 .or. k > nrayth) then - j=j+1 - k=1 - end if - - if(.not. ipol) 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),br,bz,bphi) - ! bv(i) = B_i in cartesian coordinates - bv(1)=br*csphi-bphi*snphi - bv(2)=br*snphi+bphi*csphi - bv(3)=bz - - if (norm2(bv) > 0) then - 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 - ! X/O mode are undefined if B=0 - psipol0 = 0 - chipol0 = 0 - 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 + 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 diff --git a/src/multipass.f90 b/src/multipass.f90 index 8a8e056..6a2bb30 100644 --- a/src/multipass.f90 +++ b/src/multipass.f90 @@ -13,55 +13,73 @@ module multipass contains -! ------------------------------ - subroutine plasma_in(i,xv,anv,bres,sox,cpl,psipol1,chipol1,iop,ext,eyt) ! ray enters plasma -! arguments - integer, intent(in) :: i ! ray index - real(wp_), dimension(3), intent(in) :: xv,anv - 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 - complex(wp_), dimension(:), intent(inout), pointer :: ext,eyt -! local variables - real(wp_) :: rm,csphi,snphi,bphi,br,bz - real(wp_) :: qq1,uu1,vv1,qq,uu,vv,powcpl - real(wp_), dimension(3) :: bv,xmv - complex(wp_) :: ext1,eyt1 -! - iop(i)=iop(i)+1 ! out->in - - xmv=xv*0.01_wp_ ! convert from cm to m - rm=sqrt(xmv(1)**2+xmv(2)**2) - csphi=xmv(1)/rm - snphi=xmv(2)/rm - call bfield(rm,xmv(3),br,bz,bphi) - bv(1)=br*csphi-bphi*snphi - bv(2)=br*snphi+bphi*csphi - bv(3)=bz - - call pol_limit(anv,bv,bres,sox,ext1,eyt1) - call stokes_ce(ext1,eyt1,qq1,uu1,vv1) ! stokes parameter at plasma entrance - call stokes_ce(ext(i),eyt(i),qq,uu,vv) ! stokes parameter at plasma exit/wall reflection - powcpl = half*(one + vv*vv1+uu*uu1+qq*qq1) ! coupling for incoming mode - - if(sox.eq.-one) then ! incoming mode = O - cpl=(/powcpl,one-powcpl/) - else ! incoming mode = X - cpl=(/one-powcpl,powcpl/) - end if - - if(i.eq.1) then ! polarization angles at plasma entrace for central ray - call polellipse(qq1,uu1,vv1,psipol1,chipol1) - psipol1=psipol1/degree ! convert from rad to degree - chipol1=chipol1/degree + subroutine plasma_in(i, x, N, Bres, sox, cpl, psi, chi, iop, ext, eyt, perfect) + ! Computes the ray polarisation and power couplings when it enteres the plasma + use const_and_precisions, only: cm + + ! subroutine arguments + integer, intent(in) :: i ! ray index + real(wp_), intent(in) :: x(3), N(3) ! position, refactive index + real(wp_), intent(in) :: Bres ! resonant B field + integer, intent(in) :: sox ! sign of polarisation mode: -1 ⇒ O, +1 ⇒ X + real(wp_), intent(out) :: cpl(2) ! power coupling vector (O, X) + real(wp_), intent(out) :: psi, chi ! polarisation ellipse angles + integer, intent(inout), pointer :: iop(:) ! inside/outside plasma flag + complex(wp_), intent(inout), pointer :: ext(:), eyt(:) ! ray polarisation vector (e_x, e_y) + logical, intent(in) :: perfect ! whether to assume perfect coupling + + ! local variables + real(wp_) :: R, z, cosphi, sinphi, B_phi, B_R, B_z + real(wp_) :: B(3) + real(wp_) :: c + complex(wp_) :: e_mode(2), e_ray(2) + + ! Update the inside/outside flag + iop(i) = iop(i) + 1 + + ! Compute B in cartesian coordinates + R = norm2(x(1:2)) * cm + z = x(3) * cm + cosphi = x(1)/R * cm + sinphi = x(2)/R * cm + call bfield(R, z, B_R, B_z, B_phi) + B(1) = B_R*cosphi - B_phi*sinphi + B(2) = B_R*sinphi + B_phi*cosphi + B(3) = B_z + + ! Get the polarisation vector of the given mode + call pol_limit(N, B, Bres, sox, e_mode(1), e_mode(2)) + + 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 else - psipol1 = zero - chipol1 = zero + psi = 0 + chi = 0 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 -! ------------------------------ + + subroutine plasma_out(i,xv,anv,bres,sox,iop,ext,eyt) ! ray exits plasma ! arguments integer, intent(in) :: i ! ray index @@ -239,7 +257,7 @@ contains p0ray(nray),taus(nray,nbeam_tot),tau1(nray),etau1(nray), & cpls(nray,nbeam_tot),cpl1(nray),lgcpl1(nray),jphi_beam(dim), & pins_beam(dim),currins_beam(dim),dpdv_beam(dim),jcd_beam(dim), & - psipv(nbeam_tot),chipv(nbeam_tot)) + psipv(0:nbeam_tot),chipv(0:nbeam_tot)) end subroutine alloc_multipass ! ------------------------------ diff --git a/src/polarization.f90 b/src/polarization.f90 index a7adda4..f5e8b27 100644 --- a/src/polarization.f90 +++ b/src/polarization.f90 @@ -97,58 +97,4 @@ contains ! gam = atan2(sngam,csgam)/degree end subroutine pol_limit - subroutine polarcold(anpl,anpr,xg,yg,sox,exf,eyif,ezf,elf,etf) - use const_and_precisions, only : wp_,zero,one -! 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