diff --git a/src/gray.f b/src/gray.f
index 16e3505..f23c239 100644
--- a/src/gray.f
+++ b/src/gray.f
@@ -503,8 +503,8 @@ c
       
       write(4,*)' #sst R z phi psi rhot ne Te Br Bphi Bz Jphi '//
      .'N Npl ki alpha tau Pt Jcd dIds nh iohkw index_rt'
-      write(8,*) ' #istep j k xt yt zt rt psin'
-      write(9,*) ' #istep j k xt yt zt rt psin'
+      write(8,*) ' #istep j k xt yt zt rt psin modvg modn vgdotn'
+      write(9,*) ' #istep j k xt yt zt rt psin modvg modn vgdotn'
       write(17,*) ' #sst Dr_11 Dr_Nr1 Di_Nr1'
       write(33,*) ' #i j k sst x y R z psi tauv Npl alpha index_rt'
       write(12,*) ' #i sst psi w1 w2'
@@ -2855,6 +2855,7 @@ c
       dimension ywrk(ndim,jmx,kmx),ypwrk(ndim,jmx,kmx)
       dimension grad2(jmx,kmx),dgrad2v(3,jmx,kmx)
       dimension gri(3,jmx,kmx),ggri(3,3,jmx,kmx)
+      dimension vgv(3),vgv11(3)
 c
       common/nray/nrayr,nrayth
       common/dsds/dst
@@ -2867,6 +2868,9 @@ c
       common/gr/gr2
       common/dgr/dgr2,dgr,ddgr
       common/igrad/igrad
+      
+      
+      common/vgv11/vgv,vgv11
 c
       h=dst
       hh=h*0.5d0
@@ -2906,6 +2910,7 @@ c
      .        +h6*(fk1(ieq)+2.0d0*fk2(ieq)+2.0d0*fk3(ieq)+fk4(ieq))
           end do
         end do
+        if(j.eq.1) vgv11=vgv
       end do
 c
       call updatepos
@@ -2970,9 +2975,10 @@ c
       implicit real*8 (a-h,o-z)
       parameter(ndim=6)
       dimension y(ndim),dery(ndim)
-      dimension xv(3),anv(3),vgv(3),bv(3),derbv(3,3),derxg(3),deryg(3)
+      dimension xv(3),anv(3),bv(3),derbv(3,3),derxg(3),deryg(3)
       dimension derdxv(3),danpldxv(3),derdnv(3)
       dimension dgr2(3),dgr(3),ddgr(3,3),dbgr(3)
+      dimension vgv(3),vgv11(3)
 c
       common/gr/gr2
       common/dgr/dgr2,dgr,ddgr
@@ -2991,6 +2997,7 @@ c
       common/anv/anv
       common/xv/xv
       common/idst/idst
+      common/vgv11/vgv,vgv11
 c
       xx=y(1)
       yy=y(2)
@@ -5655,12 +5662,14 @@ c     gg=F(u)/u with F(u) as in Cohen paper
       implicit real*8 (a-h,o-z)
       parameter(jmx=31,kmx=36)
       dimension ywrk(6,jmx,kmx),ypwrk(6,jmx,kmx)
+      dimension vgv(3),vgv11(3)
 c
       common/nray/nrayr,nrayth
       common/wrk/ywrk,ypwrk
       common/psinv11/psinv11
       common/istep/istep
       common/ss/st
+      common/vgv11/vgv,vgv11
 c
       rtimn=1.d+30
       rtimx=-1.d-30
@@ -5676,10 +5685,16 @@ c
           dy=ywrk(2,j,k)-ywrk(2,1,1)
           dz=ywrk(3,j,k)-ywrk(3,1,1)
 c
-          dirx=ywrk(4,j,k)
-          diry=ywrk(5,j,k)
-          dirz=ywrk(6,j,k)
+          dirxn=ywrk(4,j,k)
+          diryn=ywrk(5,j,k)
+          dirzn=ywrk(6,j,k)
+          dirn=sqrt(dirxn*dirxn+diryn*diryn+dirzn*dirzn)
+          dirx=vgv11(1)
+          diry=vgv11(2)
+          dirz=vgv11(3)
           dir=sqrt(dirx*dirx+diry*diry+dirz*dirz)
+c          
+          vgdotn=dirxn*dirx+diryn*diry+dirzn*dirz
 c
           if(j.eq.1.and.k.eq.1) then
             csth1=dirz/dir
@@ -5708,7 +5723,7 @@ c
           end if
 
           if(.not.(iproj.eq.0.and.j.eq.1))
-     .    write(nfile,111) istep,j,k,xti,yti,zti,rti,psinv11
+     .    write(nfile,111) istep,j,k,xti,yti,zti,rti,psinv11,dir,dirn,vgdotn
 c
           if(rti.ge.rtimx.and.j.eq.nrayr) rtimx=rti
           if(rti.le.rtimn.and.j.eq.nrayr) rtimn=rti