Generates a neutral beam particle trajectory
Type | Intent | Optional | Attributes | Name | ||
---|---|---|---|---|---|---|
real(kind=Float64), | intent(out), | dimension(3) | :: | vnbi | Velocity [cm/s] |
|
integer, | intent(in) | :: | efrac | Beam neutral type (1,2,3) |
||
real(kind=Float64), | intent(out), | dimension(3) | :: | rnbi | Starting position on beam_grid |
|
logical, | intent(out) | :: | err | Error Code |
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