diff --git a/src/gray.f b/src/gray.f
index 46505e2..c8ab213 100644
--- a/src/gray.f
+++ b/src/gray.f
@@ -153,6 +153,7 @@ c
       if(iprof.eq.0) write(6,*)  'ANALTYCAL PROFILES'
       if(ibeam.ge.1) write(6,*)  'LAUNCHER CASE : ',filenmbm
 
+      call print_tox
       end if
 c
 c       compute power and current density profiles for all rays
@@ -267,6 +268,7 @@ c         iov=2 first exit from plasma, iov=3 after 2nd entrance into plasma
           if(index_rt.eq.1) then
           if(j.eq.1) then
             psinv11=psinv
+            call calc_tox(i)
             if(ipolc.eq.0.and.iov(j,k).eq.1) then
               call pol_limit(qqin,uuin,vvin)
               ipolc=1
@@ -355,6 +357,77 @@ c
       end
 c
 c
+c
+      subroutine calc_tox(i)
+      implicit real*8 (a-h,o-z)
+c computation of O-X transmission function
+c calculated only for xg=(wpl/w)**2>0.5
+c values stored at maximum xg and at maximum transmission
+      parameter(pi=3.14159265358979d0)
+      dimension derxg(3),deryg(3),xv(3),anv(3)
+c
+      common/xgxg/xg
+      common/ygyg/yg
+      common/xv/xv
+      common/anv/anv
+      common/nplr/anpl,anpr
+      common/ddd/ddr,an2s,an2,fdia,bdotgr,ddi,ddr11
+      common/nprw/anprr,anpri
+      common/parwv/ak0,akinv,fhz
+      common/dxgyg/derxg,deryg
+      common/toxmxx/toxxmx,xgmax,ygxmx,anplxmx,anprxmx,thkbxmx,ixmx
+      common/toxmxt/toxmax,xgtmx,ygtmx,anpltmx,anprtmx,thkbtmx,itmx
+
+      anplopt=sqrt(yg/(yg+1.0d0))
+      thkb=1.8d2/pi*acos(anpl/sqrt(an2))
+      dlen=sqrt(derxg(1)**2+derxg(2)**2+derxg(3)**2)
+      if (dlen.gt.0.0d0) then
+        dlen=xg/dlen
+        tox=exp(-pi*ak0*dlen*sqrt(0.5d0*yg)*
+     .      (2*(1.0d0+yg)*(anplopt-abs(anpl))**2+anpr**2))
+      else
+        tox=0.0d0
+      end if
+      if (xg.gt.0.5d0.and.xg.le.1.0d0) then
+        if (xg.gt.xgmax) then
+          ixmx=i
+          toxxmx=tox
+          xgmax=xg
+          ygxmx=yg
+          anplxmx=anpl
+          anprxmx=anpr
+          thkbxmx=thkb
+        end if
+        if (tox.gt.toxmax) then
+          itmx=i
+          toxmax=tox
+          xgtmx=xg
+          ygtmx=yg
+          anpltmx=anpl
+          anprtmx=anpr
+          thkbtmx=thkb
+        end if
+        write(46,98) i,tox,xg,yg,anpl,anpr,thkb,xv/1.d2,anv,derxg,deryg
+      end if
+c
+      return
+ 98   format(1x,i8,30(1x,e16.8e3))
+      end
+c
+c
+c
+      subroutine print_tox
+      implicit real*8 (a-h,o-z)
+      common/toxmxx/toxxmx,xgmax,ygxmx,anplxmx,anprxmx,thkbxmx,ixmx
+      common/toxmxt/toxmax,xgtmx,ygtmx,anpltmx,anprtmx,thkbtmx,itmx
+      write(44,98) ixmx,toxxmx,xgmax,ygxmx,anplxmx,anprxmx,thkbxmx
+      write(45,98) itmx,toxmax,xgtmx,ygtmx,anpltmx,anprtmx,thkbtmx
+c
+      return
+ 98   format(1x,i8,30(1x,e16.8e3))
+      end
+c
+c
 c
       subroutine print_output(i,j,k)
       implicit real*8 (a-h,o-z)
@@ -2585,6 +2658,9 @@ c
       common/ierr/ierr
       common/istop/istop
       common/ipol/ipolc
+c used only for computation of O-X transmission function
+      common/toxmxx/toxxmx,xgmax,ygxmx,anplxmx,anprxmx,thkbxmx,ixmx
+      common/toxmxt/toxmax,xgtmx,ygtmx,anpltmx,anprtmx,thkbtmx,itmx
 c
       istpr=0
       istpl=1
@@ -2592,6 +2668,11 @@ c
       istep=0
       istop=0
       ipolc=0
+c
+      ixmx=0
+      itmx=0
+      xgmax=0.0d0
+      toxmax=0.0d0
 c
       return
       end