Calculate neutron emission rate using a fast-ion distribution function F(E,p,r,z)
subroutine neutron_f
!+ Calculate neutron emission rate using a fast-ion distribution function F(E,p,r,z)
integer :: ir, iz, ie, ip, iphi, nphi
type(LocalProfiles) :: plasma
type(LocalEMFields) :: fields
real(Float64) :: eb,pitch,r,z
real(Float64) :: erel, rate
real(Float64), dimension(3) :: rg, ri
real(Float64), dimension(3) :: vi
real(Float64), dimension(3) :: uvw, uvw_vi
real(Float64) :: vnet_square, factor
real(Float64) :: maxcnt, inv_maxcnt, cnt
allocate(neutron%weight(fbm%nenergy,fbm%npitch,fbm%nr,fbm%nz))
neutron%weight = 0.d0
nphi = 20
maxcnt=fbm%nr*fbm%nz
inv_maxcnt = 100.d0/maxcnt
cnt=0.0
!$OMP PARALLEL DO schedule(guided) private(fields,vi,ri,rg,pitch,eb,&
!$OMP& ir,iz,ie,ip,iphi,plasma,factor,uvw,uvw_vi,vnet_square,rate,erel)
z_loop: do iz = 1, fbm%nz
r_loop: do ir=1, fbm%nr
cnt = cnt+1
!! Calculate position
uvw(1) = fbm%r(ir)
uvw(2) = 0.d0
uvw(3) = fbm%z(iz)
call uvw_to_xyz(uvw, rg)
!! Get fields
call get_fields(fields,pos=rg)
if(.not.fields%in_plasma) cycle r_loop
factor = 2*pi*fbm%r(ir)*fbm%dE*fbm%dp*fbm%dr*fbm%dz/nphi
!! Loop over energy/pitch/phi
pitch_loop: do ip = 1, fbm%npitch
pitch = fbm%pitch(ip)
energy_loop: do ie =1, fbm%nenergy
eb = fbm%energy(ie)
gyro_loop: do iphi=1, nphi
call gyro_correction(fields,eb,pitch,ri,vi)
!! Get plasma parameters at particle position
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]
erel = v2_to_E_per_amu*inputs%ab*vnet_square ![kev]
!! Get neutron production rate
call get_neutron_rate(plasma, erel, rate)
neutron%weight(ie,ip,ir,iz) = neutron%weight(ie,ip,ir,iz) &
+ rate*factor
rate = rate*fbm%f(ie,ip,ir,iz)*factor
!! Store neutrons
call store_neutrons(rate)
enddo gyro_loop
enddo energy_loop
enddo pitch_loop
if (inputs%verbose.ge.2)then
WRITE(*,'(f7.2,"% completed",a,$)') cnt*inv_maxcnt,char(13)
endif
enddo r_loop
enddo z_loop
!$OMP END PARALLEL DO
if(inputs%verbose.ge.1) then
write(*,'(T4,A,ES14.5," [neutrons/s]")') 'Rate: ',sum(neutron%rate)
write(*,'(30X,a)') ''
endif
call write_neutrons()
end subroutine neutron_f