diff --git a/src/gray.f b/src/gray.f
index 55b8405..3e3296b 100644
--- a/src/gray.f
+++ b/src/gray.f
@@ -513,7 +513,7 @@ c
       write(33,*) ' #i j k sst x y R z psi tauv Npl alpha index_rt'
       write(12,*) ' #i sst psi w1 w2'
       write(7,*)'#Icd Pa Jphip dPdVp '//
-     .'rhotj rhotjava rhotp rhotpav drhotjava drhotpav '//
+     .'rhotj rhotjava rhotp rhotpav drhotjava drhotpav ratjamx '//
      .'ratjbmx stmx psipol chipol index_rt Jphimx dPdVmx drhotj drhotp'
       write(48,*) '#psi rhot Jphi Jcdb dPdV Icdins Pins P% index_rt'
 
@@ -543,7 +543,7 @@ c
       character*255 filenmeqq,filenmprf,filenmbm
       parameter(qe=4.8032d-10,me=9.1095d-28,vc=2.9979d+10)
       parameter(pi=3.14159265358979d0,cvdr=pi/180.0d0)
-      parameter(nmx=8000,nbb=1000)
+      parameter(nmx=8000,nbb=5000)
       real*8 rlim(nbb),zlim(nbb)
 c
       common/xgcn/xgcn
@@ -1066,7 +1066,7 @@ c
       implicit real*8 (a-h,o-z)
       parameter(nnw=501,nnh=501)
       parameter(pi=3.14159265358979d0)
-      parameter(nbb=1000)
+      parameter(nbb=5000)
 c      parameter(np=100)
       character*48 stringa
       dimension fpol(nnw),pres(nnw),qpsi(nnw)
@@ -2501,7 +2501,7 @@ c  ffhlam = Bmn/Bmx/fc integral_lambda^1 dlam/<sqrt(1-lam*B(rhop)/Bmx)>
       end do
 
       write(56,*)' #psi rhot_eq rhot_av |<B>| |Bmx| |Bmn|'//
-     .'  Area Vol |I_pl| <J_phi> qq fc'
+     .'  Area Vol |I_pl| <J_phi> qq fc ratioJa ratioJb'
 
       qqv(1)=qqv(2)
       vajphiav(1)=vajphiav(2)
@@ -2516,6 +2516,7 @@ c  ffhlam = Bmn/Bmx/fc integral_lambda^1 dlam/<sqrt(1-lam*B(rhop)/Bmx)>
         write(56,99) pstab(jp),rhot_eq,rhotqv(jp),
      .    bav(jp),bmxpsi(jp),bmnpsi(jp),varea(jp),vvol(jp), 
      .    vcurrp(jp),vajphiav(jp),qqv(jp),ffc(jp)
+     .    ,vratja(jp),vratjb(jp)
       end do 
 
       rarea=sqrt(varea(nintp)/pi)
@@ -6296,6 +6297,7 @@ c     of gaussian profile
       rhpj=frhopol(rhotjava)
       dpdvp=pabs*2.0d0/(sqrt(pi)*drhotjava*fdvdrhot(rhpp))
       ajphip=currt*2.0d0/(sqrt(pi)*drhotjava*fdadrhot(rhpj))
+      call ratioj(rhpj,ratjamx,ratjbmx,ratjplmx)
       
       call profwidth(nd,rtab,dpdv,rhotp,rhopp,dpdvmx,
      .                 drhotp,drhopp)
@@ -6303,7 +6305,7 @@ c     of gaussian profile
         call profwidth(nd,rtab,ajphiv,rhotjfi,rhopfi,ajmxfi,
      .                 drhotjfi,drhopfi)
         xps=rhopfi
-        call ratioj(rhopfi,ratjamx,ratjbmx,ratjplmx)
+c        call ratioj(rhopfi,ratjamx,ratjbmx,ratjplmx)
       else
         rhotjfi=0.0d0
         rhopfi=0.0d0
@@ -6370,17 +6372,17 @@ c     dPdV [MW/m^3], Jcd [MA/m^2]
 
       write(6,*)' '
       write(6,*)'#beta0 alpha0 Icd Pa Jphip dPdVp '//
-     .'rhotj rhotjava rhotp rhotpav drhotjava drhotpav '//
+     .'rhotj rhotjava rhotp rhotpav drhotjava drhotpav ratjamx '//
      .'ratjbmx stmx psipol chipol index_rt Jphimx dPdVmx drhotj drhotp'
       write(6,99) beta0,alpha0,currtka,pabstot,ajphip,dpdvp,
      .  rhotjfi,rhotjava,rhotp,rhotpav,
-     .  drhotjava,drhotpav,stmx,psipol,chipol,real(index_rt),
-     .  ajmxfi,dpdvmx,drhotjfi,drhotp
+     .  drhotjava,drhotpav,ratjamx,ratjbmx,stmx,psipol,chipol,
+     .  real(index_rt),ajmxfi,dpdvmx,drhotjfi,drhotp
      
       write(7,99) currtka,pabstot,ajphip,dpdvp,
      .  rhotjfi,rhotjava,rhotp,rhotpav,
-     .  drhotjava,drhotpav,ratjbmx,stmx,psipol,chipol,real(index_rt),
-     .  ajmxfi,dpdvmx,drhotjfi,drhotp
+     .  drhotjava,drhotpav,ratjamx,ratjbmx,stmx,psipol,chipol,
+     .  real(index_rt),ajmxfi,dpdvmx,drhotjfi,drhotp
 
       do i=1,nd
         if (ipec.eq.0) then
@@ -6747,7 +6749,7 @@ c     wave vector and electric field after reflection in lab frame
       implicit none
       integer*4 ivac
       integer nbb,nlim,i,imax
-      parameter(nbb=1000)
+      parameter(nbb=5000)
       real*8 st,rrm,zzm,psinv,dst,psdbnd,dstvac,smax
       real*8 anv(3),xvstart(3),xvend(3),walln(3),y(6),dery(6)
       real*8 xv0(3)