!-----------------------------------------------------------------------
!  addn0 -- add to the 3d box fast neutral density monte carlo sum
!
subroutine addn0(zlen,zprob,zmfpi,vneut,ideptrk,lsbeam_parent)
  !
  ! note, by mg 01/15/15 -- old subroutines (boxfind, boxinc) for  beam-in-box
  ! are using left-handed coordinate system. 3d halo model is using
  ! right-handed coordinate system. all related corrections (sign for an x-coordinate)
  ! are taken into consideration. also pboxn0 post-processing program is fixed to
  ! show the result in right-handed coordinate system.
  !
  !  input arguments:
  !    zlen -- how far the neutral has traveled
  !    zprog -- probability that an injected neutral would get this far
  !
  !    if zlen=0.0:  initialize
  !
  !    subsequently, expect a sequence of calls with strictly increasing
  !    zlen and flat or decreasing zprob
  !
  !
  !  using the common outputs from nbline, this routine computes the
  !  intersection of a neutral trajectory with cartesian elements of
  !  a 3d grid (box), and increments the intersected element sums by
  !  an amount proportional to zprob * (intersection length) / vion
  !  i.e. probability of reaching zone * time of flight in zone.
  !
  !-----------------------------------------------------------------------
  !
  use nbi_com
  use nbi_halo_mod_def
  use logmod
  !
  !-----------------------------------------------------------------------
  !  saved local quantities:
  !
  !============
  implicit none
  integer, intent(in) :: lsbeam_parent !will be used only for levmod_halo=2
  !to gather sorbn0_depneed to know
  !lsbeam from 'parent'
  integer :: igridx,igridy,igridl,ix,iy,il,iadr,iadra
  integer :: iadrah ! address for halo (first generation) in a boxn0
  integer :: iadrah1 ! address for halo (>first generation) in a boxn0
  !============
  real(fp) :: zprob,zlen,zlbox,zxbox,zybox,zx,zy,zl,zdl
  !============
  real(fp) :: zlenp  ! trajectory length previously reached
  real(fp) :: zprobp ! previously passed survival probability
  real(fp) :: zmfpi  ! inverse mean free path
  real(fp) :: vneut  ! velocity of neutral,[cm/sec]
  integer :: ideptrk !if =1 then track will be used for a 3d halo, call from dep
  !if =0 then add neutral to a box for a density only, call from dep
  !if >=2 then add neutral to a box for a 3d halo, call from depbox_all

  integer, save :: idtr !for ideptrk >= 2, idtr = ideptrk-2
  !for ideptrk =1, idtr= i_trk0 from nbi_com
  integer, save :: idgen !generation id for neutrals
  !!for ideptrk =< 1, idgen = 0
  integer :: jba_0 ! pitch angle index
  real(fp), save :: xz,xa !z and amu of neutral
  real(fp) :: zlap   ! l at beam aperture
  real(fp) :: zxap   ! x at beam aperture
  real(fp) :: zyap   ! y at beam aperture
  real(fp) :: zxboxp  ! previous interval endpoint coordinate: x
  real(fp) :: zyboxp  ! previous interval endpoint coordinate: y
  real(fp) :: zlboxp  ! previous interval endpoint coordinate: l
  real(fp) :: zxp
  real(fp) :: zyp
  real(fp) :: zlp
  !
  integer :: igridxp  ! x bin of previous endpoint
  integer :: igridyp  ! y bin of previous endpoint
  integer :: igridlp  ! l bin of previous endpoint
  logical insidep  ! .true. if endpoint is inside gridded region
  logical ilinside_pl       !related to 3d halo model only, if .true.
  ! then neutrals inside plasma boundary
  real(fp) :: p1,dtl,deltp,rb_u(3)
  real(fp) :: z_kappa ! decay factor [1/sec]
  real(fp), save :: dxdltrk,dydltrk
  real(fp), save :: pta(3) ! start {x,y,l} point of track
  real(fp) :: ptb(3) ! end {x,y,l} point of track
  real(fp) :: dl_chk
  integer :: n_j ! transp zone index, needed for levmod_halo=2 option
  ! to collect zone avereaged halo source
  integer :: n_iz ! 2d mc zone index, needed for levmod_halo=2 option
  ! to collect zone avereaged halo source
  real(fp), save :: isiztt !box size for beam neutrals only
  integer, save :: inzbx,inzbx_inplasma(2)
  real(fp), save :: wght_trk
  !
  save zlenp,zprobp,zlap,zxap,zyap,zxboxp,zyboxp,zlboxp
  !
  save igridxp,igridyp,igridlp,insidep
  !
  !-----------------------------------------------------------------------
  !  other local variables
  !
  logical inside,iin,idone
  logical, save :: il_first,il_iadr0,il_3dhalo_box
  !
  !-----------------------------------------------------------------------
  !
  n_j=0
  n_iz=0
  jba_0=0
  if(zlen.eq.0.0_fp) then
     isiztt=nxbox*nybox*nlbox*nbbox !box size for beam neutrals only
     il_3dhalo_box=.false.   ! gather 3d halo neutrals into a box
     if(levmod_halo.eq.2.and.nsplt.gt.0) then
        il_3dhalo_box=.true.  ! gather 3d halo neutrals into a box
     end if

     !
     !     initialization
     dxdltrk=0.0_fp !dx/dltrk
     dydltrk=0.0_fp !dy/dltrk
     !
     if(ideptrk.lt.2) then
        idtr = i_trk0
        idgen = 0
        wght_trk=1.0_fp ! track weight = xninj_box/ndepbox_use
        ! for a  ' 0' generation will be normalized in depbox
        xz=xzbeams(lsbeam) !z of neutral
        xa=abeams(lsbeam)  !amu  of neutral
        zxap=xbgrid+xlbapa(nbsbox(ldep))*dxdlbox
        zyap=ybgrid+xlbapa(nbsbox(ldep))*dydlbox
        zlap = xlbapa(nbsbox(ldep))
        dxdltrk=dlbxdtk*dxdlbox !dx/dltrk
        dydltrk=dlbxdtk*dydlbox !dy/dltrk
     elseif(ideptrk.ge.2) then
        idtr = ideptrk - 2
        ldep = itrackh(idtr)%idbox
        idgen = itrackh(idtr)%idgen
        wght_trk=itrackh(idtr)%weight
        xz=itrackh(idtr)%zspecies !z of neutral
        xa=itrackh(idtr)%aspecies   !amu  of neutral
        zxap=itrackh(idtr)%pta(1)*100.0_fp ![m] to c[m]
        zyap=itrackh(idtr)%pta(2)*100.0_fp ![m] to c[m]
        zlap=itrackh(idtr)%pta(3)*100.0_fp ![m] to c[m]
        rb_u(:)=(itrackh(idtr)%ptb(:)-itrackh(idtr)%pta(:))
        rb_u(:)=rb_u(:)/sqrt(sum(rb_u(:)*rb_u(:))) !unit vector in box coordinate
        dlbxdtk=rb_u(3)
        dxdltrk=rb_u(1)
        dydltrk=rb_u(2)
        dl_chk=dlbxdtk**2+dxdltrk**2+dydltrk**2
     end if
     !
     zlenp=zlen
     zprobp=zprob
     !     starting coordinates, [cm]
     zlboxp=zlap
     zxboxp=zxap
     zyboxp=zyap
     !
     !     starting zone
     call boxfind(zxboxp,zyboxp,zlboxp, &
          igridxp,igridyp,igridlp,insidep)

     inzbx=0                !# of box zones was intercepted by track
     inzbx_inplasma=0       !# of box zones when track inters to plasma region
     il_first=.true.        !use this logic to save first track zone in plsama region
     il_iadr0=.true.        !use this logic to save track's input coordinates
     pta(1) =zxap*0.01_fp   ![cm] to [m]
     pta(2) =zyap*0.01_fp   ![cm] to [m]
     pta(3) =zlboxp*0.01_fp ![cm] to [m]
     !
  else
     if(ideptrk.lt.2) then !for fast ptcls only
        jba_0=int((1.0_fp+xksid)*0.5_fp*real(nznbma0))+1
        jba_0=min(nznbma0,max(1,jba_0))
     end if

     !
     !     tracking
     !
     if(zlen.le.zlenp) then ! lengths must be increasing
        write(nonlin,*) &
             ' ?? nubeam addn0: strict increasing length sequence ', &
             'expected.'
        return
     end if

     if(zprob.gt.zprobp) then ! probabilities must decrease
        write(nonlin,*) &
             ' ?? nubeam addn0: non-increasing probability sequence ', &
             'expected.'
        return
     end if
     !
     !     new coordinates
     zlbox=zlap+dlbxdtk*zlen
     !
     !     mg,03/28/2012 this is a bug        zxbox=zxap+zlbox*dxdlbox
     !     mg,03/28/2012 this is a bug        zybox=zyap+zlbox*dydlbox
     ! verified few times mg 05/31/12
     ! must be dltrk**2 = dx**2 +dy**2+dlbx**2
     ! where, dltrk = zlen
     ! as consequence dlbxdtk**2*(1+dxdlbox**2+dydlbox**2)=1
     ! dlbxdtk = dl/dltrack - dl  is length along the grid-normal line,
     !                       dltrack is length along the particle trajectory
     ! dxdlbox = dx/dl !it is not dx/dltrack
     ! dydlbox = dy/dl !it is not dy/dltrack
     !
     ! dxdltrk = dlbxdtk*dxdlbox !dx/dltrk
     zxbox=zxap+zlen*dxdltrk
     ! dydltrk = dlbxdtk*dxdlbox !dy/dltrk
     zybox=zyap+zlen*dydltrk
     !
     !     new endpoint zone
     !
     call boxfind(zxbox,zybox,zlbox, &
          igridx,igridy,igridl,inside)

     !
     !     save coorinates
     !
     zx=zxboxp
     zy=zyboxp
     zl=zlboxp
     ix=igridxp
     iy=igridyp
     il=igridlp
     iin=insidep
     !
     zlenp=zlen
     p1=zprob
     zprobp=zprob
     !
     zxboxp=zxbox
     zyboxp=zybox
     zlboxp=zlbox
     igridxp=igridx
     igridyp=igridy
     igridlp=igridl
     insidep=inside
     !
     !     loop to add to zones in box
     !
     idone=.false.
     ptb(1)=zx
     ptb(2)=zy
     ptb(3)=zl
     z_kappa=zmfpi*vneut
     !
100  continue
     !
     if(.not.idone) then
        !
        zxp=zx
        zyp=zy
        zlp=zl
        call boxinc(zx,zy,zl,ix,iy,il,iin, &
             zxbox,zybox,zlbox,igridx,igridy,igridl,inside, &
             idone,ilinside_pl,n_j,n_iz,iadr,zdl,.true.)
        if(ideptrk.ge.2.and.il_first.and..not.ilinside_pl) then
           call infolog &
                (' addn0: set ilinside_pl=.true.')
           ilinside_pl=.true.
        end if
        if(zdl.eq.0.0_fp) then
           goto 100        ! this track did not evolve, perhaps on the boundary
        end if
        !     find time (dtl) in segment and decay factor (deltp)
        if(.not.ilinside_pl) then
           dtl = p_delt(zdl/vneut,0.0_fp,deltp)
        else
           dtl = p_delt(zdl/vneut,zmfpi*vneut,deltp)
        end if
        !
        if(levmod_halo.eq.2.and. n_j.gt.0.and.(ilinside_pl))then
           if(itrackh(idtr)%type.eq.'THERM') then

              call addn0_3dhalo(n_j,n_iz,lsbeam_parent, &
                   p1*dtl*wght_trk,z_kappa)

           end if
        end if
        if(iadr.gt.0)then
           iadra=nadrbox(ldep)+iadr-1
           if(ideptrk.le.1) then
              boxn0(iadra)=boxn0(iadra)+p1*dtl*wght_trk !sec*#/sec = #, fast neutrals only
           end if
           if((ideptrk.ge.2).and.(il_3dhalo_box).and.(ilinside_pl)) &
                then
              if(idgen.eq.1) then
                 iadrah=iadra+isiztt
                 boxn0(iadrah)=boxn0(iadrah)+p1*dtl*wght_trk !sec*#/sec = #, halo neutrals only
              elseif(idgen.gt.1) then
                 iadrah1=iadra+isiztt*2
                 boxn0(iadrah1)=boxn0(iadrah1)+p1*dtl*wght_trk !sec*#/sec = #, halo neutrals only

              end if
           end if
           if(ideptrk.ge.1) then !only tracks that are in plasma region not in vacuum
              ! and levmod_halo.eq.2
              if(il_iadr0) then !save track's input coordinates
                 pta(1)=zxp*0.01_fp ![cm] to [m]
                 pta(2)=zyp*0.01_fp
                 pta(3)=zlp*0.01_fp
                 il_iadr0 =.false.
              end if
              ptb(1)=zx*0.01_fp ![cm] to [m]
              ptb(2)=zy*0.01_fp
              ptb(3)=zl*0.01_fp
              inzbx=inzbx+1 !# of box zones was intercepted by track
              if(levmod_halo.eq.2)then
                 if(ilinside_pl)then
                    inzbx_inplasma(2)=inzbx
                    if(il_first)then
                       inzbx_inplasma(1)=inzbx
                       if(p1.ne.1.0_fp) then
                          call errorlog &
                               (' addn0: iinplasma(1), p1<1')
                          call errorlog &
                               (' addn0: idtr=',idtr)
                       end if
                       il_first=.false.
                    end if
                 end if
                 !     gather  neutral tracks
                 ! note: '0' generation traks have weight = xninj_box/ndepbox_use, it will be taken into account later in depbox
                 !        track's weight from generation >'0'  already stored in depbox_all/nbi_halo_new_tr
                 call fill_halo(idtr,ldep,zdl*0.01_fp,z_kappa, &
                      inzbx,inzbx_inplasma,idgen,xz,xa, &
                      vneut*0.01_fp,jba_0,p1,pta,ptb) !track's  weight is already set in depbox_all for idgen >1
                 !for idgen = 0 will be set as  xninj_box/ndepbox_use
              end if
           end if            !ideptrk >= 1
           p1=p1*deltp      !note: if zmfpi=0 then deltp=1

        end if               !iadr.gt.0
        !
        goto 100
        !
     end if
     !
     !     exit
     !
  end if                     ! initialization vs. subsequent call
  !
  return
contains
  !
  ! ------------------ pdelt -------------------
  ! returns the average amount of time a particle
  ! exists in a time interval delt given a decay rate of kappa
  !
  real(fp)function p_delt(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.e-4_fp) then
       z = (1.0_fp-0.5_fp*x*(1.-x/3.0_fp))
       p_delt = delt*z
       decay = 1.-x*z
    else if (x>40) then
       p_delt = 1./kappa
       decay = 0.0_fp
    else
       decay = exp(-x)
       p_delt = (1.-decay)/kappa
    end if
  end function p_delt

  subroutine addn0_3dhalo(j_tr,j_nz,lsbeam_p, wght_add, kappa)
    use iso_c_binding, only: fp => c_double

    !     "#/s; thermal neutral source broken down by fast neutral specie
    use physconst_mod, only : evptocm_sec
    use nbi_com, only : lsbeam, nsbeam,sn0s3dh
    use  nbi_com, only : ti,omegag,voln0,evol0,ovol0
    use  nbi_com, only : voln02,evol02,ovol02
    implicit none
    integer, intent(in) :: j_tr ! transp zone index
    integer, intent(in) :: j_nz ! 2d mc zone index
    integer, intent(in) :: lsbeam_p !lsbeam from 'parent'
    real(fp), intent(in)  :: kappa ! decay factor
    integer :: ig,igi            ! thermal particle index, [1/sec]
    real(fp), intent(in) :: wght_add !weight to be added, [#]
    real(fp) :: e_neut ![kev/amu]
    real(fp) :: evptocm_sec2
    evptocm_sec2=evptocm_sec*evptocm_sec 
    if(lsbeam_p.le.nsbeam.and.lsbeam.gt.nsbeam) then
       !get thermal neutral source from fast ptcls.
       ig = lsbeam-nsbeam

       sn0s3dh(j_tr,ig,lsbeam_p)=sn0s3dh(j_tr,ig, &
            lsbeam_p) + wght_add*kappa  !#/sec
    end if
    if(il_3dhalo_den) then  ! eligible to use halo density, vol av halo e and
       ! vol av halo omega from 3d halo model
       ! (not frantic one, set in parms.for)
       ! all beams have to be set in 3d-boxes

       if(lsbeam.gt.nsbeam) then !only thermal neutrals
          e_neut=vneut*vneut/(evptocm_sec2)*1.e-3_fp !this is simplified formula
          !(vneut*0.01_fp)**2*10._fp/(1.384e6_fp**2)
          ! cm**2/s**2 to kev/amu
          igi=lsbeam-nsbeam
          ig=lsbeam_p
          if(lsbeam_p.gt.nsbeam)then
             ig=lsbeam_p-nsbeam
          end if

          nbi_den0mn(igi,j_tr,ig)=nbi_den0mn(igi,j_tr,ig) + wght_add !#
          nbi_en0mn(igi,j_tr,ig)=nbi_en0mn(igi,j_tr,ig) + &
               e_neut*wght_add !#
          voln0(igi,j_tr)=voln0(igi,j_tr) + wght_add
          evol0(igi,j_tr)=evol0(igi,j_tr) + wght_add*ti(j_tr,2)
          ovol0(igi,j_tr)=ovol0(igi,j_tr) + wght_add* &
               0.5_fp*(omegag(j_tr)+omegag(j_tr+1)) !rotation,[1/sec]
          voln02(igi,j_nz)=voln02(igi,j_nz) + wght_add
          evol02(igi,j_nz)=evol02(igi,j_nz) + wght_add*ti(j_tr,2)
          ovol02(igi,j_nz)=ovol02(igi,j_nz) + wght_add* &
               0.5_fp*(omegag(j_tr)+omegag(j_tr+1)) !rotation,[1/sec]
       end if
    end if
  end subroutine addn0_3dhalo
end subroutine addn0
