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