From ae9f6051119f16422be00c36f400413127ce5c45 Mon Sep 17 00:00:00 2001 From: Michele Guerini Rocco Date: Sun, 17 Dec 2023 23:55:45 +0100 Subject: [PATCH] src/gray_core.f90: set initial ray position to z=0 MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit 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. --- src/gray_core.f90 | 14 +++++++++----- 1 file changed, 9 insertions(+), 5 deletions(-) diff --git a/src/gray_core.f90 b/src/gray_core.f90 index 8da0ded..c1edae8 100644 --- a/src/gray_core.f90 +++ b/src/gray_core.f90 @@ -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