!
! this module is used to deposit particles in tubes ("beam in tube") in a similar

! way that depbox is used to deposit particles in a box.  the variables defined in ! this module will be built or rebuilt when bmzini is called on the first invocation
! of nbstart.
!
module tube_mod
  !
  ! explicitly import common variables
  !
  use iso_c_binding, only: fp => c_double
  use nbi_dimensions, only: mitubem, mituben

  use nbi_com, only: nbtube,  nbstube,   nbetube, rtube,   ytube,   xzetatube
  use nbi_com, only: phitube, thetatube, lmidtube, xl1tube, xl2tube, nsegtube
  use nbi_com, only: rhotube, ndeptube
  use nbi_com, only: tuben0, tubedot, tubes0, itube_try
  use nbi_com, only: nonlin, degrad, rsta, ysta, esta, phia, thea, xbzeta
  use nbi_com, only: nbeam, nsjdotb, nsnccw, vion, nlsym2b, plhgt
  use nbi_com, only: ffulla, fhalfa, xninja, nsbeam, abeama, abeams, lsbeama
  use nbi_com, only: xzbeama, xzbeams, lsbeam, nbshapa, bmwidza, bmwidra, xzbeami
  use nbi_com, only: ldep, xpi, nlheavy, nbapsha, xzpedga, rapedga
!mg  use nbi_com, only: xbgrid,ybgrid,dxdlbox,dydlbox,dlbxdtk
  use nbi_com, only: ntrace,nbsbox,nbebox

  use dep_mpi_split
  use logmod

  implicit none
  save

  ! --- state ---
  logical              :: isetup_tube=.false. ! .true. when tube_setup has been called

  real(fp), allocatable  :: loctube(:,:)       ! (x,y,z) location of tube start [3,nbtube] (cm)
  real(fp), allocatable  :: dirtube(:,:)       ! unit vector in direction of tube [3,nbtube] 

  real(fp), allocatable  :: deltaseg(:)        ! length of each segment
  real(fp), allocatable ::  count_int(:)      ! hit count

  integer :: ipccw                           ! +1 for ip ccw, -1 for ip cw

  !
  ! --- mpi buffering ---
  ! buffer organization repeated for each tube
  !   itube_try
  !   s0*itube_try
  !   tuben0(1)*itube_try
  !   tuben0(2)*itube_try
  !    ...
  !   tuben0(nsegtube)*itube_try
  !   tubedot(1)*itube_try
  !   tubedot(2)*itube_try
  !    ...
  !   tubedot(nsegtube)*itube_try
  !   
  real(fp), allocatable :: deptubebuf(:)       ! buffer filled on client processors with client data
  real(fp), allocatable :: deptubebuf1(:)      ! buffer containing reduced data

  ! --- parameters ---
  real(fp), parameter :: nueps = .01d0  ! trajectory considered parallel when the fractional
                                    ! radial change along the entire sectioned length of the tube
                                    ! is less then this number

  real(fp), parameter :: segeps = .02d0 ! fraction of a segment for which a track may be included in
                                    ! an adjacent segment

  real(fp), parameter :: minrho = .01d0 ! minimum allowable tube radius (cm)

  real(fp), parameter :: fatness = .1d0 ! minimum length/radius of semented part of tube

  character*12, parameter :: deptube_xplasma_name = "deptube_data"  ! name of blackbox data in xplasma object

  ! --- transfer of data between procedures ---
  integer :: itube_int              ! tube used in last tube_intersect_nbline
  real(fp)  :: s0_int                 ! tube distance from last tube_intersect_nbline
  real(fp)  :: slope_int              ! slope from last tube_intersect_nbline
  real(fp)  :: t1_int                 ! entrance intercept time from last tube_intersect_nbline
  real(fp)  :: t2_int                 ! exit intercept time from last tube_intersect_nbline
  real(fp)  :: hitdot_int             ! v_hat.b_hat = cos(angle_particle_to_tube)
  integer :: ihit_int               ! hit flag from last tube_intersect_nbline


  ! ---- interfaces -----
  interface tube_addn0
     module procedure tube_addn0_r4, tube_addn0_r8
  end interface

  private allocate_deptube_buf
contains

  !
  ! ------------- allocate_deptube_buf -----------
  ! create the buffers for mpi
  !
  subroutine allocate_deptube_buf(ier)
    integer, intent(out) :: ier    ! nonzero on error

    integer :: iallot   ! size of allocated buffer
    integer :: ineed    ! required size of buffer

    ier=0

    iallot = -1
    if (allocated(deptubebuf)) then
       iallot = size(deptubebuf)
    end if

    ineed = 2*(nbtube + sum(nsegtube(1:nbtube)))
#ifdef __DEBUG
    ineed = 2*nbtube + 8*sum(nsegtube(1:nbtube))
#endif
    if (iallot/=ineed) then
       if (allocated(deptubebuf)) deallocate(deptubebuf,deptubebuf1)
       allocate(deptubebuf(ineed),deptubebuf1(ineed))
       deptubebuf = 0
       deptubebuf1 = 0
    end if
  end subroutine allocate_deptube_buf

  !
  ! ------------- deptube_buf_put -----------
  ! create the buffers for mpi and fill with the current tube data
  !
  subroutine deptube_buf_put(ier)
  
  use nbi_mpi_com
#ifdef __DEBUG
  use nbi_halo_mod_def, only : tuben0_g
#endif

    integer, intent(out) :: ier    ! nonzero on error

    integer :: k     ! index of next free slot in buffer
    integer :: it    ! number of trajectories attempted
    integer :: ks    ! number of segments on tube
    integer :: itube ! loop variable
    integer :: iigen ! loop variable

    call allocate_deptube_buf(ier) ; if (ier/=0) return

    k=1
    do itube=1,nbtube
       it = itube_try(itube)
       ks = nsegtube(itube)
       deptubebuf(k) = it                              ; k=k+1
!mg 09/20/12 i removed it*..., as it is now consistent with deptube
!mg       deptubebuf(k) = tubes0(itube)                ; k=k+1
       deptubebuf(k) = count_int(itube)                ; k=k+1 !only if tubes are used 
                                                               !for beam neutrals in tube
       deptubebuf(k:(k+ks-1)) = tuben0(1:ks,itube)  ; k=k+ks ; 
       deptubebuf(k:(k+ks-1)) = tubedot(1:ks,itube) ; k=k+ks ; 
#ifdef __DEBUG
       do iigen=0,5
          deptubebuf(k:(k+ks-1)) = tuben0_g(1:ks,itube,iigen) ; k=k+ks ; 
       enddo
#endif

    end do

    if (k /= (size(deptubebuf)+1)) then
       call errorlog( 'deptube_buf_put: logic error')
       ier=1
    end if
  end subroutine deptube_buf_put
  
  !
  ! ------------- deptube_buf_get -----------
  ! create the buffers for mpi and fill with the mpi reduced data
  !
  subroutine deptube_buf_get(ier)

  use nbi_mpi_com
#ifdef __DEBUG
  use nbi_halo_mod_def, only : tuben0_g
#endif

    integer, intent(out) :: ier    ! nonzero on error

    integer :: k     ! index of next free slot in buffer
    integer :: it    ! number of trajectories attempted but nonzero
    integer :: ks    ! number of segments on tube
    integer :: itube ! loop variable
    integer :: iigen ! loop variable
!mg    call allocate_deptube_buf(ier) ; if (ier/=0) return
! sanitary check
    if(.not.allocated(deptubebuf1)) then
       call errorlog( 'deptube_buf_get: buffer does not allocated')
       ier = 1
       return
    endif

    k=1
    do itube=1,nbtube
       ks = nsegtube(itube)

       itube_try(itube) = int(deptubebuf1(k)) ; k=k+1
!mg 09/20/12 i removed .../it, as it is now consistent with deptube, deptube_buf_put
!mg       tubes0(itube)       = deptubebuf1(k)          ; k=k+1
       count_int(itube)       = deptubebuf1(k)          ; k=k+1
       tuben0(1:ks,itube)  = deptubebuf1(k:(k+ks-1)) ; k=k+ks ; 
       tubedot(1:ks,itube) = deptubebuf1(k:(k+ks-1)) ; k=k+ks ; 
#ifdef __DEBUG
       do iigen=0,5
          tuben0_g(1:ks,itube,iigen)  = deptubebuf1(k:(k+ks-1)) ; k=k+ks ; 
       enddo
#endif

    end do

    if (k /= (size(deptubebuf1)+1)) then
       call errorlog( 'deptube_buf_get: logic error')
       ier=1
    end if
     if(allocated(deptubebuf)) deallocate(deptubebuf)
     if(allocated(deptubebuf1)) deallocate(deptubebuf1)
   
  end subroutine deptube_buf_get
  
  !
  ! ------------- deptube_xplasma_create -----------
  ! fill xplasma with the current tube data
  !
  subroutine deptube_xplasma_create(s, ier)
    use xplasma_obj_instance, only: xplasma_create_blackbox, xplasma

    type(xplasma), pointer     :: s      ! xplasma object 
    integer,       intent(out) :: ier    ! nonzero on error

    integer :: k     ! index of next free slot in array
    integer :: ks    ! number of segments on tube
    integer :: itube ! loop variable

    integer :: ixt(nbtube+1)  ! index into data aray for each tube + 1 for indexing

    real(fp), allocatable :: zd(:) ! packed data array

    ier=0

    k=1
    do itube=1, nbtube
       ixt(itube) = k
       k = k+2*nsegtube(itube)+1
    end do
    ixt(nbtube+1)=k

    allocate(zd(max(1,k-1))) ! funny netcdf message if use zero size array
    zd = 0.d0

    do itube=1,nbtube
       ks = nsegtube(itube)
       k  = ixt(itube)
       
       zd(k)             = tubes0(itube)
       zd(k+1:k+ks)      = tuben0(1:ks,itube)
       zd(k+ks+1:k+2*ks) = tubedot(1:ks,itube)
    end do

    call xplasma_create_blackbox(s, deptube_xplasma_name, 1000, k, ier, &
         iarray = ixt, r8array = zd, label = "first generation beam data along tubes")

    if (ier/=0) then
       call errorlog('deptube_xplasma_create: error building blackbox')
       ier=1
    end if

100 continue
    deallocate(zd)
  end subroutine deptube_xplasma_create
  

  !
  ! ----------------------- tube_check ------------------------
  ! check tube input data for consistency -- allocation not necessary at this point
  !
  subroutine tube_check(ier)
    integer, intent(out) :: ier   ! error flag

    integer  :: i                 ! loop variable
    integer  :: ibs, ibe          ! beam source and energy index
    integer  :: imaxe             ! maximum energy index
    real(fp)   :: zdtube            ! length of sectioned tube

    ier=0
    if (rhotube<minrho) then
       call errorlog('tube_setup: tube radius is too small, minimum = ', minrho)
       ier=1
    end if
    do i = 1, nbtube
       zdtube = xl2tube(i)-xl1tube(i)

       if (zdtube<fatness*rhotube) then
          call errorlog('tube_setup: tube is too fat or xl1tube>=xl2tube for tube ',i)
          ier=1 ; cycle
       end if

       if (nsegtube(i)<1) then
          call errorlog('tube_setup: number of segments must be greater then 0 for tube ',i)
          ier=1 ; cycle
       end if

       if (nsegtube(i)>mituben) then
          call errorlog('tube_check: tube ', i)
          call errorlog(' has too many segments = ', nsegtube(i))
          call errorlog('      max= ', mituben)
          ier=1 ; cycle
       endif

       ibs = nbstube(i)
       ibe = nbetube(i)
       if (ibs<=0 .or. ibs>nbeam) then
          call errorlog('tube_check: tube ', i)
          call errorlog(' refers to nonexistent beam ', ibs)
          ier=1 ; cycle
       end if

       imaxe=1
       if (abs(xzbeama(ibs)-1.0)<=.05) imaxe=3  ! hydrogen isotope
       if ((ibe<=0) .or. (ibe>imaxe)) then
          call errorlog( ' ?tube_check: energy fraction input is out of bounds for tube ', i)
          call errorlog( '  beam ', ibs)
          ier=1 ; cycle
       end if
    end do

  end subroutine tube_check

  !
  ! ---------------------- tube_allocate ----------------------
  ! allocate local module arrays
  !
  subroutine tube_allocate(ier)
    integer, intent(out) :: ier   ! error flag
    !integer :: is                 ! maximum number of segments

    ier=0

    call tube_check(ier)
    if (ier/=0) return            ! bad_input

    !
    ! allocate based on maximum number of segments and number of tubes
    !
    !is = 0
    !if (nbtube>0) is = maxval(nsegtube(1:nbtube))

    if (allocated(loctube)) deallocate(loctube)
    allocate(loctube(3,nbtube))

    if (allocated(dirtube)) deallocate(dirtube)
    allocate(dirtube(3,nbtube))

    if (allocated(deltaseg)) deallocate(deltaseg)
    allocate(deltaseg(nbtube))

    if (allocated(count_int)) deallocate(count_int)
    allocate(count_int(nbtube))
    count_int=0
  end subroutine tube_allocate

  !
  ! --------------------- tube_setup ------------------------
  ! allocate and setup the local data
  !
  subroutine tube_setup(ier)
    integer, intent(out) :: ier   ! error flag

    integer :: i          ! loop variable
    integer :: ib         ! beam index
    integer :: ihit       ! hit flag

    real(fp)  :: zphi       ! beam sightline corrected for ip direction
    real(fp)  :: zx(3)      ! location of aperture
    real(fp)  :: zvec(3)    ! direction of beam sightline
    real(fp)  :: zdist      ! nearest distance to tube from sightline
    real(fp)  :: ztd        ! distance from aperture to tube nearest distance
    real(fp)  :: zt1        ! if hit, distance to entrance
    real(fp)  :: zt2        ! if hit, distance to exit
    real(fp)  :: zs0,zslope ! s = zs0 + zslope*t, formula constants

    ier = 0
    call infolog('tube_setup: setting up tube data')

    isetup_tube = .true.

    call tube_allocate(ier)
    if (ier/=0) return       ! bad input

    ipccw = nsnccw*nsjdotb   ! ip direction -- phi variables referenced from positive ip

    ! 
    ! find rectangular coordinates of tubes
    !
    do i=1, nbtube
       loctube(1,i) = rtube(i)*cos(degrad*xzetatube(i))  ! a vector
       loctube(2,i) = rtube(i)*sin(degrad*xzetatube(i))
       loctube(3,i) = ytube(i)

       zphi = ipccw*phitube(i)  ! positive direction in positive zeta direction

       dirtube(1,i) = -cos(degrad*thetatube(i))*cos(degrad*(zphi-xzetatube(i)))  ! b-hat vector
       dirtube(2,i) =  cos(degrad*thetatube(i))*sin(degrad*(zphi-xzetatube(i)))
       dirtube(3,i) =  sin(degrad*thetatube(i))

       deltaseg(i) = (xl2tube(i)-xl1tube(i))/nsegtube(i)

       ib = nbstube(i)
       call nbline(ib,1)        ! find coordinates down center of beam at aperture

       tubes0(i)=0.0                 ! avoid undefined error in tube_intersect
       call part_nbline(zx, zvec)
       call tube_intersect(i, zx, zvec, ihit, zdist, ztd, zt1, zt2, zs0, zslope)
       
       if (lmidtube(i)==0) then
          tubes0(i) = xl1tube(i)
       else
          if (ihit==1) then
             call errorlog( '?tube_setup: can not center tube with tube parallel to beam')
             ier=1 ; return
          end if
          tubes0(i) = xl1tube(i) + zs0 + ztd*zslope  ! center wrt to beam line
       end if
     
    end do
  end subroutine tube_setup
  
  !
  ! --------------------- part_nbline ----------------------
  ! returns the location and direction vectors for the particle trajectory produced 
  ! from the last nbline invocation assuming a velocity of 1.0
  !
  subroutine part_nbline(zx, zvec)
    real(fp),  intent(out) :: zx(3)       ! (x,y,z) coordinates of particle start (cm)
    real(fp),  intent(out) :: zvec(3)     ! (vx,vy,vz) velocity coordinates of particle (cm/sec)

    real(fp)  :: zeta     ! toroidal angle of start point in radians
    real(fp)  :: zphia    ! trajectory relative to positive zeta direction in radians
    real(fp)  :: ztheta   ! elevation angle of trajectory in radians

    if (.not. isetup_tube) call not_setup_exit

    zeta = degrad*esta  ! start angle

    zx(1) = rsta*cos(zeta)
    zx(2) = rsta*sin(zeta)
    if (nlsym2b) then
       zx(3) = ysta + plhgt  ! up/down symmetric -> ysta is relative to plasma center
    else
       zx(3) = ysta
    end if

    zphia  = degrad*ipccw*phia
    ztheta = degrad*thea

    zvec(1) = -cos(ztheta)*cos(zphia-zeta)
    zvec(2) =  cos(ztheta)*sin(zphia-zeta)
    zvec(3) =  sin(ztheta)
  end subroutine part_nbline

  !
  ! -------------------- not_setup_exit ----------------------
  !
  subroutine not_setup_exit
    call errorlog('?tube_mod::not_setup_exit: tube data is not setup')
    call nbi_exit    
  end subroutine not_setup_exit

  !
  ! ------------------------ tube_intersect ------------------
  ! evaluate the intersection of a trajectory with a tube
  ! zt1,zt2 are not valid if tube is not hit
  !
  subroutine tube_intersect(itube, zx, zvec, ihit, zdist, ztd, zt1, zt2, zs0, zslope)

    integer, intent(in) :: itube       ! tube index

    real(fp),  intent(in) :: zx(3)       ! (x,y,z) coordinates of particle start (cm)
    real(fp),  intent(in) :: zvec(3)     ! (vx,vy,vz) velocity coordinates of particle (cm/sec)
    
    integer, intent(out) :: ihit       ! 0  -> trajectory does not intersect tube
                                       ! 1  -> trajectory intersects and is parallel to tube
                                       ! 2  -> trajectory intersects and is not parallel to tube
    real(fp),  intent(out) :: zdist      ! closest distance trajectory comes to tube sightline 
    real(fp),  intent(out) :: ztd        ! time of particle at closest distance to tube
    real(fp),  intent(out) :: zt1,zt2    ! for ihit=2, this is the time of entrance and exit
                                       ! from the tube, for ihit=1 this is the time to enter 
                                       ! and exit the ends of the tube
    real(fp),  intent(out) :: zs0,zslope ! provides the linear transformation from t to s
                                       !  s = zs0 + zslope*t

    real(fp)   :: delta_tube             ! segmented length of tube
    real(fp)   :: rel(3)                 ! from tube start to trajectory start
    real(fp)   :: y(3), nu(3),nu2,z1(3)  ! junk
    real(fp)   :: delta(3)               ! vector orthogonally from tube to trajectory
    real(fp)   :: zdist2                 ! zdist**2
    real(fp)   :: dt                     ! time of flight in tube/2.
    real(fp)   :: center                 ! center of tube
    real(fp)   :: zz                     ! temporary

    if (.not. isetup_tube) call not_setup_exit

    zt1=0.0
    zt2=0.0
    delta_tube = xl2tube(itube)-xl1tube(itube)

    rel    = zx-loctube(:,itube)            ! x-a
    zs0    = zzdot(rel,  dirtube(:,itube))  ! (x-a).b_hat
    zslope = zzdot(zvec, dirtube(:,itube))  ! v.b_hat

    y  = rel  - zs0*dirtube(:,itube)        ! y  = (x-a)-b_hat*(x-a).b_hat
    nu = zvec - zslope*dirtube(:,itube)     ! nu = v - b_hat*v.b_hat

    nu2 = zzdot(nu,nu)
    if (sqrt(nu2) < nueps*rhotube*zslope/delta_tube) then
       !
       ! almost parallel, deviation along length of segmented tube
       ! is less then nueps*rho
       !
       !print *, '%tube_intersect: trajectory is parallel to tube'
       center = tubes0(itube)+(delta_tube/2)           ! approximate center of tube in s
       ztd    = (center-zs0)/zslope                 ! not really but good enough
       delta  = y + ztd*nu
       zdist  = sqrt(zzdot(delta,delta))
       if (zdist<=rhotube) then
          ihit=1
          zt1 = (tubes0(itube)-zs0)/zslope
          zt2 = zt1+delta_tube/zslope
          if (zt1>zt2) then
             zz =zt1    ! swap so zt1<=zt2
             zt1=zt2
             zt2=zz
          end if
       else
          ihit=0        ! missed
       end if
       return
    end if

    ztd = - zzdot(y,nu)/nu2
    call zzcross(y,nu,z1)

    zdist2 = zzdot(z1,z1)/nu2
    zdist  = sqrt(zdist2)

    if (rhotube<=zdist) then
       ihit=0           ! missed the tube
    else
       dt = sqrt(max(0.0_fp,(rhotube*rhotube-zdist2)/nu2))
       zt1 = ztd-dt
       zt2 = ztd+dt
       ihit=2
    end if
  end subroutine tube_intersect

  !
  ! --------------------- tube_intersect_nbline ------------------------
  ! find intersection of argument tube with the last nbline trajectory and
  ! store results in local module variables
  !
  subroutine tube_intersect_nbline(itube, ihit)
    integer, intent(in)  :: itube       ! tube index
    integer, intent(out) :: ihit        ! hit flag

    real(fp) :: zx(3), zvec(3)   ! location and velocity vectors for trajectory
    real(fp) :: zdist, ztd       ! misc

    call part_nbline(zx, zvec)         ! get nbline data
    zvec = vion*zvec                   ! convert from unit vector to velocity vector
    call tube_intersect(itube, zx, zvec, ihit, zdist, ztd, t1_int, t2_int, &
         s0_int, slope_int)
    itube_int  = itube
    ihit_int   = ihit
    hitdot_int = slope_int/vion
    !print '(a,i2,a,f8.3)', 'hit = ', ihit, '  dist = ', zdist
    if (ihit/=0) then
       count_int(itube) = count_int(itube)+1.0
       !print '(a,e18.6,a,e18.6)', 't1=',t1_int, '  t2=',t2_int
    end if
  end subroutine tube_intersect_nbline

  !
  ! ---------------------- findseg -------------------------
  ! return the segment index containing the segment length argument
  !
  integer function findseg(itube, s)
    integer, intent(in) :: itube  ! tube indes
    real(fp),  intent(in) :: s      ! length from start point
    findseg = 1+int((s-tubes0(itube))/deltaseg(itube))
  end function findseg


  !
  ! ------------------ pdelt -------------------
  ! returns the average amount of time a particle
  ! exists in a time interval delt given a decay rate of kappa
  !
  real(fp) function pdelt(delt,kappa,decay)
    real(fp), intent(in)  :: delt    ! advancement time
    real(fp), intent(in)  :: kappa   ! decay factor
    real(fp), intent(out) :: decay   ! exp(-delt*kappa), ugly
    real(fp) :: x,z

    x = kappa*delt
    if (abs(x)<1.d-4) then
       z = (1.0d0-0.5d0*x*(1.-x/3.0d0))
       pdelt = delt*z
       decay = 1.-x*z
    else if (x>40) then
       pdelt = 1./kappa
       decay = 0.0d0
    else
       decay = exp(-x)
       pdelt = (1.-decay)/kappa
    end if
  end function pdelt

  !
  ! ---------------------- tube_addn0_r8 ------------------------
  ! add to tuben0(i,itube_int) the amount of time a particle spends in
  ! each tube segment from time index tstart to tstart+deltat
  !
  subroutine tube_addn0_from_box_r8(idtrk,wght)
    !tstar=t1_int
    !tend=t2_int
    use nbi_halo_mod_def

    implicit none	
    integer, intent(in) :: idtrk     ! neutral track id#
    real(fp), intent(in) :: wght !neutral track weight
    integer :: idir      ! increment to get to next segment
    integer :: nseg      ! number of sements
    real(fp)  :: seg0      ! distance to first segment
    real(fp)  :: dseg      ! segment size
    real(fp)  :: dsmin     ! minimum increment along tube
    real(fp)  :: tend      ! end of time interval to add


    real(fp)  :: t1,t2, t2_end,t_zone,t_zone0    ! increasing start point
    real(fp)  :: tfrc
    real(fp)  :: p1,p2 ,pend,p2_end       ! survival probability at t1,t2,pend
    real(fp)  :: dt        ! amount of time left until the end is reached

    integer :: i1        ! segment index containing t1
    integer :: ie        ! segment index containing tstart+deltat where tstart isbeginning of time index
    real(fp)  :: ds        ! increment along tube
    integer :: inext     ! next segment along trajectory
    real(fp)  :: sb,tb     ! tube boundary in s and t
    real(fp)  :: dtb       ! time increment to get to tube boundary
    real(fp)  :: x,z       ! decay factor and expected time in segment
    
    integer :: i_zone_in, i_zone_out,i_zone1, i_zone2    ! zone index for a given idtrk
    integer :: iz,intrcpt,i_zone0
    real(fp), allocatable, dimension(:) :: pdep_wrk
    real(fp) :: z_kappa,x_wrk
    integer :: id_gen

    if (.not. isetup_tube) call not_setup_exit
    id_gen=itrackh(idtrk)%idgen

    if (ihit_int==0) return                         ! did not hit tube

    intrcpt=itrackh(idtrk)%maxpts ! # of box zones intercept
    if(allocated(pdep_wrk)) deallocate(pdep_wrk)
    allocate(pdep_wrk(intrcpt))

    pdep_wrk(1:intrcpt)=itrackh(idtrk)%pdep(1:intrcpt)
!    wght=itrackh(idtrk)%weight*ftarg
    idir  = int(sign(1.0_fp,slope_int)) 
    nseg  = nsegtube(itube_int)
    seg0  = tubes0(itube_int)           
    dseg  = deltaseg(itube_int)     
    dsmin = segeps*dseg

    t1 = t1_int
    i_zone_in=itrk_zone(idtrk,t1_int,p1)
    i_zone_out=itrk_zone(idtrk,t2_int,pend)
    i_zone1=i_zone_in
    i_zone2=i_zone_out

    if(i_zone_in.eq.0) then
       call infolog('tube_addn0_from_box_r8: track did not hit tube')
       return
    endif
    tend = t2_int

    dt = tend-t1
    i1 = findseg(itube_int,s0_int+slope_int*t1)      ! segment where this point is located
    ie = findseg(itube_int,s0_int+slope_int*tend)    ! segment containing last time point

    if (i1<1    .and. ie<1)    return   ! interval entirely before tube segments
    if (i1>nseg .and. ie>nseg) return   ! interval entirely after tube segments

    do while (dt>0.0 .and. p1>1.d-20)
       if ((idir>0 .and. i1>nseg) .or. (idir<0 .and. i1<1)) exit  ! past the end
       ds = abs(slope_int*dt)           ! increment in s if taking the full dt step
       if (ds<dsmin) then
          ! add to one segment only
          if (i1>=1 .and. i1<=nseg) then
             p2=max(pend,itrackh(idtrk)%pdep(i_zone1))
             t_zone=itrackh(idtrk)%xlngtk(i_zone1)/itrackh(idtrk)%vneut
             i_zone0=i_zone1-1
             if(i_zone0.eq.0) then
                t_zone0=0.0_fp
             else
                t_zone0=itrackh(idtrk)%xlngtk(i_zone0)/itrackh(idtrk)%vneut
             endif
             t2=min(tend,t_zone)
             z_kappa=itrackh(idtrk)%dkappa(i_zone1)
             tfrc=pdelt(t2-t1,z_kappa,x_wrk)
             tuben0(i1,itube_int)  = tuben0(i1,itube_int) +p2*tfrc*wght
             tubedot(i1,itube_int) = tubedot(i1,itube_int)+p2*tfrc*hitdot_int*wght
#ifdef __DEBUG
             if(id_gen.le.5) then
                tuben0_g(i1,itube_int,id_gen)  = tuben0_g(i1,itube_int,id_gen) +p2*tfrc*wght
             endif
#endif
             t1=t2
             p1=p2
             if(t2.lt.tend) then !more to add
                do iz=i_zone1,i_zone2-1
                   t_zone0=t_zone
                   p2=max(pend,itrackh(idtrk)%pdep(iz+1))
                   t_zone=itrackh(idtrk)%xlngtk(iz+1)/itrackh(idtrk)%vneut
                   t2=min(tend,t_zone)
                   z_kappa=itrackh(idtrk)%dkappa(iz+1)
                   tfrc=pdelt(t2-t1,z_kappa,x_wrk)
                   tuben0(i1,itube_int)  = tuben0(i1,itube_int) +p2*tfrc*wght
                   tubedot(i1,itube_int) = tubedot(i1,itube_int)+p2*tfrc*hitdot_int*wght
#ifdef __DEBUG
                   if(id_gen.le.5) then
                      tuben0_g(i1,itube_int,id_gen)  = tuben0_g(i1,itube_int,id_gen) +p2*tfrc*wght
                   endif
#endif
                   t1=t2
                   p1=p2
                enddo
             endif
          end if
          dt=0.0   ! done
       else
          inext = i1+idir                   ! next segment in direction of increasing t
          do while ((idir>0 .and. inext<1) .or. (idir<0 .and. inext>nseg))
             inext = inext+idir             ! skip to end of the tube if heading toward segments
          end do
          sb = seg0 + dseg*(2*inext-idir-1)/2   ! next tube boundary in s, boundary of inext,inext-idir
          tb = (sb-s0_int)/slope_int            ! tube boundary in time
          dtb= max(0.0_fp,tb-t1)                ! distance in time to next tube boundary
          dtb = min(dtb,dt)                     ! step to boundary or end of time interval
          t2 = t1+dtb
          t2_end = t2
          i_zone2=itrk_zone(idtrk,t2_end,p2_end)
          if (i1>=1 .and. i1<=nseg) then
             p2=max(p2_end,itrackh(idtrk)%pdep(i_zone1))
             !find box zone length in time 
             t_zone=itrackh(idtrk)%xlngtk(i_zone1)/itrackh(idtrk)%vneut
             i_zone0=i_zone1-1
             if(i_zone0.eq.0) then
                t_zone0=0.0_fp
             else
                t_zone0=itrackh(idtrk)%xlngtk(i_zone0)/itrackh(idtrk)%vneut
             endif
             t2=min(t2_end,t_zone)
             z_kappa=itrackh(idtrk)%dkappa(i_zone1)
             tfrc=pdelt(t2-t1,z_kappa,x_wrk)
             tuben0(i1,itube_int)  = tuben0(i1,itube_int) +p2*tfrc*wght
             tubedot(i1,itube_int) = tubedot(i1,itube_int)+p2*tfrc*hitdot_int*wght
#ifdef __DEBUG
             if(id_gen.le.5) then
                tuben0_g(i1,itube_int,id_gen)  = tuben0_g(i1,itube_int,id_gen) +p2*tfrc*wght
             endif
#endif
             t1=t2
             p1=p2
             if(t2.lt.t2_end) then !more to add
                do iz=i_zone1,i_zone2-1
                   t_zone0=t_zone
                   p2=max(p2_end,itrackh(idtrk)%pdep(iz+1))
                   t_zone=itrackh(idtrk)%xlngtk(iz+1)/itrackh(idtrk)%vneut
                   t2=min(t2_end,t_zone)
                   z_kappa=itrackh(idtrk)%dkappa(iz+1)
                   tfrc=pdelt(t2-t1,z_kappa,x_wrk)
                   tuben0(i1,itube_int)  = tuben0(i1,itube_int) +p2*tfrc*wght
                   tubedot(i1,itube_int) = tubedot(i1,itube_int)+p2*tfrc*hitdot_int*wght
#ifdef __DEBUG
                   if(id_gen.le.5) then
                      tuben0_g(i1,itube_int,id_gen)  = tuben0_g(i1,itube_int,id_gen) +p2*tfrc*wght
                   endif
#endif
                  t1=t2
                  p1=p2
                enddo
             endif
          end if
          p1 = p2_end
          dt = dt-dtb
          i1 = inext     ! if boundary point was not used then dt==0.0 and loop will exit
          i_zone1=i_zone2
       end if
    end do
    
  end subroutine tube_addn0_from_box_r8

  !
  ! ---------------------- tube_addn0_r8 ------------------------
  ! add to tuben0(i,itube_int) the amount of time a particle spends in
  ! each tube segment from time index tstart to tstart+deltat
  !
  subroutine tube_addn0_r8(tstart, deltat, pstart, kappa)
    real(fp), intent(in) :: tstart    ! beginning of time index
    real(fp), intent(in) :: deltat    ! time increment
    real(fp), intent(in) :: pstart    ! survival probability at tstart
    real(fp), intent(in) :: kappa     ! decay rate

    integer :: idir      ! increment to get to next segment
    integer :: nseg      ! number of sements
    real(fp)  :: seg0      ! distance to first segment
    real(fp)  :: dseg      ! segment size
    real(fp)  :: dsmin     ! minimum increment along tube
    real(fp)  :: tend      ! end of time interval to add

    real(fp)  :: t1        ! increasing start point
    real(fp)  :: p1        ! survival probability at t1
    real(fp)  :: dt        ! amount of time left until the end is reached

    integer :: i1        ! segment index containing t1
    integer :: ie        ! segment index containing tstart+deltat
    real(fp)  :: ds        ! increment along tube
    integer :: inext     ! next segment along trajectory
    real(fp)  :: sb,tb     ! tube boundary in s and t
    real(fp)  :: dtb       ! time increment to get to tube boundary
    real(fp)  :: x,z       ! decay factor and expected time in segment

    if (.not. isetup_tube) call not_setup_exit

    !
    ! restrict time range to tube intercept time ranges
    !
    tend  = tstart+deltat

    if (ihit_int==0) return                         ! did not hit tube
    if (tstart>=t2_int .or. tend<=t1_int) return    ! outside of intercept range

    idir  = int(sign(1.0_fp,slope_int)) 
    nseg  = nsegtube(itube_int)
    seg0  = tubes0(itube_int)           
    dseg  = deltaseg(itube_int)     
    dsmin = segeps*dseg

    if (tstart<t1_int) then
       x = kappa*(t1_int-tstart)  ! decay
       if (x>40) return           ! decayed before it got to the tube
       p1 = pstart*exp(-x)        ! set start to first intercept
       t1 = t1_int
    else
       t1 = tstart
       p1 = pstart
    end if
    tend = min(tend, t2_int)

    dt = tend-t1
    i1 = findseg(itube_int,s0_int+slope_int*t1)      ! segment where this point is located
    ie = findseg(itube_int,s0_int+slope_int*tend)    ! segment containing last time point

    if (i1<1    .and. ie<1)    return   ! interval entirely before tube segments
    if (i1>nseg .and. ie>nseg) return   ! interval entirely after tube segments

    do while (dt>0.0 .and. p1>1.d-20)
       if ((idir>0 .and. i1>nseg) .or. (idir<0 .and. i1<1)) exit  ! past the end
       ds = abs(slope_int*dt)           ! increment in s if taking the full dt step
       if (ds<dsmin) then
          ! add to one segment only
          if (i1>=1 .and. i1<=nseg) then
             tuben0(i1,itube_int)  = tuben0(i1,itube_int) +p1*pdelt(dt,kappa,x)
             tubedot(i1,itube_int) = tubedot(i1,itube_int)+p1*pdelt(dt,kappa,x)*hitdot_int
          end if
          dt=0.0   ! done
       else
          inext = i1+idir                   ! next segment in direction of increasing t
          do while ((idir>0 .and. inext<1) .or. (idir<0 .and. inext>nseg))
             inext = inext+idir             ! skip to end of the tube if heading toward segments
          end do
          sb = seg0 + dseg*(2*inext-idir-1)/2   ! next tube boundary in s, boundary of inext,inext-idir
          tb = (sb-s0_int)/slope_int            ! tube boundary in time
          dtb= max(0.0_fp,tb-t1)                ! distance in time to next tube boundary
          dtb = min(dtb,dt)                     ! step to boundary or end of time interval
          z = pdelt(dtb,kappa,x)                ! expected time in segment and decay factor
          if (i1>=1 .and. i1<=nseg) then
             tuben0(i1,itube_int)  = tuben0(i1,itube_int) +p1*z
             tubedot(i1,itube_int) = tubedot(i1,itube_int)+p1*z*hitdot_int
          end if
          p1 = p1*x
          t1 = t1+dtb
          dt = dt-dtb
          i1 = inext     ! if boundary point was not used then dt==0.0 and loop will exit
       end if
    end do
    
  end subroutine tube_addn0_r8

  !
  ! ---------------------- tube_addn0_r4 ---------------------
  !
  subroutine tube_addn0_r4(tstart, deltat, pstart, kappa)
    real, intent(in) :: tstart    ! beginning of time index
    real, intent(in) :: deltat    ! time increment
    real, intent(in) :: pstart    ! survival probability at tstart
    real, intent(in) :: kappa     ! decay rate

    real(fp) :: tstart_r8    ! beginning of time index
    real(fp) :: deltat_r8    ! time increment
    real(fp) :: pstart_r8    ! survival probability at tstart
    real(fp) :: kappa_r8     ! decay rate
    tstart_r8 = tstart
    deltat_r8 = deltat
    pstart_r8 = pstart
    kappa_r8  = kappa

    call tube_addn0_r8(tstart_r8, deltat_r8, pstart_r8, kappa_r8)
  end subroutine tube_addn0_r4

  !
  ! ------------------ deptube_from_box -----------------
  ! main routine for looping over tubes and depositing particles
  ! to calculate the average neutral density, tuben0, in each tube
  ! using neutral density from beam_in_box
  !
  subroutine deptube_from_box(ierr)

    use nbi_mpi_com
    use nbi_halo_mod_def

    implicit none	
    integer,intent(out) :: ierr   ! status code: 0=ok
    integer :: i,j,itube,id_box           ! misc
    integer :: isize_trk
    real(fp), dimension(3) :: r_start,r_end,v_vec
    real(fp) :: v_mod
    integer :: ihit,mpierr
    integer, dimension(nbtube) :: itrk_fast,itrk_therm,iwrk
    real(fp) :: zdist, ztd,t1_box,t2_box       ! misc
    real(fp)  :: scpu,fcpu    ! cpu times
    real(fp) :: dvtube,wght
    integer :: intrcpt
#ifdef __DEBUG
    integer :: it1, it2
#endif

    call enterlog("deptube_from_box") ! all returns through exitlog()
    ierr=0

    if (.not. isetup_tube) then
       call tube_setup(j)
       if (j/=0) then
          call errorlog('?deptube: error setting up tube data')
          ierr=1
          return
       end if
    end if
    
    call cpu_time(scpu)
    tuben0(:,1:nbtube)  = 0.0  ! start from scratch
    tubedot(:,1:nbtube) = 0.0
    itube_try(1:nbtube) = 0

#ifdef __DEBUG
    it1=size(tuben0,dim=1)
    it2=size(tuben0,dim=2)
    if(allocated(tuben0_g)) deallocate(tuben0_g)
    allocate(tuben0_g(it1,it2,0:5))
    tuben0_g=0.0_fp 
#endif

    isize_trk=size(itrackh_rqinit) !how many tracks
    itrk_fast=0
    itrk_therm=0
    do i=1,isize_trk !loop over all neutral tracks
       if(itrackh_rqinit(i)) cycle !track hasn't been initialized
       id_box=itrackh(i)%idbox
       !map start point of neutral track from box cartesian to tokamak cartesian 
       call box2car(id_box,itrackh(i)%pta*100.0_fp,r_start) !start point [cm]

       call box2car(id_box,itrackh(i)%ptb*100.0_fp,r_end) !end point[cm]
       v_vec=r_end-r_start
       v_mod=sqrt(zzdot(v_vec,v_vec))
       if(v_mod.gt.0.0_fp)then
          v_vec=v_vec/v_mod !unit velocity vector
          else
             cycle
       endif
       if(itrackh(i)%vneut.le.0.0_fp) cycle
       
       v_vec=v_vec*itrackh(i)%vneut*100.0_fp ![m/sec] to [cm/sec]

       do itube=1,nbtube !loop over all tubes
          if(nbsbox(id_box).ne.nbstube(itube).or.nbebox(id_box).ne.nbetube(itube)) cycle !this track is coming from a 
                                                                                          !box that is not a target 
                                                                                          !of this tube
          call tube_intersect(itube,r_start,v_vec,ihit,zdist,ztd,t1_int,t2_int,s0_int,slope_int)
          if(ihit.eq.0) cycle ! neutral track did not hit tube
          intrcpt=itrackh(i)%maxpts
          itube_int  = itube
          ihit_int   = ihit
          hitdot_int = slope_int/(itrackh(i)%vneut*100.0_fp)
          t1_box =0.0_fp
          t2_box=itrackh(i)%xlngtk(intrcpt)/itrackh(i)%vneut
          if(t1_int.gt.t2_box.or.t2_int.lt.t1_box) cycle !this track started/stopped after/before reaching the tube
          t1_int=max(t1_int,t1_box)
          t2_int=min(t2_int,t2_box)
          wght=itrackh(i)%weight
          call tube_addn0_from_box_r8(i,wght)
          itube_try(itube)=itube_try(itube)+1
          if(itrackh(i)%type.eq.'FASTB')then
             itrk_fast(itube)=itrk_fast(itube)+1
          elseif(itrackh(i)%type.eq.'THERM')then
             itrk_therm(itube)=itrk_therm(itube)+1
          endif
       enddo


    enddo            !end of loop over all neutral tracks

        !
        !  renormalize the monte carlo sum (tuben0)
        !  sum currently contains sum of flight times per tube element (sec),
        !  summed over ntry mc neutral tracks
        !
#ifdef __MPI
    mpierr=0
    iwrk=0
    call mpi_allreduce(itrk_fast,iwrk,nbtube,MPI_INTEGER,MPI_SUM,nbi_comm,mpierr)
    itrk_fast=iwrk
    iwrk=0
    call mpi_allreduce(itrk_therm,iwrk,nbtube,MPI_INTEGER,MPI_SUM,nbi_comm,mpierr)
    itrk_therm=iwrk
    call deptube_buf_reduce
#endif
    do itube=1,nbtube !loop over all tubes
       dvtube = xpi*(rhotube**2)*deltaseg(itube)
       j = nsegtube(itube)
       where(tuben0(1:j,itube)>0.0)
          tubedot(1:j,itube) = tubedot(1:j,itube)/tuben0(1:j,itube)
          tuben0(1:j,itube)  = tuben0(1:j,itube)/dvtube
       end where
#ifdef __DEBUG
       tuben0_g(1:j,itube,0:5)  = tuben0_g(1:j,itube,0:5)/dvtube
#endif
       if(myidd.eq.0) then
!          call errorlog(' not an error tube_mod.deptube_from_box: # of tracks',itube_try(itube))
!          call errorlog(' not an error                          : hit tube = ',itube)
          call errorlog(' not an error tube_mod.deptube_from_box: tube #',itube)
          call errorlog(' not an error                          : fastb hits = ',itrk_fast(itube))
          call errorlog(' not an error                          : therm hits = ',itrk_therm(itube))
       endif
    enddo
#ifdef __DEBUG
    if(allocated(tuben0_g)) deallocate(tuben0_g)
#endif

    call cpu_time(fcpu)
    if (nbtube>0) call infolog('tube_mod.deptube_from_box : cpu time(sec) = ', (fcpu-scpu))
    call exitlog("deptube_from_box") ! all returns through exitlog()
  end subroutine deptube_from_box  
  ! ------------------ deptube -----------------
  ! main routine for looping over tubes and depositing particles
  ! to calculate the average neutral density, tuben0, in each tube
  !
  subroutine deptube(ierr)

    use nbi_mpi_com
    use nbi_random
    implicit none	
    integer :: itube        ! tube index
    integer :: ib           ! beam index
    integer :: ntry         ! number of particles to deposit
    integer :: isb          ! beam species loop
    integer :: idep0        ! storage of ndeptube
    integer :: idep         ! loop variable over particle deposition
    integer :: j            ! misc
    integer :: iscreen1,ind1,iseed,idum ! mpi related

    real(fp)  :: area         ! area of beam
    real(fp)  :: zfrac        ! fraction of beam at selected energy
    real(fp)  :: zinj         ! flux
    real(fp)  :: dvtube       ! voluem of a tube segment
    real(fp)  :: znorm        ! misc
    real(fp)  :: scpu,fcpu    ! cpu times

    integer,intent(out) :: ierr   ! status code: 0=ok
    integer :: nusey
 
    ierr=0

    if (.not. isetup_tube) then
       call tube_setup(j)
       if (j/=0) then
         call errorlog('?deptube: error setting up tube data')
          ierr=1
          return
       end if
    end if

    call cpu_time(scpu)
    tuben0(:,1:nbtube)  = 0.0  ! start from scratch
    tubedot(:,1:nbtube) = 0.0
    itube_try(1:nbtube) = 0

    do itube = 1, nbtube
       ldep = itube           ! pass tube index into dep()

       !
       !  get fraction of flux at selected energy
       !
       ib=nbstube(ldep)
       if(nbetube(ldep).eq.1) then
          zfrac=ffulla(ib)
       else if(nbetube(ldep).eq.2) then
          zfrac=fhalfa(ib)
       else
          zfrac=1.0-ffulla(ib)-fhalfa(ib)
       endif

       !
       ! beam area
       !
       if (nbshapa(ib)==1) then
          area = 4.0_fp*bmwidza(ib)*bmwidra(ib)  ! rectangular beam
       else
          area = xpi*bmwidra(ib)**2          ! circular beam
       end if
       if (nbapsha(ib)==1) then
          area = min(area, 4.0_fp*xzpedga(ib)*rapedga(ib))  ! rectangular beam
       else
          area = min(area,1.0_fp*xpi*rapedga(ib)**2)        ! circular beam
       end if
       ntry = abs(ndeptube)*sqrt(area)/(2*rhotube)

       !
       !  the net flux at the selected energy:
       !
       zinj=xninja(ib)*zfrac   ! #/sec
       if(zinj.le.0.0) cycle   ! no particles so go to next tube

       !
       !  ok there is some flux -- find the species
       !
        do isb=1,nsbeam
           if(lsbeama(ib).eq.isb) then
              lsbeam=isb
              go to 10
           endif
        enddo

! cannot pass this way -- nbstube must pt to a beam
        call errorlog(' ?depbox:  nbstube(ldep)=',ib)
        call errorlog(' ldep = ',ldep)
        call errorlog('   does not match any beam in the simulation.')
        ierr=1
        return

10      continue
        
        call inispc
        
        idep0    =  ndeptube
        ndeptube = -ndeptube  ! this signals dep that calls are from deptube

        call depall_mpi_split(1,ntry,nusey)

        iseed=abs(12111000*nbrandm(idum))

        count_int=0
!       do idep=1,ntry
        do idep=1,nusey

           ind1=ninddumn(idep)
           call nbi_random_init_dep(iseed,ind1,iscreen1)

           itube_try(itube) = itube_try(itube)+1
           call nbidep(0,ierr)
           if(ierr.ne.0) then
              call errorlog('deptube: nbidep error')
              return
           endif
        enddo

        call infolog('tube_mod.deptube: tube ',itube)
        call infolog(' hit count=', count_int(itube))
        call infolog(' desired=', idep0)
        call infolog(' attempt=', ntry)

        ndeptube=idep0     ! restore ndeptube

        !
        !  renormalize the monte carlo sum (tuben0)
        !  sum currently contains sum of flight times per tube element (sec),
        !  summed over ntry mc neutral tracks
        !
        dvtube = xpi*(rhotube**2)*deltaseg(itube)
        znorm=zinj/dvtube         ! n*cm**-3*sec**-1, shouldn't be normalized by *itube_try(itube) here, as it wasn't
                                  ! reduced yet for mpi case
        j = nsegtube(itube)
        tubedot(1:j,itube) = znorm*tubedot(1:j,itube)
        tuben0(1:j,itube)  = znorm*tuben0(1:j,itube)
     end do

#ifdef __MPI
    call deptube_buf_reduce
#endif
    ! comment - itube_try(itube) 
    ! has been already reduced for mpi case
    do itube = 1, nbtube
       if (itube_try(itube).le.0) cycle
       where(tuben0(1:j,itube)>0.0)
          tubedot(1:j,itube) = tubedot(1:j,itube)/tuben0(1:j,itube)
          tuben0(1:j,itube)  = tuben0(1:j,itube)/itube_try(itube)
       end where
    enddo
    call cpu_time(fcpu)
    if (nbtube>0) call infolog('tube_mod.deptube: cpu time(sec) = ', (fcpu-scpu))

    !call boxtube_print(1,1,1)
  end subroutine deptube
  !
  ! ------------------ zzdot -------------------
  ! dot product
  !
  real(fp) function zzdot(a,b)
    real(fp), intent(in) :: a(3),b(3)
    zzdot =  a(1)*b(1)+a(2)*b(2)+a(3)*b(3)
  end function zzdot

  !
  ! ------------------ zzcross -----------------
  ! cross product  c = a x b
  !
  subroutine zzcross(a,b,c)
    real(fp), intent(in)  :: a(3),b(3)
    real(fp), intent(out) :: c(3)
    c(1) = a(2)*b(3)-b(2)*a(3)
    c(2) = a(3)*b(1)-b(3)*a(1)
    c(3) = a(1)*b(2)-b(1)*a(2)
  end subroutine zzcross

  !
  ! --------------- boxcheck -----------------
  ! debugging routine for comparing with a box
  !
  subroutine boxcheck(ix,r,den)
    use nbi_com, only: xybapa, xybsca, xlbapa, xlbtna, rtcena
    use nbi_com, only: nbsbox,nlco, boxn0, nadrbox, nxbox, nybox

    integer, intent(in)  :: ix        ! box index
    real(fp),  intent(in)  :: r(3)      ! coordinates of point in tokamak
    real(fp),  intent(out) :: den       ! box density at this point

    integer :: ib,ipccw
    real(fp)  :: thetal, b, phil
    real(fp)  :: zz, b0(3), that(3), xhat(3), yhat(3)
    real(fp)  :: d(3)
    real(fp)  :: l, x, y

    integer :: jx,jy,jz
    logical :: inside

    ib    = nbsbox(ix)        ! beam index for box
    if(ntrace(ib).ne.0) ib = abs(ntrace(ib))

    ipccw = nsnccw*nsjdotb  

    thetal = asin((xybapa(ib)-xybsca(ib))/xlbapa(ib))  ! sightline elevation
    b      = xlbtna(ib)*cos(thetal)         ! equatorial distance from source to tangent
                      
    phil   = ipccw*atan(rtcena(ib)/b)       ! angle from source to axis assuming positive ip
    
    if (.not.nlco(ib)) phil=-phil           ! now relative to positive zeta

    zz = xbzeta(ib)*degrad                  ! zeta of source
    !mg comment: here b supposed to be sqrt(b*b+rtcena(ib)**2)
    b0(1) = b*cos(zz)                       ! x,y,z of source
    b0(2) = b*sin(zz)
    b0(3) = xybsca(ib)

    that(1) = -cos(thetal)*cos(phil-zz)     ! sightline unit vector
    that(2) =  cos(thetal)*sin(phil-zz)
    that(3) =  sin(thetal)

    call zzcross((/0.0d0,0.0d0,1.d0/),that,xhat)  ! as defined in nbline as rhat direction
    if (phil<0.0) xhat=-xhat                     ! point toward the axis
    xhat = xhat/sqrt(zzdot(xhat,xhat))          ! make unit vector
    call zzcross(that,xhat,yhat)                ! up like direction
    if (yhat(3)<0.0) yhat=-yhat

    d = r-b0            ! vector relative to source
    l = zzdot(d,that)   ! along sightline
    x = zzdot(d,xhat)   
    y = zzdot(d,yhat)

    call boxfind(x,y,l,jx,jy,jz,inside)
    if (inside) then
       den = boxn0(nadrbox(ix)+(jx-1)+nxbox*(jy-1)+nxbox*nybox*(jz-1))
    else
       den = 0.0d0
    end if
  end subroutine boxcheck

  !
  ! -------------------- boxtube_print ---------------------
  ! compares the tube density with the box density
  !
  subroutine boxtube_print(itube, ibox, istop) 
    integer, intent(in) :: itube    ! tube to loop over
    integer, intent(in) :: ibox     ! compare with this box
    integer, intent(in) :: istop    ! nonzero to stop
    integer :: i, io
    real(fp)  :: s, r(3), dbox

    if (.not. isetup_tube) call not_setup_exit
    call find_io_unit(io)

    do i = 1, nsegtube(itube)
       s = tubes0(itube) + (i-0.5d0)*deltaseg(itube)        ! center of segment
       r = loctube(1:3,itube)+s*dirtube(1:3,itube)    ! x,y,z coordinates of point in tokamak
       call boxcheck(ibox,r,dbox)                     ! get box density

       write(io,'(i5,2x,f12.5,2x,e16.8,2x,e16.8)') i, s, dbox, tuben0(i,itube)
    end do
    if (istop /= 0) then
       close(io)
       call infolog('boxtube_print: exiting, istop = ',istop)
       call good_exit
    end if
  end subroutine boxtube_print
  !  
  ! function itrk_zone returns zone# for a given 3d halo neutral track with 'id' and time 't'
  !
  integer function itrk_zone(id,t,p)
    use nbi_halo_mod_def
    implicit none	
    integer, intent(in) :: id !neutral track id#
    real(fp), intent(in) :: t ! time
    real(fp), intent(out) :: p !survival probabily*decay
    real(fp) :: t1,t2,t1_box,t2_box
    integer :: i,intrcpt

    itrk_zone=0
    p=0.0_fp
    intrcpt=itrackh(id)%maxpts
    t1_box =0.0_fp
    t2_box=itrackh(id)%xlngtk(intrcpt)/itrackh(id)%vneut
    if(t.gt.t2_box.or.t.lt.t1_box)then
       call infolog('function itrk_zone: given time is out of track time range')
       return
    endif
    t1=0.0_fp
    t2=0.0_fp
    do i=1,intrcpt
       itrk_zone= itrk_zone+1
       t1=t2 
       t2=itrackh(id)%xlngtk(i)/itrackh(id)%vneut
       if(t.lt.t2) exit
    enddo
    if(itrk_zone.gt.0.and.itrk_zone.le.intrcpt) then
       p=itrackh(id)%pdep(itrk_zone)
    endif
    
  end function itrk_zone
  
end module tube_mod
