mc_nbi Subroutine

subroutine mc_nbi(vnbi,efrac,rnbi,err)
    !+ Generates a neutral beam particle trajectory
    integer, intent(in)                      :: efrac
        !+ Beam neutral type (1,2,3)
    real(Float64), dimension(3), intent(out) :: vnbi
        !+ Velocity [cm/s]
    real(Float64), dimension(3), intent(out) :: rnbi
        !+ Starting position on [[libfida:beam_grid]]
    logical, intent(out)                     :: err
        !+ Error Code

    real(Float64), dimension(3) :: r_exit
    real(Float64), dimension(3) :: uvw_src    !! Start position on ion source
    real(Float64), dimension(3) :: xyz_src    !! Start position on ion source
    real(Float64), dimension(3) :: uvw_ray    !! NBI velocity in uvw coords
    real(Float64), dimension(3) :: xyz_ray    !! NBI velocity in xyz coords
    real(Float64), dimension(3) :: xyz_ape    !! Aperture plane intersection point
    real(Float64), dimension(2) :: randomu    !! uniform random numbers
    real(Float64), dimension(2) :: randomn    !! normal random numbers
    real(Float64) :: length, sqrt_rho, theta
    integer :: i, j
    logical :: inp, valid_trajectory

    err = .False.
    valid_trajectory = .False.
    rejection_loop: do i=1,1000
        call randu(randomu)
        select case (nbi%shape)
            case (1)
                ! Uniformally sample in rectangle
                xyz_src(1) =  0.d0
                xyz_src(2) =  nbi%widy * 2.d0*(randomu(1)-0.5d0)
                xyz_src(3) =  nbi%widz * 2.d0*(randomu(2)-0.5d0)
            case (2)
                ! Uniformally sample in ellipse
                sqrt_rho = sqrt(randomu(1))
                theta = 2*pi*randomu(2)
                xyz_src(1) = 0.d0
                xyz_src(2) = nbi%widy*sqrt_rho*cos(theta)
                xyz_src(3) = nbi%widz*sqrt_rho*sin(theta)
        end select

        !! Create random velocity vector
        call randn(randomn)
        xyz_ray(1)= 1.d0
        xyz_ray(2)=(-xyz_src(2)/nbi%focy + tan(nbi%divy(efrac)*randomn(1)))
        xyz_ray(3)=(-xyz_src(3)/nbi%focz + tan(nbi%divz(efrac)*randomn(2)))

        aperture_loop: do j=1,nbi%naperture
            xyz_ape = xyz_ray*nbi%adist(j) + xyz_src
            select case (nbi%ashape(j))
                case (1)
                    if ((abs(xyz_ape(2) - nbi%aoffy(j)).gt.nbi%awidy(j)).or.&
                        (abs(xyz_ape(3) - nbi%aoffz(j)).gt.nbi%awidz(j))) then
                        cycle rejection_loop
                    endif
                case (2)
                    if ((((xyz_ape(2) - nbi%aoffy(j))*nbi%awidz(j))**2 + &
                         ((xyz_ape(3) - nbi%aoffz(j))*nbi%awidy(j))**2).gt. &
                         (nbi%awidy(j)*nbi%awidz(j))**2) then
                        cycle rejection_loop
                    endif
            end select
        enddo aperture_loop
        valid_trajectory = .True.

        !! Convert to beam centerline coordinates to beam grid coordinates
        uvw_src = matmul(nbi%basis,xyz_src) + nbi%src
        uvw_ray = matmul(nbi%basis,xyz_ray)
        exit rejection_loop
    enddo rejection_loop

    !Set Default trajectory in case rejection sampling fails
    if(.not.valid_trajectory)then
        if(inputs%verbose.ge.0) then
            write(*,'(a)') "MC_NBI: Failed to find trajectory though aperture(s). Using beam centerline."
        endif
        uvw_src = nbi%src
        uvw_ray = nbi%axis
    endif
    vnbi = uvw_ray/norm2(uvw_ray)

    !! Determine start position on beam grid
    call grid_intersect(uvw_src,vnbi,length,rnbi,r_exit)
    if(length.le.0.0)then
        err = .True.
        nbi_outside = nbi_outside + 1
    endif

    !! Check if start position is in the plasma
    call in_plasma(rnbi,inp)
    if(inp)then
        if(inputs%verbose.ge.0) then
            write(*,'(a)') "MC_NBI: A beam neutral has started inside the plasma."
            write(*,'(a)') "Move the beam grid closer to the source to fix"
        endif
        stop
    endif

    !! Determine velocity of neutrals corrected by efrac
    vnbi = vnbi*nbi%vinj/sqrt(real(efrac))
end subroutine mc_nbi