Calculate neutron flux using a Monte Carlo Fast-ion distribution
subroutine neutron_mc
!+ Calculate neutron flux using a Monte Carlo Fast-ion distribution
integer :: iion, ngamma, igamma
type(FastIon) :: fast_ion
type(LocalProfiles) :: plasma
type(LocalEMFields) :: fields
real(Float64) :: eb, rate
real(Float64), dimension(3) :: ri
real(Float64), dimension(3) :: vi
real(Float64), dimension(3) :: uvw, uvw_vi
real(Float64) :: vnet_square
real(Float64) :: phi, s, c, factor, delta_phi
if(inputs%verbose.ge.1) then
write(*,'(T6,"# of markers: ",i10)') particles%nparticle
endif
!! Correct neutron rate when equilibrium is 3D and MC distribution is 4D
if(particles%axisym.and.(inter_grid%nphi.gt.1)) then
delta_phi = inter_grid%phi(inter_grid%nphi)-inter_grid%phi(1)
delta_phi = delta_phi + delta_phi/(inter_grid%nphi-1)/2 !Add half a cell
factor = delta_phi/(2*pi) * 2 !Riemann sum below assumes coord's are at midpoint of cell
else
factor = 1
endif
rate=0.0
ngamma = 20
!$OMP PARALLEL DO schedule(guided) private(iion,fast_ion,vi,ri,s,c, &
!$OMP& plasma,fields,uvw,uvw_vi,vnet_square,rate,eb,igamma,phi)
loop_over_fast_ions: do iion=istart,particles%nparticle,istep
fast_ion = particles%fast_ion(iion)
if(fast_ion%vabs.eq.0.d0) cycle loop_over_fast_ions
!! Calculate position in machine coordinates
if(particles%axisym) then
s = 0.d0
c = 1.d0
else
phi = fast_ion%phi
s = sin(phi)
c = cos(phi)
endif
uvw(1) = fast_ion%r*c
uvw(2) = fast_ion%r*s
uvw(3) = fast_ion%z
if(inputs%dist_type.eq.2) then
!! Get electomagnetic fields
call get_fields(fields, pos=uvw, input_coords=1)
if(.not.fields%in_plasma) cycle loop_over_fast_ions
gyro_loop: do igamma=1,ngamma
!! Correct for Gyro-motion
call gyro_correction(fields, fast_ion%energy, fast_ion%pitch, ri, vi)
!! Get plasma parameters
call get_plasma(plasma,pos=ri)
if(.not.plasma%in_plasma) cycle gyro_loop
!! Calculate effective beam energy
vnet_square=dot_product(vi-plasma%vrot,vi-plasma%vrot) ![cm/s]
eb = v2_to_E_per_amu*inputs%ab*vnet_square ![kev]
!! Get neutron production rate
call get_neutron_rate(plasma, eb, rate)
rate = rate*fast_ion%weight/ngamma*factor
!! Store neutrons
call store_neutrons(rate, fast_ion%class)
enddo gyro_loop
else
!! Get plasma parameters
call get_plasma(plasma,pos=uvw,input_coords=1)
if(.not.plasma%in_plasma) cycle loop_over_fast_ions
!! Calculate effective beam energy
uvw_vi(1) = fast_ion%vr
uvw_vi(2) = fast_ion%vt
uvw_vi(3) = fast_ion%vz
vi = matmul(beam_grid%inv_basis,uvw_vi)
vnet_square=dot_product(vi-plasma%vrot,vi-plasma%vrot) ![cm/s]
eb = v2_to_E_per_amu*inputs%ab*vnet_square ![kev]
!! Get neutron production rate
call get_neutron_rate(plasma, eb, rate)
rate = rate*fast_ion%weight*factor
!! Store neutrons
call store_neutrons(rate, fast_ion%class)
endif
enddo loop_over_fast_ions
!$OMP END PARALLEL DO
#ifdef _MPI
call parallel_sum(neutron%rate)
#endif
if(inputs%verbose.ge.1) then
write(*,'(T4,A,ES14.5," [neutrons/s]")') 'Rate: ',sum(neutron%rate)
write(*,'(30X,a)') ''
write(*,*) 'write neutrons: ' , time(time_start)
endif
#ifdef _MPI
if(my_rank().eq.0) call write_neutrons()
#else
call write_neutrons()
#endif
end subroutine neutron_mc