Generates a neutral beam velocity vector
that passes through cell at ind with weight weight
| Type | Intent | Optional | Attributes | Name | ||
|---|---|---|---|---|---|---|
| integer, | intent(in), | dimension(3) | :: | ind | Cell index |
|
| integer, | intent(in) | :: | neut_type | Neutral Type (1=Full,2=Half,3=Third) |
||
| real(kind=Float64), | intent(out), | dimension(3) | :: | vnbi | Normalized Velocity |
|
| real(kind=Float64), | intent(out) | :: | weight | Weigth/probability of trajectory |
subroutine mc_nbi_cell(ind, neut_type, vnbi, weight)
!+ Generates a neutral beam velocity vector
!+ that passes through cell at `ind` with weight `weight`
integer, dimension(3), intent(in) :: ind
!+ Cell index
integer, intent(in) :: neut_type
!+ Neutral Type (1=Full,2=Half,3=Third)
real(Float64), dimension(3), intent(out) :: vnbi
!+ Normalized Velocity
real(Float64), intent(out) :: weight
!+ Weigth/probability of trajectory
real(Float64), dimension(3) :: rc !! Center of cell in uvw coords
real(Float64), dimension(3) :: uvw_rf !! End position in xyz coords
real(Float64), dimension(3) :: xyz_rf !! End position in xyz coords
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(3) :: randomu !! uniform random numbers
real(Float64) :: sqrt_rho, theta, vy, vz, theta_y, theta_z, py, pz
integer :: i, j
logical :: valid_trajectory
rc = [beam_grid%xc(ind(1)), beam_grid%yc(ind(2)), beam_grid%zc(ind(3))]
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 position in the cell
call randu(randomu)
uvw_rf = rc + (randomu - 0.5)*beam_grid%dr
xyz_rf = matmul(nbi%inv_basis, uvw_rf - nbi%src)
xyz_ray = xyz_rf - xyz_src
xyz_ray = xyz_ray/norm2(xyz_ray)
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)
vnbi = nbi%vinj*uvw_ray/norm2(uvw_ray)/sqrt(real(neut_type))
exit rejection_loop
enddo rejection_loop
!Set Default trajectory in case rejection sampling fails
if(.not.valid_trajectory)then
call randu(randomu)
uvw_rf = rc + (randomu - 0.5)*beam_grid%dr
uvw_ray = uvw_rf - nbi%src
vnbi = nbi%vinj*uvw_ray/norm2(uvw_ray)/sqrt(real(neut_type))
endif
!! Find probability of trajectory
vy = xyz_ray(2)/xyz_ray(1)
vz = xyz_ray(3)/xyz_ray(1)
theta_y = atan(vy + xyz_src(2)/nbi%focy)
theta_z = atan(vz + xyz_src(3)/nbi%focz)
py = (1.0/(1.0 + (vy + xyz_src(2)/nbi%focy)**2)) * &
exp(-(theta_y**2)/(2*nbi%divy(neut_type)**2)) / &
sqrt(2*nbi%divy(neut_type)**2)
pz = (1.0/(1.0 + (vz + xyz_src(3)/nbi%focy)**2)) * &
exp(-(theta_z**2)/(2*nbi%divz(neut_type)**2)) / &
sqrt(2*nbi%divz(neut_type)**2)
weight = py*pz
end subroutine mc_nbi_cell