added the calculation of the reference ray arclength and some modifications in projxyzt subroutine

src/gray.f

@@ -72,9 +72,16 @@ c     second pass into plasma
       subroutine gray_integration
       implicit real*8 (a-h,o-z)
       
+      parameter(jmx=31,kmx=36)
+      dimension dery0(3)
+      dimension ywrk(6,jmx,kmx),ypwrk(6,jmx,kmx)
+      
       common/ss/st
       common/dsds/dst
       common/istep/istep
+      common/s11/s11
+      common/dery0/dery0,dery0mod
+      common/wrk/ywrk,ypwrk
       common/nstep/nstep
       common/istop/istop
       common/strfl11/strfl11
@@ -93,6 +100,14 @@ c       advance one step
 
 c       calculations after one step:
 
+        ds11=0
+        do j=1,3
+          ds11=ds11+dery0(j)*ywrk(3+j,1,1) 
+        enddo
+        ds11=ds11/dery0mod      
+        ds11=dst/ds11
+        s11=s11+ds11
+
         call after_onestep(istep,istop)
         if(istop.eq.1) exit
 c
@@ -2587,6 +2602,7 @@ c
       implicit real*8(a-h,o-z)
 c
       common/istep/istep
+      common/s11/s11
       common/istgr/istpr,istpl
       common/ierr/ierr
       common/istop/istop
@@ -2596,6 +2612,7 @@ c
       istpl=1
       ierr=0
       istep=0
+      s11=0
       istop=0
       ipolc=0
 c
@@ -5749,7 +5766,9 @@ c
       common/derip11/dpsi11dr,dpsi11dz,psinv11,bvx11,bvy11,bvz11
       common/istep/istep
       common/rwmax/rwmax
+      common/dsds/dst
       common/ss/st
+      common/s11/s11
       common/dnpar/dnpara
       common/atjki/tauv,alphav
       common/waist/w0csi,w0eta
@@ -5923,6 +5942,10 @@ c
           dx=ywrk(1,j,k)-ywrk(1,1,1)
           dy=ywrk(2,j,k)-ywrk(2,1,1)
           dz=ywrk(3,j,k)-ywrk(3,1,1)
+c
+c          dirx=dery0(1)
+c          diry=dery0(2)
+c          dirz=dery0(3)
 c
           dirx=ywrk(4,j,k)
           diry=ywrk(5,j,k)
@@ -6193,15 +6216,11 @@ c     call fast Fourier transform 2D
         do i=1,nxgrid
           ij=nxgrid*(i-1)+j
           aecompl(i,j)=exp(ui*(ak0*srint(ij)+ui*zwint(ij)))
-c     .                 *0.25*(1-cos((2*pi*(i-1))/(nxgrid-1)))*
-c     .                       (1-cos((2*pi*(j-1))/(nxgrid-1)))
           aemodel(i,j)=exp(ui*(aac*xgridv(i)**2+bbc*ygridv(j)**2+
      .                 ccc*xgridv(i)*ygridv(j)))    
           argcos=pi*(max(rho0,min(sqrt(xgridv(i)**2+ygridv(j)**2)
      .           ,rhomax))-rho0)/(rhomax-rho0)
           adiff(i,j)=(aecompl(i,j)-aemodel(i,j))*0.5*(1+cos(argcos))
-c     .              *0.25*(1-cos((2*pi*(i-1))/(nxgrid-1)))
-c     .                   *(1-cos((2*pi*(j-1))/(nxgrid-1)))
         end do
       end do
       do j=1,2*nxgrid-1
@@ -6256,8 +6275,9 @@ c controllare offset griglia
       write(82,*) ''
       adev2=sqrt(adev2/nindex)
       adev=adev/nindex
-      write(83,*) istep,adev2,adev,aar,bbr,ccr,aaw,bbw,ccw,akw,bkw,ckw,
-     .            dk1,dk2,wcsi,weta,rcicsi,rcieta
+      write(83,*) istep,s11/100.,istep*dst/100.,adev2,adev,
+     .            aar,bbr,ccr,aaw,bbw,ccw,akw,bkw,ckw,dk1,dk2,
+     .            wcsi,weta,rcicsi,rcieta
       return
  99   format(i5,22(1x,e16.8e3))
 111   format(3i5,30(1x,e16.8e3))