corrected psi grid leading dimension in JINTRAC interface. Committed new branch for updated wall reflection algorithm.

src/gray.f

@@ -2,7 +2,6 @@
       common/istop/istop
       common/ierr/ierr
       common/igrad/igrad
-      common/iovmin/iopmin,iowmin
       common/mode/sox
       common/p0/p0mw
       common/powrfl/powrfl
@@ -35,7 +34,7 @@ c     postprocessing
       currtott=currtot
       powtr=p0mw-pabstot
 
-      if (iowmin.eq.2.and.ipass.gt.1) then
+      if (ipass.gt.1) then
 c     second pass into plasma
         p0mw1=p0mw
         igrad=0
@@ -169,14 +168,18 @@ c
 c
       subroutine after_onestep(i,istop)
       implicit real*8 (a-h,o-z)
+      complex*16 ext,eyt,extr,eytr,exin2,eyin2
       parameter(jmx=31,kmx=36,nmx=8000)
       parameter(taucr=12.0d0,pi=3.14159265358979d0,cvdr=pi/180.0d0)
       dimension psjki(jmx,kmx,nmx),ppabs(jmx,kmx,nmx),ccci(jmx,kmx,nmx)
       dimension iiv(jmx,kmx),tauv(jmx,kmx,nmx),alphav(jmx,kmx,nmx)
-      dimension iop(jmx,kmx),iow(jmx,kmx)
+      dimension iop(jmx,kmx),iow(jmx,kmx),ihcd(jmx,kmx),istore(jmx,kmx)
-      dimension tau1v(jmx,kmx),yyrfl(jmx,kmx,6)
+      dimension tau1v(jmx,kmx),yyrfl(jmx,kmx,6),y(6),dery(6)
-      dimension xv(3),anv(3),xvrfl(3),anvrfl(3)
+      dimension ext(jmx,kmx,0:3),eyt(jmx,kmx,0:3)
-      
+      dimension xv(3),anv(3),xvrfl(3),anvrfl(3),xvvac(3),anw(3),anwcl(3)
+      dimension xvjk(jmx,kmx),anvjk(jmx,kmx)
+      dimension ywrk(6,jmx,kmx),ypwrk(6,jmx,kmx)
+c
       common/pcjki/ppabs,ccci
       common/atjki/tauv,alphav
       common/tau1v/tau1v
@@ -187,7 +190,8 @@ c
       common/istgr/istpr,istpl
 c
       common/iiv/iiv
-      common/iov/iop,iow
+      common/iov/iop,iow,ihcd,istore
+      common/refln/anwcl,jclosest
       common/psjki/psjki
       common/psival/psinv
       common/psinv11/psinv11
@@ -199,8 +203,8 @@ c
 c
       common/p0/p0mw
       common/pol0/psipol0,chipol0
-      common/ipol/ipolc
+      common/polcof/psipol,chipol
-      common/iovmin/iopmin,iowmin
+      common/evt/ext,eyt
       common/densbnd/psdbnd
       common/yyrfl/yyrfl
       common/powrfl/powrfl
@@ -209,9 +213,12 @@ c
       common/dsds/dst
       common/index_rt/index_rt
       common/ipass/ipass
-      common/rwallm/rwallm
       common/bound/zbmin,zbmax
-
+      common/limiter/rlim,zlim,nlim
+c
+      common/igrad/igrad
+      common/wrk/ywrk,ypwrk
+c
       pabstot=0.0d0
       currtot=0.0d0
       taumn=1d+30
@@ -219,6 +226,7 @@ c
       psinv11=1.0d0
       iopmin=100
       iowmin=100
+      iowmax=0
 c
       do j=1,nrayr
         kkk=nrayth
@@ -242,7 +250,7 @@ c              exit
           zzm=xv(3)/100.d0
           if(j.eq.1) rrm11=rrm
 
-          if (iwarm.gt.0.and.i.gt.1) then
+          if (iwarm.gt.0.and.i.gt.1.and.ihcd(j,k).gt.0) then
             if(psinv.ge.0.and.psinv.le.1.0d0.and.
      .         zzm.ge.zbmin.and.zzm.le.zbmax) then
               call pabs_curr(i,j,k)
@@ -257,87 +265,90 @@ c              exit
           end if 
           call print_output(i,j,k)
 
-          if(i.gt.1.and.psinv.ge.0.and.psinv.lt.psdbnd) 
+cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
-     .     iop(j,k)=1
+
-          if(iop(j,k).eq.1.and.
+          if (mod(iop(j,k),2).eq.0 .and. inside_plasma(rrm,zzm)) then
-     .       (psinv.ge.psdbnd.or.
+            iop(j,k)=iop(j,k)+1
-     .        (psinv.lt.1.0d0.and.(zzm.lt.zbmin.or.zzm.gt.zbmax))))
+            call pol_limit(ext(j,k,iop(j,k)),eyt(j,k,iop(j,k)))
-     .     iop(j,k)=2
+
-c         iov=0 initially, iov=1 first entrance in plasma, 
+            if (ipass.gt.1 .and. index_rt.eq.1 .and. 
-c         iov=2 first exit from plasma, iov=3 after 2nd entrance into plasma
+     .          iowmax.gt.1 .and. istore(j,k).eq.0) then
+              istore(j,k)=istore(j,k)+1
+              yyrfl(j,k,1:3)=xv
+              yyrfl(j,k,4:6)=anv
+              ihcd(j,k)=0
+            end if
+          else if (mod(iop(j,k),2).eq.1.and..not.inside_plasma(rrm,zzm))
+            iop(j,k)=iop(j,k)+1
+            call pol_limit(ext(j,k,iop(j,k)),eyt(j,k,iop(j,k)))
+          end if
           
-          if(index_rt.eq.1) then
+          if (ipass.gt.1) then
-          if(j.eq.1) then
+            if (iow(j,k).eq.0 .and. inside(rlim,zlim,nlim,rrm,zzm)) then
-            psinv11=psinv
+              iow(j,k)=1
-            if(ipolc.eq.0.and.iop(j,k).eq.1) then
+            else if (iow(j,k).eq.1 .and. 
-              call pol_limit(qqin,uuin,vvin)
+     .               .not.inside(rlim,zlim,nlim,rrm,zzm))
-              ipolc=1
+              iow(j,k)=2
-               qqa=cos(2.0d0*psipol0*cvdr)*cos(2.0d0*chipol0*cvdr)
+              if (inside_plasma(rrm,zzm)) then
-               uua=sin(2.0d0*psipol0*cvdr)*cos(2.0d0*chipol0*cvdr)
+                iop(j,k)=iop(j,k)+1
-               vva=sin(2.0d0*chipol0*cvdr)
+                call pol_limit(ext(j,k,iop(j,k)),eyt(j,k,iop(j,k)))
-               powa=0.5d0*(1.0d0+vva*vvin+uua*uuin+qqa*qqin)
-c              p0mw=p0mw*powa
-c              print*,' '
-c              print*,'Coupled power fraction =',powa
-c              print*,' '
-c              print*,'Input coupled power (MW) =',p0mw
-c              print*,' '
-            end if
-            if (iop(j,k).eq.2.and.rrm.le.rcen.and.ipass.gt.1
-     .          .and.ipolc.eq.1) then
-              call pol_limit(qqout,uuout,vvout)
-              ipolc=2
-              call wall_refl(xvrfl,anvrfl,qqrfl,uurfl,vvrfl,irfl)
-              strfl11=i*dst+dstvac
-              call pol_limit(qqin2,uuin2,vvin2)
-              if(irfl.gt.0) then
-                powrfl=0.5d0*(1.0d0+vvrfl*vvin2+uurfl*uuin2+qqrfl*qqin2)
-              else
-                powrfl=0.5d0*(1.0d0+vvout*vvin2+uuout*uuin2+qqout*qqin2)
               end if
-              write(6,*) 'Reflected power fraction =',powrfl
+              call wall_refl(xv-dst*anv,anv,ext(j,k,iop(j,k)),
-              iop(j,k)=3
+     .                eyt(j,k,iop(j,k)),xvrfl,anvrfl,extr,eytr,anw,irfl)
-              yyrfl(j,k,1)=xvrfl(1)
+              istore(j,k)=istore(j,k)+1
-              yyrfl(j,k,2)=xvrfl(2)
+              yyrfl(j,k,1:3)=xvrfl
-              yyrfl(j,k,3)=xvrfl(3)
+              yyrfl(j,k,4:6)=anvrfl
-              yyrfl(j,k,4)=anvrfl(1)
-              yyrfl(j,k,5)=anvrfl(2)
-              yyrfl(j,k,6)=anvrfl(3)
               tau1v(j,k)=tauv(j,k,iiv(j,k))
+              ext(j,k,iop(j,k))=extr
+              eyt(j,k,iop(j,k))=eytr
+              if (j.lt.jclosest) then
+                jclosest=j
+                anwcl=anw
+              end if
+              xv=xvrfl
+              anv=anvrfl
+              rrm=1.d-2*sqrt(xv(1)**2+xv(2)**2)
+              zzm=1.d-2*xv(3)
+              ywrk(1:3,j,k)=xv
+              ywrk(4:6,j,k)=anv
+              igrad=0
+              call gwork(j,k)
+              if (inside_plasma(rrm,zzm)) then
+                iop(j,k)=iop(j,k)+1
+                call pol_limit(ext(j,k,iop(j,k)),eyt(j,k,iop(j,k)))
+                if (index_rt.eq.1) ihcd(j,k)=0
+                end if
+              end if
             end if
-          else
-            if(iop(j,k).eq.2.and.rrm.le.rcen.and.ipass.gt.1)  then
-              call wall_refl(xvrfl,anvrfl,qqrfl,uurfl,vvrfl,irfl)
-              iop(j,k)=3
-              yyrfl(j,k,1)=xvrfl(1)
-              yyrfl(j,k,2)=xvrfl(2)
-              yyrfl(j,k,3)=xvrfl(3)
-              yyrfl(j,k,4)=anvrfl(1)
-              yyrfl(j,k,5)=anvrfl(2)
-              yyrfl(j,k,6)=anvrfl(3)
-              tau1v(j,k)=tauv(j,k,iiv(j,k))
-            end if  
-          end if   
           end if
 
+          if(index_rt.eq.1 .and. j.eq.1) psinv11=psinv
           if(iop(j,k).lt.iopmin) iopmin=iop(j,k)
+          if(iow(j,k).lt.iowmin) iowmin=iow(j,k)
+          if(iow(j,k).gt.iowmax) iowmax=iow(j,k)
+
+          xvjk(j,k)=xv
+          anvjk(j,k)=anv
 
         end do
       end do
 
-      if(iwarm.gt.0.and.taumn.lt.1d+30.and.taumn.gt.taucr) istop=1
+      if(jclosest.le.nrayr) then
+        aknmin=1.0d0
+        do j=1,nrayr
+          kkk=nrayth
+          if(j.eq.1) kkk=1
+          do k=1,kkk
+            akdotn=dot_product(anvjk(j,k),anwcl)
+            if(akdotn.lt.aknmin) aknmin=akdotn
+          end do
+        end do
+      else
+        aknmin=-1.0d0
+      end if
 
-      psimin=psjki(1,1,i)
+cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
-      if(nrayr.gt.1)
-     .    psimin=min(psimin,minval(psjki(2:nrayr,1:nrayth,i)))
-      if(psimin.gt.1.0d0.and.rrm11.gt.rcen.and.index_rt.gt.1)
-     .    istop=1
 
-      if(rrm11.lt.rwallm.and.ipass.eq.1) istop=1
-      if(iopmin.eq.2.and.ipass.eq.1) istop=1
-      
-      if(iopmin.eq.3) istop=1
-     
 c     print ray positions for j=nrayr in local reference system
 
       istpr=istpr+1
@@ -353,7 +364,82 @@ c        print*,'istep = ',i
 c
       if (istpl.eq.istpl0) istpl=0
       istpl=istpl+1
-      
+
+
+cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
+
+c     single pass is stopped when all the rays have crossed the plasma
+c         or complete absorption has occurred
+c     same for successive passes of multi-pass simulations (here exit
+c         from vessel is detected too
+c     first pass in multi-pass simulation is stopped when at least one
+c         ray has reflected and all rays are directed away from
+c         reflection point, or when no reflection has occurred and
+c         central ray re-enters the plasma
+      if((ipass.eq.1 .and. ((iopmin.gt.1) .or.
+     .                      (taumn.lt.1d+30.and.taumn.gt.taucr)))
+     .   .or.(index_rt.gt.1 .and. (iopmin.gt.1 .or. iowmin.gt.1 .or.
+     .                     (taumn.lt.1d+30.and.taumn.gt.taucr)))) then
+        istop=1
+      else if(ipass.gt.1 .and. index_rt.eq.1 .and.
+     .   ((iowmin.gt.1 .and. aknmin.gt.0) .or.
+     .    (iowmax.le.1 .and. iop(1,1).gt.2)) then
+c     flag second pass mode coupling as unset
+        powrfl=-1.0d0
+        qqout=0.0d0
+        uuout=0.0d0
+        vvout=0.0d0
+        do j=1,nrayr
+          kkk=nrayth
+          if(j.eq.1) kkk=1
+          do k=1,kkk
+c     store missing initial conditions for the second pass
+            if (istore(j,k).eq.0) then
+              istore(j,k)=istore(j,k)+1
+              yyrfl(j,k,1:3)=xvjk(j,k)
+              yyrfl(j,k,4:6)=anvjk(j,k)
+              tau1v(j,k)=tauv(j,k,iiv(j,k))
+            end if
+c     determine mode coupling at the plasma boundary
+            if (powrfl.lt.0.0d0) then
+              call vacuum_rt(xvjk(j,k),anvjk(j,k),xvvac,ivac)
+c     look for first ray hitting the plasma, starting from the central
+c     and evaluate polarization
+              if (ivac.eq.1) then
+                y(1:3)=xvjk(j,k)
+                y(4:6)=anvjk(j,k)
+                call fwork(y,dery)
+                call pol_limit(exin2,eyin2)
+                call stokes(exin2,eyin2,qqin2,uuin2,vvin2)
+       powloop: do j1=1,nrayr
+                  kkkk=nrayth
+                  if(j1.eq.1) kkkk=1
+                  do k1=1,kkkk
+c     look for first ray which completed the first pass in the plasma
+                    if (iop(j1,k1).gt.1) then
+c     if found, use its polarization state to compute mode coupling
+                      call stokes(ext(j1,k1,2),eyt(j1,k1,2),
+     .                            qqout,uuout,vvout)
+                      exit powloop
+                    end if
+                  end do
+                end do
+c     if no ray completed a first pass in the plasma, use central ray
+c     initial polarization (possibly reflected)
+                if (qqout.le.0.0d0) then
+                  call stokes(ext(1,1,0),eyt(1,1,0),qqout,uuout,vvout)
+                end if
+                powrfl=0.5d0*(1.0d0+vvout*vvin2+
+     .                              uuout*uuin2+qqout*qqin2)
+              end if
+            end if
+          end do
+        end do
+        strfl11=i*dst
+        write(6,*) 'Reflected power fraction =',powrfl
+        istop=1
+      end if
+
       return
       end
 c
@@ -2605,16 +2691,18 @@ c   spline interpolation of H(lambda,rhop) and dH/dlambda
       parameter(jmx=31,kmx=36,nmx=8000)
       dimension psjki(jmx,kmx,nmx)
       dimension tauv(jmx,kmx,nmx),alphav(jmx,kmx,nmx)
-      dimension pdjki(jmx,kmx,nmx),ppabs(jmx,kmx,nmx)
+      dimension pdjki(jmx,kmx,nmx),ppabs(jmx,kmx,nmx),tau1v(jmx,kmx)
       dimension currj(jmx,kmx,nmx),didst(jmx,kmx,nmx),ccci(jmx,kmx,nmx)
-      dimension iiv(jmx,kmx),iop(jmx,kmx),iow(jmx,kmx),tau1v(jmx,kmx)
+      dimension iiv(jmx,kmx),iop(jmx,kmx),iow(jmx,kmx),ihcd(jmx,kmx)
+      dimension istore(jmx,kmx),anwcl(3)
       parameter(tmax=5,npts=500)
       real*8 ttv(npts+1),extv(npts+1)
+
       common/ttex/ttv,extv
 c
       common/warm/iwarm,ilarm
       common/iiv/iiv
-      common/iov/iop,iow
+      common/iov/iop,iow,ihcd,istore
       common/psjki/psjki
       common/atjki/tauv,alphav
       common/dpjjki/pdjki,currj
@@ -2623,11 +2711,13 @@ c
       common/nray/nrayr,nrayth
       common/nstep/nstep
       common/tau1v/tau1v
+      common/refln/anwcl,jclosest
 c
       if(nstep.gt.nmx) nstep=nmx
       if(nrayr.gt.jmx) nrayr=jmx
       if(nrayth.gt.kmx) nrayth=kmx
-      
+      jclosest=nrayr+1
+      anwcl(1:3)=0.0d0
 c
       do i=1,nstep
         do k=1,nrayth
@@ -2643,6 +2733,8 @@ c
             iiv(j,k)=1
             iop(j,k)=0
             iow(j,k)=0
+            ihcd(j,k)=1
+            istore(j,k)=0
             tau1v(j,k)=0.0d0
           end do
         end do
@@ -2668,10 +2760,11 @@ c
       dimension tauv(jmx,kmx,nmx),alphav(jmx,kmx,nmx)
       dimension pdjki(jmx,kmx,nmx),ppabs(jmx,kmx,nmx)
       dimension currj(jmx,kmx,nmx),didst(jmx,kmx,nmx),ccci(jmx,kmx,nmx)
-      dimension iiv(jmx,kmx),iop(jmx,kmx),iow(jmx,kmx)
+      dimension iiv(jmx,kmx),iop(jmx,kmx),iow(jmx,kmx),ihcd(jmx,kmx)
+      dimension istore(jmx,kmx)
 
       common/iiv/iiv
-      common/iov/iop,iow
+      common/iov/iop,iow,ihcd,istore
       common/psjki/psjki
       common/atjki/tauv,alphav
       common/dpjjki/pdjki,currj
@@ -2694,6 +2787,7 @@ c
             iiv(j,k)=1
             iop(j,k)=0
             iow(j,k)=0
+            ihcd(j,k)=1
           end do
         end do
       end do
@@ -2710,14 +2804,12 @@ c
       common/istgr/istpr,istpl
       common/ierr/ierr
       common/istop/istop
-      common/ipol/ipolc
+
-c
       istpr=0
       istpl=1
       ierr=0
       istep=0
       istop=0
-      ipolc=0
 c
       return
       end
@@ -3841,6 +3933,7 @@ c
       dimension dffiu(jmx),ddffiu(jmx)
       dimension grad2(jmx,kmx),dgrad2v(ndimm,jmx,kmx)
       dimension gri(3,jmx,kmx),ggri(3,3,jmx,kmx)
+      complex*16 ext(jmx,kmx,0:3),eyt(jmx,kmx,0:3)
       complex*16 ui,sss,ddd,phic,qi1,qi2,tc,ts,qqxx,qqxy,qqyy
       complex*16 dqi1,dqi2,dqqxx,dqqyy,dqqxy
       complex*16 d2qi1,d2qi2,d2qqxx,d2qqyy,d2qqxy
@@ -3862,6 +3955,7 @@ c
       common/gradjk/gri
       common/ggradjk/ggri
       common/ddd/dd,an2s,an2,fdia,bdotgr,ddi,ddr11
+      common/evt/ext,eyt
 c
       ui=(0.0d0,1.0d0)
       csth=anz0c
@@ -4060,6 +4154,8 @@ c
           ypwrk0(4,j,k) = dgr2x/an0/2.0d0
           ypwrk0(5,j,k) = dgr2y/an0/2.0d0
           ypwrk0(6,j,k) = dgr2z/an0/2.0d0
+c
+          call pol_limit(ext(j,k,0),eyt(j,k,0))
 c
           grad2(j,k)=gr2
           dgrad2v(1,j,k)=dgr2x
@@ -4121,6 +4217,7 @@ c
       dimension dffiu(jmx),ddffiu(jmx)
       dimension grad2(jmx,kmx),dgrad2v(ndimm,jmx,kmx)
       dimension gri(3,jmx,kmx),ggri(3,3,jmx,kmx)
+      complex*16 ext(jmx,kmx,0:3),eyt(jmx,kmx,0:3)
 c
       common/nray/nrayr,nrayth
       common/rwmax/rwmax
@@ -4135,6 +4232,7 @@ c
       common/gradjk/gri
       common/ggradjk/ggri
       common/ddd/dd,an2s,an2,fdia,bdotgr,ddi,ddr11
+      common/evt/ext,eyt
 c
       csth=anz0c
       snth=sqrt(1.0d0-csth**2)
@@ -4217,6 +4315,8 @@ c
           ypwrk0(4,j,k) = 0.0d0
           ypwrk0(5,j,k) = 0.0d0
           ypwrk0(6,j,k) = 0.0d0
+c
+          call pol_limit(ext(j,k,0),eyt(j,k,0))
 c
           do iv=1,3
             gri(iv,j,k)=0.0d0
@@ -4277,6 +4377,7 @@ c
       dimension dffiu(jmx),ddffiu(jmx)
       dimension grad2(jmx,kmx),dgrad2v(ndimm,jmx,kmx)
       dimension gri(3,jmx,kmx),ggri(3,3,jmx,kmx)
+      complex*16 ext(jmx,kmx,0:3),eyt(jmx,kmx,0:3)
 c
       common/nray/nrayr,nrayth
       common/wrk/ywrk0,ypwrk0
@@ -4287,6 +4388,7 @@ c
       common/ggradjk/ggri
       common/ddd/dd,an2s,an2,fdia,bdotgr,ddi,ddr11
       common/yyrfl/yyrfl
+      common/evt/ext,eyt
 
       do j=1,nrayr
         do k=1,nrayth
@@ -4316,6 +4418,8 @@ c
           ypwrk0(4,j,k) = 0.0d0
           ypwrk0(5,j,k) = 0.0d0
           ypwrk0(6,j,k) = 0.0d0
+c
+          call pol_limit(ext(j,k,0),eyt(j,k,0))
 c
           do iv=1,3
             gri(iv,j,k)=0.0d0
@@ -6564,17 +6668,14 @@ c
       return
       end
 
-      subroutine pol_limit(qq,uu,vv)
+      subroutine pol_limit(ext,eyt)
       implicit none
-      integer*4 ipolc
       real*8 bv(3),anv(3)
       real*8 anx,any,anz,anpl,anpr,yg,xe2om,ye2om,xe2xm,ye2xm
       real*8 an2,an,anxy,sngam,csgam,csg2,sng2,ffo,ffx,ffo2,ffx2
       real*8 deno,denx,anpl2,dnl,del0
-      real*8 uuom,vvom,qqom,uuxm,vvxm,qqxm,ellom,ellxm,qq,uu,vv
-      real*8 aaom,bbom,llmom,aaxm,bbxm,llmxm,psiom,psixm,chiom,chixm
       real*8 pi,beta0,alpha0,gam
-      real*8 sox,psipol,chipol
+      real*8 sox
       complex*16 ui,exom,eyom,exxm,eyxm,ext,eyt
       parameter(ui=(0.0d0,1.0d0),pi=3.14159265358979d0)
 c
@@ -6582,11 +6683,7 @@ c
       common/nplr/anpl,anpr
       common/ygyg/yg
       common/bb/bv
-      common/angles/alpha0,beta0
       common/mode/sox
-      common/polcof/psipol,chipol
-      common/ipol/ipolc
-      common/evt/ext,eyt
 c
       anx=anv(1)
       any=anv(2)
@@ -6614,103 +6711,139 @@ c
 c
       exom=(ffo*csgam-ui*anpl*sngam)/sqrt(deno)
       eyom=(-ffo*sngam-ui*anpl*csgam)/sqrt(deno)
-      qqom=abs(exom)**2-abs(eyom)**2
-      uuom=2.0d0*dble(exom*dconjg(eyom))
-      vvom=2.0d0*dimag(exom*dconjg(eyom))
-      llmom=sqrt(qqom**2+uuom**2)
-      aaom=sqrt((1+llmom)/2.0d0)
-      bbom=sqrt((1-llmom)/2.0d0)
-      ellom=bbom/aaom
-      psiom=0.5d0*atan2(uuom,qqom)*180.d0/pi
-      chiom=0.5d0*asin(vvom)*180.d0/pi
 c
       exxm=(ffx*csgam-ui*anpl*sngam)/sqrt(denx)
       eyxm=(-ffx*sngam-ui*anpl*csgam)/sqrt(denx)
-      qqxm=abs(exxm)**2-abs(eyxm)**2
-      uuxm=2.0d0*dble(exxm*dconjg(eyxm))
-      vvxm=2.0d0*dimag(exxm*dconjg(eyxm))
-      llmxm=sqrt(qqxm**2+uuxm**2)
-      aaxm=sqrt((1+llmxm)/2.0d0)
-      bbxm=sqrt((1-llmxm)/2.0d0)
-      ellxm=bbxm/aaxm
-      psixm=0.5d0*atan2(uuxm,qqxm)*180.d0/pi
-      chixm=0.5d0*asin(vvxm)*180.d0/pi
 c
       if (sox.lt.0.0d0) then
-        psipol=psiom
-        chipol=chiom
         ext=exom
         eyt=eyom
-        qq=qqom
-        vv=vvom
-        uu=uuom
       else
-        psipol=psixm
-        chipol=chixm
         ext=exxm
         eyt=eyxm
-        qq=qqxm
-        vv=vvxm
-        uu=uuxm
       endif
+
       gam=atan(sngam/csgam)*180.d0/pi
 
       return
  111  format(20(1x,e12.5))      
       end
 
-
+      subroutine stokes(ext,eyt,qq,uu,vv)
-
-      subroutine wall_refl(xvrfl,anvrfl,qqtr,uutr,vvtr,irfl)
       implicit none
-      integer*4 ivac,irfl
+      complex*16 ext,eyt
-      real*8 anv(3),xv(3),xvrfl(3)
+      real*8 qq,uu,vv
-      real*8 walln(3),anvrfl(3),vv1(3),vv2(3),vv3(3),xvout(3)
+      qq=abs(ext)**2-abs(eyt)**2
-      real*8 uutr,vvtr,qqtr,qq,uu,vv
+      uu=2.0d0*dble(ext*dconjg(eyt))
+      vv=2.0d0*dimag(ext*dconjg(eyt))
+      end subroutine stokes
+      
+      subroutine polellipse(qq,uu,vv,psipol,chipol)
+      implicit none
+      real*8 qq,uu,vv,psipol,chipol
+c      real*8 llm,aa,bb,ell
       real*8 pi
-      real*8 psipol,chipol,psitr,chitr
+      parameter(pi=3.14159265358979d0)
+c      llm=sqrt(qq**2+uu**2)
+c      aa=sqrt((1+llm)/2.0d0)
+c      bb=sqrt((1-llm)/2.0d0)
+c      ell=bb/aa
+      psipol=0.5d0*atan2(uu,qq)*180.d0/pi
+      chipol=0.5d0*asin(vv)*180.d0/pi
+      end subroutine polellipse
+
+      
+      logical function inside_plasma(rrm,zzm)
+      implicit none
+      real*8 rrm,zzm,psdbnd,psinv,zbmin,zbmax
+      integer iequil
+      
+      common/densbnd/psdbnd
+      common/bound/zbmin,zbmax
+      common/psival/psinv
+      common/iieq/iequil
+
+      if(iequil.eq.1) then
+        call equian(rrm,zzm)
+      else
+        call equinum(rrm,zzm)
+      end if
+
+      if (psinv.ge.0.0d0.and.psinv.lt.psdbnd) then
+        if (psinv.lt.1.0d0.and.zzm.lt.zbmin.or.zzm.gt.zbmax) then
+          inside_plasma=.false.
+        else
+          inside_plasma=.true.
+        end if
+      else
+        inside_plasma=.false.
+      end if
+      
+      end function inside_plasma
+
+
+      subroutine wall_refl(xv,anv,ext,eyt,xvrfl,anvrfl,extr,eytr,walln,
+     .                     irfl)
+      implicit none
+      integer*4 irfl
+      real*8 anv(3),anv0(3),xv(3),xvrfl(3)
+      real*8 walln(3),anvrfl(3),vv1(3),vv2(3),vv3(3)
+      real*8 pi,smax
       complex*16 ui,extr,eytr,eztr,ext,eyt
       complex*16 evin(3),evrfl(3) 
       parameter(ui=(0.0d0,1.0d0),pi=3.14159265358979d0)
+      integer nbb,nlim
+      parameter(nbb=5000)
+      real*8 rlim(nbb),zlim(nbb)
 
-      common/xv/xv
+      common/limiter/rlim,zlim,nlim
-      common/anv/anv
-      common/polcof/psipol,chipol
-      common/evt/ext,eyt
-      common/wrefl/walln
       
-      anv=anv/sqrt(anv(1)**2+anv(2)**2+anv(3)**2)
+      anv0=anv/sqrt(anv(1)**2+anv(2)**2+anv(3)**2)
+      rrm=1.d-2*sqrt(xv(1)**2+xv(2)**2)
+      zzm=1.d-2*xv(3)
       
 c     computation of reflection coordinates and normal to the wall
+      call inters_linewall(xv/1.d2,anv0,rlim(1:nlim),zlim(1:nlim),
+     .                     nlim,smax,walln)
+      smax=smax*1.d2
+      xvrfl=xv+smax*anv0
       irfl=1
-      ivac=1
+      if (.not.inside(rlim,zlim,nlim,rrm,zzm)) then
-      call vacuum_rt(xv,xvout,ivac)
+        ! first wall interface is outside-inside
-      
+        if (dot_product(walln,walln)<tiny(walln)) then
-      if(ivac.lt.0) then
+          ! wall never hit
-        irfl=0
+          xvrfl=xv
-        xvrfl=xvout
+          anvrfl=anv0
-        xv=xvout
+          extr=ext
-        anvrfl=anv
+          eytr=eyt
-        return
+          irfl=0
-      end if  
+          return
+        end if
+        ! search second wall interface (inside-outside)
+        call inters_linewall(xvrfl/1.d2,anv0,rlim(1:nlim),zlim(1:nlim),
+     .                       nlim,smax,walln)
+        smax=smax*1.d2
+        xvrfl=xvrfl+smax*anv0
+        irfl=2
+      end if
       
 c     rotation matrix from local to lab frame
-      vv1(1)=anv(2)
+      vv1(1)=anv0(2)
-      vv1(2)=-anv(1)
+      vv1(2)=-anv0(1)
       vv1(3)=0.0d0
-      vv2(1)=anv(1)*anv(3)
+      vv2(1)=anv0(1)*anv0(3)
-      vv2(2)=anv(2)*anv(3)
+      vv2(2)=anv0(2)*anv0(3)
-      vv2(3)=-anv(1)*anv(1)-anv(2)*anv(2)
+      vv2(3)=-anv0(1)*anv0(1)-anv0(2)*anv0(2)
       vv1=vv1/sqrt(vv1(1)**2+vv1(2)**2+vv1(3)**2)
       vv2=vv2/sqrt(vv2(1)**2+vv2(2)**2+vv2(3)**2)
-      vv3=anv
+      vv3=anv0
       
       evin=ext*vv1+eyt*vv2
 c     wave vector and electric field after reflection in lab frame
-      anvrfl=anv-2.0d0*
+      anvrfl=anv0-2.0d0*
-     .       (anv(1)*walln(1)+anv(2)*walln(2)+anv(3)*walln(3))*walln
+     .       (anv0(1)*walln(1)+anv0(2)*walln(2)+anv0(3)*walln(3))*walln
       evrfl=-evin+2.0d0*
-     .     (evin(1)*walln(1)+evin(2)*walln(2)+evin(3)*walln(3))*walln
+     .       (evin(1)*walln(1)+evin(2)*walln(2)+evin(3)*walln(3))*walln
 
       vv1(1)=anvrfl(2)
       vv1(2)=-anvrfl(1)
@@ -6726,25 +6859,11 @@ c     wave vector and electric field after reflection in lab frame
       eytr=dot_product(vv2,evrfl)
       eztr=dot_product(vv3,evrfl)
 
-      qqtr=abs(extr)**2-abs(eytr)**2
-      uutr=2.0d0*dble(extr*dconjg(eytr))
-      vvtr=2.0d0*dimag(extr*dconjg(eytr))
-      psitr=0.5d0*atan2(uutr,qqtr)*180.d0/pi
-      chitr=0.5d0*asin(vvtr)*180.d0/pi
-
-      ivac=2
-      anv=anvrfl
-      xvrfl=xvout
-      xv=xvout
-      
-      call vacuum_rt(xv,xvout,ivac)
-      xv=xvout
-      
       return
  111  format(20(1x,e12.5))      
       end
       
-      subroutine vacuum_rt(xvstart,xvend,ivac)
+      subroutine vacuum_rt(xvstart,anv,xvend,ivac)
       use reflections, only : inters_linewall,inside
       implicit none
       integer*4 ivac
@@ -6752,23 +6871,24 @@ c     wave vector and electric field after reflection in lab frame
       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)
+      real*8 xv0(3),anv0(3)
       real*8 rlim(nbb),zlim(nbb)
+      logical plfound
       
-      common/wrefl/walln
       common/limiter/rlim,zlim,nlim
-      common/anv/anv
       common/dsds/dst      
       common/psival/psinv
       common/densbnd/psdbnd
       common/dstvac/dstvac
-c     ivac=1 first interface plasma-vacuum
+
-c     ivac=2 second interface vacuum-plasma after wall reflection
+c     ivac=1 plasma hit before wall reflection
-c     ivac=-1 second interface vacuum-plasma WITHOUT wall reflection
+c     ivac=2 wall hit before plasma
+c     ivac=-1 vessel (and thus plasma) never crossed
       
 !     the real argument associated to xvstart is in a common block
 !     used by fwork!!!
       xv0=xvstart
+      anv0=anv
       call inters_linewall(xv0/1.d2,anv,rlim(1:nlim),zlim(1:nlim),
      .                     nlim,smax,walln)
       smax=smax*1.d2
@@ -6793,29 +6913,27 @@ c     ivac=-1 second interface vacuum-plasma WITHOUT wall reflection
       i=0
       do
         st=i*dst
-        xvend=xv0+st*anv
+        xvend=xv0+st*anv0
-        y(1)=xvend(1)
+        rrm=1.d-2*sqrt(xvend(1)**2+xvend(2)**2)
-        y(2)=xvend(2)
+        zzm=1.d-2*xvend(3)
-        y(3)=xvend(3)
+        plfound=inside_plasma(rrm,zzm)
-        y(4)=anv(1)
+        if (st.ge.smax.or.plfound) exit
-        y(5)=anv(2)
-        y(6)=anv(3)
-        call fwork(y,dery)
-        if (st.ge.smax.or.(psinv.gt.0.0d0.and.psinv.lt.psdbnd)) exit
         i=i+1
       end do
       
-      if (st.lt.smax) then
+      if (plfound) then
-        ivac=-1
+        ivac=1
         dstvac=st
       else
+        ivac=2
         dstvac=smax
-        xvend=xv0+smax*anv
+        xvend=xv0+smax*anv0
       end if
 
 !     the real argument associated to xvstart is in a common block
 !     used by fwork!!!
       xvstart=xv0
+      anv=anv0
       
       return
  111  format(20(1x,e12.5))