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)

Calculates a particles intersection length with the beam_grid

Arguments

Type IntentOptional AttributesName
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 beam_grid
real(kind=Float64), intent(out), dimension(3):: r_exit

Point where particle exits beam_grid
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]

Called by

proc~~grid_intersect~~CalledByGraph proc~grid_intersect 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 proc~read_chords read_chords proc~read_chords->proc~grid_intersect proc~read_npa read_npa proc~read_npa->proc~grid_intersect proc~mc_nbi mc_nbi proc~mc_nbi->proc~grid_intersect proc~read_mc read_mc proc~read_mc->proc~circle_grid_intersect proc~ndmc ndmc proc~ndmc->proc~mc_nbi program~fidasim fidasim program~fidasim->proc~fida_weights_los program~fidasim->proc~read_chords program~fidasim->proc~read_npa program~fidasim->proc~ndmc proc~read_distribution read_distribution program~fidasim->proc~read_distribution 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)
    !+ 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 [[libfida:beam_grid]]
    real(Float64), dimension(3), intent(out)          :: r_exit
        !+ Point where particle exits [[libfida:beam_grid]]
    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]

    real(Float64), dimension(3,6) :: ipnts
    real(Float64), dimension(3) :: 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

    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

end subroutine grid_intersect