Calculates a particles intersection length with the beam_grid
| Type | Intent | Optional | Attributes | Name | ||
|---|---|---|---|---|---|---|
| real(kind=Float64), | intent(in), | dimension(3) | :: | r0 | Initial position of particle [cm] | |
| real(kind=Float64), | intent(in), | dimension(3) | :: | v0 | Velocity of particle [cm/s] | |
| real(kind=Float64), | intent(out) | :: | length | Intersection length [cm] | ||
| real(kind=Float64), | intent(out), | dimension(3) | :: | r_enter | Point where particle enters | |
| real(kind=Float64), | intent(out), | dimension(3) | :: | r_exit | Point where particle exits | |
| real(kind=Float64), | intent(in), | optional | dimension(3) | :: | center_in | Alternative grid center | 
| real(kind=Float64), | intent(in), | optional | dimension(3) | :: | lwh_in | Alternative grid [length,width,height] | 
| logical, | intent(in), | optional | :: | passive | Calculates a particles intersection length with the pass_grid | 
subroutine grid_intersect(r0, v0, length, r_enter, r_exit, center_in, lwh_in, passive)
    !+ Calculates a particles intersection length with the [[libfida:beam_grid]]
    real(Float64), dimension(3), intent(in)           :: r0
        !+ Initial position of particle [cm]
    real(Float64), dimension(3), intent(in)           :: v0
        !+ Velocity of particle [cm/s]
    real(Float64), intent(out)                        :: length
        !+ Intersection length [cm]
    real(Float64), dimension(3), intent(out)          :: r_enter
        !+ Point where particle enters
    real(Float64), dimension(3), intent(out)          :: r_exit
        !+ Point where particle exits
    real(Float64), dimension(3), intent(in), optional :: center_in
        !+ Alternative grid center
    real(Float64), dimension(3), intent(in), optional :: lwh_in
        !+ Alternative grid [length,width,height]
    logical, intent(in), optional                     :: passive
        !+ Calculates a particles intersection length with the [[libfida:pass_grid]]
    real(Float64), dimension(3,6) :: ipnts
    real(Float64), dimension(3) :: ri, vi
    real(Float64), dimension(3) :: center
    real(Float64), dimension(3) :: lwh
    integer, dimension(6) :: side_inter
    integer, dimension(2) :: ind
    integer :: i, j, nunique, ind1, ind2
    real(Float64) :: dlength, max_length
    logical :: ing, pas
    pas = .False.
    if(present(passive)) pas = passive
    if(pas) then
        ing = in_passive_grid(r0)
        if (ing) then
            length = 0.d0
            return
        endif
        ri = r0 ; vi = v0
        dlength = 2.0 !cm
        max_length=0.0
        do while (.not.ing)
            ri = ri + vi*dlength ! move dlength
            ing = in_passive_grid(ri)
            max_length = max_length + dlength
            if(max_length.gt.1d3) then
                length = 0.d0
                return
            endif
        enddo
        r_enter = ri
        do while (ing)
            ri = ri + vi*dlength
            ing =  in_passive_grid(ri)
        enddo
        r_exit = ri
        length = sqrt(sum((r_exit-r_enter)**2))
    else
        if (present(center_in)) then
            center = center_in
        else
            center = beam_grid%center
        endif
        if (present(lwh_in)) then
            lwh = lwh_in
        else
            lwh = beam_grid%lwh
        endif
        side_inter = 0
        ipnts = 0.d0
        do i=1,6
            j = int(ceiling(i/2.0))
            if (j.eq.1) ind = [2,3]
            if (j.eq.2) ind = [1,3]
            if (j.eq.3) ind = [1,2]
            if (abs(v0(j)).gt.0.d0) then
                ipnts(:,i) = r0 + v0*( ( (center(j) + &
                             (mod(i,2) - 0.5)*lwh(j)) - r0(j))/v0(j) )
                if ((abs(ipnts(ind(1),i) - center(ind(1))).le.(0.5*lwh(ind(1)))).and. &
                    (abs(ipnts(ind(2),i) - center(ind(2))).le.(0.5*lwh(ind(2))))) then
                    side_inter(i) = 1
                endif
            endif
        enddo
        length = 0.d0
        r_enter = r0
        r_exit  = r0
        ind1=0
        ind2=0
        if (sum(side_inter).ge.2) then
            ! Find first intersection side
            i=1
            do while (i.le.6)
                if(side_inter(i).eq.1) exit
                i=i+1
            enddo
            ind1=i
            !Find number of unique points
            nunique = 0
            do i=ind1+1,6
                if (side_inter(i).ne.1) cycle
                if (sqrt( sum( ( ipnts(:,i)-ipnts(:,ind1) )**2 ) ).gt.0.001) then
                    ind2=i
                    nunique = 2
                    exit
                endif
            enddo
            if(nunique.eq.2) then
                vi = ipnts(:,ind2) - ipnts(:,ind1)
                if (dot_product(v0,vi).gt.0.0) then
                    r_enter = ipnts(:,ind1)
                    r_exit  = ipnts(:,ind2)
                else
                    r_enter = ipnts(:,ind2)
                    r_exit  = ipnts(:,ind1)
                endif
                length = sqrt(sum((r_exit - r_enter)**2))
            endif
        endif
    endif
end subroutine grid_intersect