track Subroutine

subroutine track(rin, vin, tracks, ncell, los_intersect)
    !+ Computes the path of a neutral through the [[libfida:beam_grid]]
    real(Float64), dimension(3), intent(in)          :: rin
        !+ Initial position of particle
    real(Float64), dimension(3), intent(in)          :: vin
        !+ Initial velocity of particle
    type(ParticleTrack), dimension(:), intent(inout) :: tracks
        !+ Array of [[ParticleTrack]] type
    integer(Int32), intent(out)                      :: ncell
        !+ Number of cells that a particle crosses
    logical, intent(out), optional                   :: los_intersect
        !+ Indicator whether particle intersects a LOS in [[libfida:spec_chords]]

    integer :: cc, i, ii, mind
    integer, dimension(3) :: ind
    logical :: in_plasma1, in_plasma2, in_plasma_tmp, los_inter
    real(Float64) :: dT, dt1, inv_50
    real(Float64), dimension(3) :: dt_arr, dr
    real(Float64), dimension(3) :: vn, inv_vn
    real(Float64), dimension(3) :: ri, ri_tmp, ri_cell
    integer, dimension(3) :: sgn
    integer, dimension(3) :: gdims
    integer, dimension(1) :: minpos

    vn = vin ;  ri = rin ; sgn = 0 ; ncell = 0

    if(dot_product(vin,vin).eq.0.0) then
        return
    endif

    gdims(1) = beam_grid%nx
    gdims(2) = beam_grid%ny
    gdims(3) = beam_grid%nz

    !! define actual cell
    call get_indices(ri,ind)
    ri_cell = [beam_grid%xc(ind(1)), &
               beam_grid%yc(ind(2)), &
               beam_grid%zc(ind(3))]

    do i=1,3
        if (vn(i).gt.0.0) sgn(i) = 1
        if (vn(i).lt.0.0) sgn(i) =-1
        if (vn(i).eq.0.0) vn(i)  = 1.0d-3
    enddo

    dr = beam_grid%dr*sgn
    inv_vn = 1/vn
    inv_50 = 1.0/50.0
    cc=1
    los_inter = .False.
    tracks%time = 0.d0
    tracks%flux = 0.d0
    call in_plasma(ri,in_plasma1)
    track_loop: do i=1,beam_grid%ntrack
        if(cc.gt.beam_grid%ntrack) exit track_loop

        if((spec_chords%inter(ind(1),ind(2),ind(3))%nchan.ne.0) &
            .and.(.not.los_inter))then
            los_inter = .True.
        endif
        dt_arr = abs(( (ri_cell + 0.5*dr) - ri)*inv_vn)
        minpos = minloc(dt_arr)
        mind = minpos(1)
        dT = dt_arr(mind)
        ri_tmp = ri + dT*vn
        call in_plasma(ri_tmp,in_plasma2)
        if(in_plasma1.neqv.in_plasma2) then
            dt1 = 0.0
            track_fine: do ii=1,50
                dt1 = dt1 + dT*inv_50
                ri_tmp = ri + vn*dt1
                call in_plasma(ri_tmp,in_plasma_tmp)
                if(in_plasma2.eqv.in_plasma_tmp) exit track_fine
            enddo track_fine
            tracks(cc)%in_plasma = in_plasma2
            tracks(cc+1)%in_plasma = .not.in_plasma2
            tracks(cc+1)%pos = ri + 0.5*(dt1 + dT)*vn
            tracks(cc)%pos = ri + 0.5*dt1*vn
            tracks(cc+1)%pos = ri + 0.5*(dt1 + dT)*vn
            tracks(cc)%time = dt1
            tracks(cc+1)%time = dT - dt1
            tracks(cc)%ind = ind
            tracks(cc+1)%ind = ind
            cc = cc + 2
        else
            tracks(cc)%in_plasma = in_plasma2
            tracks(cc)%pos = ri + 0.5*dT*vn
            tracks(cc)%time = dT
            tracks(cc)%ind = ind
            cc = cc + 1
        endif
        in_plasma1 = in_plasma2

        ri = ri + dT*vn
        ind(mind) = ind(mind) + sgn(mind)
        ri_cell(mind) = ri_cell(mind) + dr(mind)

        if (ind(mind).gt.gdims(mind)) exit track_loop
        if (ind(mind).lt.1) exit track_loop
    enddo track_loop
    ncell = cc-1
    if(present(los_intersect)) then
        los_intersect = los_inter
    endif

end subroutine track