grid_intersect Subroutine

  1. fidasim.f90
  2. libfida
  3. grid_intersect

public subroutine grid_intersect(r0, v0, length, r_enter, r_exit, center_in, lwh_in, passive)

Calculates a particles intersection length with the beam_grid

Arguments

TypeIntentOptionalAttributesName
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

Calls

proc~~grid_intersect~~CallsGraph proc~grid_intersect grid_intersect proc~in_passive_grid in_passive_grid proc~grid_intersect->proc~in_passive_grid proc~uvw_to_cyl uvw_to_cyl proc~in_passive_grid->proc~uvw_to_cyl proc~approx_le approx_le proc~in_passive_grid->proc~approx_le proc~approx_ge approx_ge proc~in_passive_grid->proc~approx_ge proc~approx_eq approx_eq proc~approx_le->proc~approx_eq proc~approx_ge->proc~approx_eq

Called by

proc~~grid_intersect~~CalledByGraph proc~grid_intersect grid_intersect proc~make_diagnostic_grids make_diagnostic_grids proc~make_diagnostic_grids->proc~grid_intersect proc~mc_nbi mc_nbi proc~mc_nbi->proc~grid_intersect proc~fida_weights_los fida_weights_los proc~fida_weights_los->proc~grid_intersect proc~circle_grid_intersect circle_grid_intersect proc~circle_grid_intersect->proc~grid_intersect program~fidasim fidasim program~fidasim->proc~make_diagnostic_grids program~fidasim->proc~fida_weights_los proc~ndmc ndmc program~fidasim->proc~ndmc proc~read_distribution read_distribution program~fidasim->proc~read_distribution proc~read_mc read_mc proc~read_mc->proc~circle_grid_intersect proc~ndmc->proc~mc_nbi proc~read_distribution->proc~read_mc

Contents

Source Code

grid_intersect

Source Code

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