!-----------------------------------------------------------------------
!  boxinc -- find next intersection segment in 3d box grid
!
subroutine boxinc(zx,zy,zl,ix,iy,il,iin, &
     zxf,zyf,zlf,ixf,iyf,ilf,iinf, &
     idone,ilinside_pl,n_j,n_iz,iadr,zdlseg,il_halo)
  !-----------------------------------------------------------------------
  !
  use nbi_com
  use logmod
  !
  !============
  implicit none
  integer :: idirx,ixbdy,idiry,iybdy,idirl,ilbdy
  logical :: il_halo !.true.if 3d halo model 'on'
  !============
  real(fp) :: zincx,zincy,zincl,zdx,zxbdy,zdy,zybdy,zdl,zlbdy,zinc
  real(fp) :: zxnew,zynew,zlnew
  !============
  real(fp) :: zx,zy,zl,zx_in,zy_in,zl_in
  integer :: ix,iy,il
  logical iin
  !
  real(fp) :: zxf,zyf,zlf
  integer :: ixf,iyf,ilf
  logical iinf
  !
  logical idone,ilinside_pl,ilinside_pl_prev
  integer :: iadr,ii,inum
  real(fp) :: zdlseg,zll
  real(fp) :: rb(3),rb1(3),r(3),r_u(3),zr,zphi
  integer :: n_j,n_iz
  logical :: l_inside
  !
  !  input/output -- (zx,zy,zl) (real space coords in box cartesian coord. system) and
  !    (ix,iy,il) (index coords) and iin (inside/outside flag)
  !    get advanced towards the target.
  !
  !  input -- (zxf,zyf,zlf) and (ixf,iyf,ilf) are the target point
  !    towards which we move (iinf flags if this is inside or outside
  !    the box grid).
  !
  !  output -- idone=.true. when the target point has been reached.
  !            iadr=0 means this segment is outside the grid
  !                =<value> means this segment falls in the indicated
  !                 grid address.
  !            zdlseg = length of segment (cm)
  !
  !-----------------------------------------------------------------------
  !
  !  1.  compute address of starting zone if it is inside
  !

  zx_in=zx
  zy_in=zy
  zl_in=zl
  ilinside_pl=.false.
  if(iin) then
     iadr=ix + nxbox*(iy-1) + nxbox*nybox*(il-1)
  else
     iadr=0
  end if
  l_inside=(levmod_halo.eq.2.and.iadr.gt.0.and.il_halo)
  !
  !  2.  there is nothing to do if both start and finish zones are
  !      outside the box on the same side
  !
  zdlseg=0.0_fp
  idone=.false.
  !
  if((ix.eq.0).and.(ixf.eq.0)) idone=.true.
  if((ix.eq.(nxbox+1)).and.(ixf.eq.(nxbox+1))) idone=.true.
  if((iy.eq.0).and.(iyf.eq.0)) idone=.true.
  if((iy.eq.(nybox+1)).and.(iyf.eq.(nybox+1))) idone=.true.
  if((il.eq.0).and.(ilf.eq.0)) idone=.true.
  if((il.eq.(nlbox+1)).and.(ilf.eq.(nlbox+1))) idone=.true.
  if(idone) return
  !
  !  3.  we are done with one more zone to compute if the source and
  !      stop zones are the same
  !
  if((ix.eq.ixf).and.(iy.eq.iyf).and.(il.eq.ilf)) then
     idone=.true.
     zdlseg=sqrt((zxf-zx)**2+(zyf-zy)**2+(zlf-zl)**2)
     rb(1)=zxf
     rb(2)=zyf
     rb(3)=zlf
     ilinside_pl=.false.
     if(l_inside) then
        ! rb [cm]
        call box2car(ldep,rb,r) ! map  x,y,z coordinates of point in box coordinates to tokamak, [m] ->[cm]
        r=r*0.01_fp          !cm -> m
        call eq_rcyl(nsnccw*nsjdotb,r(1),r(2),zr,zphi) ![m,m] -> [m,rad]
        ! find transp zone n_j, and beam zone n_iz
        call nbi_zone_ind(zr,zphi,r(3),n_j,n_iz)
        if(n_j.lt.lep1.and.n_j.ge.lcentr.and.n_iz.gt.0.and. &
             n_iz.le.mimxbzf)then
           ilinside_pl=.true. !end point inside plasma

           zx=rb(1)
           zy=rb(2)
           zl=rb(3)
        end if
        !mg           return
        if(ilinside_pl) return
        !else end point is outside plasma
        rb(1)=zx
        rb(2)=zy
        rb(3)=zl
        call find_point_inside
     end if

     return
  end if

  !
  !  4.  we are in different zones -- look for next boundary to cross
  !
  zincx=epsinv
  zincy=epsinv
  zincl=epsinv
  !
  !  how far to an x bdy?
  if(ix.ne.ixf) then
     if(ix.lt.ixf) then
        idirx=1
        ixbdy=ix
     else
        idirx=-1
        ixbdy=ix-1
     end if
     zdx=2.0_fp*xboxhw/nxbox
     zxbdy=-xboxhw+ixbdy*zdx
     zincx=(zxbdy-zx)/(zxf-zx)
  end if
  !  how far to an y bdy?
  if(iy.ne.iyf) then
     if(iy.lt.iyf) then
        idiry=1
        iybdy=iy
     else
        idiry=-1
        iybdy=iy-1
     end if
     zdy=2.0_fp*yboxhw/nybox
     zybdy=-yboxhw+iybdy*zdy
     zincy=(zybdy-zy)/(zyf-zy)
  end if
  !  how far to an l bdy?
  if(il.ne.ilf) then
     if(il.lt.ilf) then
        idirl=1
        ilbdy=il
     else
        idirl=-1
        ilbdy=il-1
     end if
     zdl=(xlbox2-xlbox1)/nlbox
     zlbdy=xlbox1+ilbdy*zdl
     zincl=(zlbdy-zl)/(zlf-zl)
  end if
  !
  !  find minimum increment
  zinc=min(zincx,zincy,zincl)
  !
  !  increment zone indices
  if(zinc.eq.zincx) then
     ix=ix+idirx
  end if
  !
  if(zinc.eq.zincy) then
     iy=iy+idiry
  end if
  !
  if(zinc.eq.zincl) then
     il=il+idirl
  end if
  !
  !  test if still inside after increment
  iin=.true.
  if((ix.eq.0).or.(ix.eq.(nxbox+1))) iin=.false.
  if((iy.eq.0).or.(iy.eq.(nybox+1))) iin=.false.
  if((il.eq.0).or.(il.eq.(nlbox+1))) iin=.false.
  !
  !  compute new coordinates
  !
  zxnew=zx+zinc*(zxf-zx)
  zynew=zy+zinc*(zyf-zy)
  zlnew=zl+zinc*(zlf-zl)
  !
  zdlseg=sqrt((zxnew-zx)**2+(zynew-zy)**2+(zlnew-zl)**2)
  !
  zx=zxnew
  zy=zynew
  zl=zlnew
  rb(1)=zx
  rb(2)=zy
  rb(3)=zl
  !
  if(l_inside) then
     ! rb [cm]
     call box2car(ldep,rb,r) ! map  x,y,z coordinates of point in box coordinates to tokamak, [m] ->[cm]
     r=r*0.01_fp            !cm -> m
     call eq_rcyl(nsnccw*nsjdotb,r(1),r(2),zr,zphi) ![m,m] -> [m,rad]
     ! find transp zone n_j, and beam zone n_iz
     call nbi_zone_ind(zr,zphi,r(3),n_j,n_iz)
     if(n_j.lt.lep1.and.n_j.ge.lcentr.and.n_iz.gt.0.and. &
          n_iz.le.mimxbzf)then
        ilinside_pl=.true.
     end if
     !else end point is outside plasma
     !mg         return
     if(ilinside_pl) return
     rb(1)=zx_in
     rb(2)=zy_in
     rb(3)=zl_in
     call find_point_inside
  end if

contains

  subroutine find_point_inside
    ! rb [cm]
    call box2car(ldep,rb,r)   ! map  x,y,z coordinates of point in box coordinates to tokamak, [m] ->[cm]
    r=r*0.01_fp               !cm -> m
    call eq_rcyl(nsnccw*nsjdotb,r(1),r(2),zr,zphi) ![m,m] -> [m,rad]
    ! find transp zone n_j, and beam zone n_iz
    call nbi_zone_ind(zr,zphi,r(3),n_j,n_iz)
    if(n_j.lt.lep1.and.n_j.ge.lcentr.and.n_iz.gt.0.and. &
         n_iz.le.mimxbzf)then
       ilinside_pl=.true.     !start point is inside plasma
    end if

    if(ilinside_pl.and.zdlseg.gt.0.0_fp) then
       !     start point is inside plasma but end point is outside plasma
       !     find end point that inside plasma
       !     unit vector
       r_u(1)=(zxf-zx)/zdlseg
       r_u(2)=(zyf-zy)/zdlseg
       r_u(3)=(zlf-zl)/zdlseg
       inum=10
       ii=1
       zll=zdlseg/2.0_fp
       rb1(:)=rb(:)
       ilinside_pl_prev=ilinside_pl
       do while(ii.le.inum)
          if(ilinside_pl_prev.neqv.ilinside_pl) then
             if(ii.ne.2) then
                if(.not.ilinside_pl) rb(:)= rb1(:) !restore prev value, which is in plasma
                ilinside_pl=.true.
                exit
             end if
          end if
          ilinside_pl_prev=ilinside_pl !save prev value and find new one for ilinside_pl
          !     rb [cm]
          rb1(:)=rb(:)        !save prev. value
          rb(:)= rb1(:)+zll*r_u(:)
          call box2car(ldep,rb,r) ! map  x,y,z coordinates of point in box coordinates to tokamak, [m] ->[cm]
          r=r*0.01_fp         !cm -> m
          call eq_rcyl(nsnccw*nsjdotb,r(1),r(2),zr,zphi) ![m,m] -> [m,rad]
          !     find transp zone n_j, and beam zone n_iz
          call nbi_zone_ind(zr,zphi,r(3),n_j,n_iz)

          ilinside_pl=.false.
          zll=-abs(zll)/2

          if(n_j.lt.lep1.and.n_j.ge.lcentr.and.n_iz.gt.0.and. &
               n_iz.le.mimxbzf)then
             ilinside_pl=.true.
             zll=-zll
          end if
          ii=ii+1
       end do
       !     if inum itteration has been done and track is still inside plasma
       !     then assume this point as an last point in
       !     plasma region
       zdlseg=sqrt((rb(1)-zx)**2+(rb(2)-zy)**2+(rb(3)-zl)**2)
       zx=rb(1)
       zy=rb(2)
       zl=rb(3)
    else
       rb(1)=zxf
       rb(2)=zyf
       rb(3)=zlf
    end if


  end subroutine find_point_inside

end subroutine boxinc
