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

fix coupling for subsequent beams

Browse Source 
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.
This commit is contained in:
Michele Guerini Rocco 2024-03-02 17:32:03 +01:00
parent 38a8edd439
commit f82f91bc8d
Signed by: rnhmjoj
GPG Key ID: BFBAF4C975F76450
3 changed files with 117 additions and 202 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

147
src/gray_core.f90
View File

@ -49,7 +49,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
@ -59,8 +59,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
@ -154,21 +154,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(.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 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
! | | | | ! | | | |
@ -196,9 +190,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(index_rt,iroff,iboff,iwait,stv,jphi_beam, & call initbeam(index_rt,iroff,iboff,iwait,stv,jphi_beam, &
pins_beam,currins_beam,dpdv_beam,jcd_beam) pins_beam,currins_beam,dpdv_beam,jcd_beam)
@ -223,10 +217,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, stv, xc, du1, gri, ggri, index_rt) ! * initial conditions
stv, 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_
@ -305,14 +296,15 @@ 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=.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(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 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
@ -325,6 +317,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
@ -347,13 +342,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]
@ -395,7 +390,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)
@ -406,13 +402,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
@ -532,12 +528,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(&
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) 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
@ -566,6 +563,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
@ -1840,74 +1838,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
! 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
logical, intent(in) :: 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
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 else
! X/O mode are undefined if B=0 e_x = 1
psipol0 = 0 e_y = 0
chipol0 = 0
end if 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

102
src/multipass.f90
View File

@ -13,55 +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
! arguments use const_and_precisions, only: cm
! subroutine arguments
integer, intent(in) :: i ! ray index integer, intent(in) :: i ! ray index
real(wp_), dimension(3), intent(in) :: xv,anv real(wp_), intent(in) :: x(3), N(3) ! position, refactive index
real(wp_), intent(in) :: bres real(wp_), intent(in) :: Bres ! resonant B field
integer, intent(in) :: sox integer, intent(in) :: sox ! sign of polarisation mode: -1 ⇒ O, +1 ⇒ X
real(wp_), dimension(2), intent(out) :: cpl ! coupling (O/X) real(wp_), intent(out) :: cpl(2) ! power coupling vector (O, X)
real(wp_), intent(out) :: psipol1,chipol1 real(wp_), intent(out) :: psi, chi ! polarisation ellipse angles
integer, dimension(:), intent(inout), pointer :: iop ! in/out plasma flag integer, intent(inout), pointer :: iop(:) ! inside/outside plasma flag
complex(wp_), dimension(:), intent(inout), pointer :: ext,eyt complex(wp_), intent(inout), pointer :: ext(:), eyt(:) ! ray polarisation vector (e_x, e_y)
! local variables logical, intent(in) :: perfect ! whether to assume perfect coupling
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 ! local variables
rm=sqrt(xmv(1)**2+xmv(2)**2) real(wp_) :: R, z, cosphi, sinphi, B_phi, B_R, B_z
csphi=xmv(1)/rm real(wp_) :: B(3)
snphi=xmv(2)/rm real(wp_) :: c
call bfield(rm,xmv(3),br,bz,bphi) complex(wp_) :: e_mode(2), e_ray(2)
bv(1)=br*csphi-bphi*snphi
bv(2)=br*snphi+bphi*csphi
bv(3)=bz
call pol_limit(anv,bv,bres,sox,ext1,eyt1) ! Update the inside/outside flag
call stokes_ce(ext1,eyt1,qq1,uu1,vv1) ! stokes parameter at plasma entrance iop(i) = iop(i) + 1
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 ! Compute B in cartesian coordinates
cpl=(/powcpl,one-powcpl/) R = norm2(x(1:2)) * cm
else ! incoming mode = X z = x(3) * cm
cpl=(/one-powcpl,powcpl/) cosphi = x(1)/R * cm
end if 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
if(i.eq.1) then ! polarization angles at plasma entrace for central ray ! Get the polarisation vector of the given mode
call polellipse(qq1,uu1,vv1,psipol1,chipol1) call pol_limit(N, B, Bres, sox, e_mode(1), e_mode(2))
psipol1=psipol1/degree ! convert from rad to degree
chipol1=chipol1/degree 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 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
! arguments ! arguments
integer, intent(in) :: i ! ray index integer, intent(in) :: i ! ray index
@ -239,7 +257,7 @@ contains
p0ray(nray),taus(nray,nbeam_tot),tau1(nray),etau1(nray), & p0ray(nray),taus(nray,nbeam_tot),tau1(nray),etau1(nray), &
cpls(nray,nbeam_tot),cpl1(nray),lgcpl1(nray),jphi_beam(dim), & cpls(nray,nbeam_tot),cpl1(nray),lgcpl1(nray),jphi_beam(dim), &
pins_beam(dim),currins_beam(dim),dpdv_beam(dim),jcd_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 end subroutine alloc_multipass
! ------------------------------ ! ------------------------------

54
src/polarization.f90
View File

@ -97,58 +97,4 @@ 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
! 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