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

src/gray_core.f90: set initial ray position to z=0

Browse Source 
Previously the ray initial positions were set to the wavefront
S_R(x,y,z) = 0, with (x,y) chosen such that S_I(x,y,0) = const.
The wavefront itself, however, was determined using the value of the
beam parameters (k_ξ, k_η, w_ξ, w_η, etc.) fixed at z=0, which is valid
only when the initial wavefront is approximately flat.
Moreover, since the ray are distributed according to S_I(z=0), this
choice creates an inconsistency between the phase (from S_R at z≠0) and
the power (from S_I at z=0) assigned to the rays.

Satisfying both conditions on S_R and S_I exactly is really hard;
however, given we do not really care about the phase and we want to
precisely track the power, so it's more sensible to simply set z=0.

This means that when integrating in the phase (idst=2), gray will no
longer construct wavefronts, but merely transport the initial phase
on the z=0 plane.
Note that k₀⋅s will still give the correct phase. so if necessary, a
wavefront could be reconstructed by interpolating the (s, x̅, y̅, z̅)
points.
This commit is contained in:
Michele Guerini Rocco 2023-12-17 23:55:45 +01:00
parent 9a0f19e9e1
commit ae9f605111
Signed by: rnhmjoj
GPG Key ID: BFBAF4C975F76450
1 changed files with 9 additions and 5 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

14
src/gray_core.f90
View File

@ -223,7 +223,7 @@ contains
p0ray = p0jk ! * initial beam power
call ic_gb(params, anv0, ak0, yw, ypw, & ! * initial conditions
xc, du1, gri, ggri, index_rt)
stv, xc, du1, gri, ggri, index_rt)
call set_pol(yw, bres, sox, params%raytracing%ipol, & ! * initial polarization
psipol, chipol, ext, eyt)
@ -649,7 +649,7 @@ contains
subroutine ic_gb(params, anv0c, ak0, ywrk0, ypwrk0, &
xc0, du10, gri, ggri, index_rt)
stv, xc0, du10, gri, ggri, index_rt)
! beam tracing initial conditions igrad=1
! !!!!!! check ray tracing initial conditions igrad=0 !!!!!!
use const_and_precisions, only : zero,one,pi,half,two,degree,ui=>im
@ -663,6 +663,7 @@ contains
real(wp_), dimension(3), intent(in) :: anv0c
real(wp_), intent(in) :: ak0
real(wp_), dimension(6, nray), intent(out) :: ywrk0, ypwrk0
real(wp_), dimension(nrayth * nrayr), intent(out) :: stv
real(wp_), dimension(3, nrayth, nrayr), intent(out) :: xc0, du10
real(wp_), dimension(3, nray), intent(out) :: gri
real(wp_), dimension(3, 3, nray), intent(out) :: ggri
@ -848,7 +849,7 @@ contains
dy0t = dcsiw*snphiw + detaw*csphiw
x0t = uj(j)*dx0t
y0t = uj(j)*dy0t
z0t = -(half*(rcixx*x0t**2 + rciyy*y0t**2) + rcixy*x0t*y0t)
z0t = 0
dx0 = x0t*csps + snps*(y0t*csth + z0t*snth)
dy0 = -x0t*snps + csps*(y0t*csth + z0t*snth)
@ -940,12 +941,15 @@ contains
select case(params%raytracing%idst)
case(0)
! optical path: s
stv(jk) = 0
denom = an0
case(1)
! "time": c⋅t
stv(jk) = 0
denom = one
case(2)
! real eikonal: S_R
stv(jk) = half*(rcixx*x0t**2 + rciyy*y0t**2) + rcixy*x0t*y0t
denom = an20
end select
ypwrk0(1,jk) = anx/denom
@ -957,8 +961,8 @@ contains
ddr = anx**2 + any**2 + anz**2 - an20
ddi = 2*(anxt*gxt + anyt*gyt + anzt*gzt)
! st=0, index_rt=1, B=0, N=0, psin=-1, Xg=0, Yg=0, gradXg=0
call print_output(0,jk,zero,one,xc0(:,k,j),-one,zero,[zero,zero,zero], &
! i=0, index_rt=1, B=0, N=0, psin=-1, Xg=0, Yg=0, gradXg=0
call print_output(0,jk,stv(jk),one,xc0(:,k,j),-one,zero,[zero,zero,zero], &
ak0,zero,zero,[zero,zero,zero],zero,zero,zero,zero,zero,zero, &
0,0,0,index_rt,ddr,ddi,zero,zero,[zero,zero,zero])
end do