From 5966ee8e9b7c829e6d0c7dce5e5e392c0e3e347a Mon Sep 17 00:00:00 2001 From: Michele Guerini Rocco Date: Sat, 2 Mar 2024 17:32:03 +0100 Subject: [PATCH] fix incorrect coupling for X mode beams --- src/gray_core.f90 | 145 +++++++++++++++---------------------------- src/multipass.f90 | 116 +++++++++++++++++++--------------- src/polarization.f90 | 54 ---------------- 3 files changed, 116 insertions(+), 199 deletions(-) diff --git a/src/gray_core.f90 b/src/gray_core.f90 index 3ec99b0..cb93646 100644 --- a/src/gray_core.f90 +++ b/src/gray_core.f90 @@ -64,7 +64,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 @@ -74,8 +74,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 @@ -180,21 +180,15 @@ contains iroff,yynext,yypnext,yw0,ypw0, & 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 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 ! | | | | @@ -222,9 +216,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(params, index_rt, iroff, iboff, iwait, stv, dst, & jphi_beam, pins_beam, currins_beam, dpdv_beam, jcd_beam) @@ -249,10 +243,7 @@ contains lgcpl1 = zero p0ray = p0jk ! * initial beam power - call ic_gb(params, anv0, ak0, yw, ypw, & ! * initial conditions - 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, xc, du1, gri, ggri, index_rt) ! * initial conditions do jk=1,params%raytracing%nray zzm = yw(3,jk)*0.01_wp_ @@ -341,14 +332,13 @@ 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=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(params%raytracing%ipol == 0) 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 @@ -361,6 +351,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 @@ -383,13 +376,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] @@ -431,7 +424,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) @@ -442,13 +436,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 @@ -601,12 +595,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 @@ -635,6 +630,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 @@ -2149,68 +2145,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 - - implicit none + 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, 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(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 + 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 3cba3d5..b1c4bba 100644 --- a/src/multipass.f90 +++ b/src/multipass.f90 @@ -13,56 +13,73 @@ module multipass contains -! ------------------------------ - subroutine plasma_in(i,xv,anv,bres,sox,cpl,psipol1,chipol1,iop,ext,eyt) ! ray enters plasma - implicit none -! 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),bphi,br,bz) - 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_phi, B_R, B_z) + 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 implicit none ! arguments @@ -110,7 +127,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 ! @@ -293,8 +309,8 @@ contains allocate(dpdv_beam(dim)) allocate(jcd_beam(dim)) - allocate(psipv(nbeam_tot)) - allocate(chipv(nbeam_tot)) + allocate(psipv(0:nbeam_tot)) + allocate(chipv(0:nbeam_tot)) end subroutine alloc_multipass diff --git a/src/polarization.f90 b/src/polarization.f90 index 2ed4055..db46042 100644 --- a/src/polarization.f90 +++ b/src/polarization.f90 @@ -98,59 +98,5 @@ 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 - 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