Calculate Passive NPA flux using a Monte Carlo fast-ion distribution
subroutine pnpa_mc
!+ Calculate Passive NPA flux using a Monte Carlo fast-ion distribution
integer :: iion,igamma,ngamma,npart
type(FastIon) :: fast_ion
real(Float64) :: phi,theta,dtheta
real(Float64), dimension(3) :: ri, rf, rg, vi
integer :: det,j,ichan,ir,nrange,it
type(LocalEMFields) :: fields
type(LocalProfiles) :: plasma
type(GyroSurface) :: gs
real(Float64), dimension(nlevs) :: rates
real(Float64), dimension(nlevs) :: states
real(Float64) :: flux
integer, dimension(3) :: ind
real(Float64), dimension(3) :: uvw, uvw_vi
real(Float64), dimension(2,4) :: gyrange
real(Float64) :: s,c
real(Float64), dimension(1) :: randomu
ngamma = 1
if(particles%axisym.or.(inputs%dist_type.eq.2)) then
ngamma = ceiling(dble(inputs%n_pnpa)/particles%nparticle)
endif
if(inputs%verbose.ge.1) then
write(*,'(T6,"# of markers: ",i10)') int(particles%nparticle*ngamma,Int64)
endif
!$OMP PARALLEL DO schedule(guided) private(iion,igamma,ind,fast_ion,vi,ri,rf,phi,s,c,ir,it,plasma, &
!$OMP& randomu,rg,fields,uvw,uvw_vi,rates,states,flux,det,ichan,gs,nrange,gyrange,theta,dtheta)
loop_over_fast_ions: do iion=istart,particles%nparticle,istep
fast_ion = particles%fast_ion(iion)
if(fast_ion%vabs.eq.0) cycle loop_over_fast_ions
gamma_loop: do igamma=1,ngamma
if(particles%axisym) then
!! Pick random toroidal angle
call randu(randomu)
phi = pass_grid%phi(1)+pass_grid%nphi*pass_grid%dphi*randomu(1)
else
phi = fast_ion%phi
endif
s = sin(phi)
c = cos(phi)
!! Calculate position in machine coordinates
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)
!! Correct for gyro motion and get position and velocity
call gyro_surface(fields, fast_ion%energy, fast_ion%pitch, gs)
detector_loop: do ichan=1,npa_chords%nchan
call npa_gyro_range(ichan, gs, gyrange, nrange)
if(nrange.eq.0) cycle detector_loop
gyro_range_loop: do ir=1,nrange
dtheta = gyrange(2,ir)
theta = gyrange(1,ir) + 0.5*dtheta
call gyro_trajectory(gs, theta, ri, vi)
!! Check if particle hits a NPA detector
call hit_npa_detector(ri, vi ,det, rf, det=ichan)
if(det.ne.ichan) then
if (inputs%verbose.ge.0)then
write(*,*) "PNPA_MC: Missed Detector ",ichan
endif
cycle gyro_range_loop
endif
!! Calculate CX probability with beam and halo neutrals
call get_plasma(plasma, pos=ri)
call bt_cx_rates(plasma, plasma%denn, vi, beam_ion, rates)
if(sum(rates).le.0.) cycle gyro_range_loop
!! Weight CX rates by ion source density
if(particles%axisym) then
states=rates*fast_ion%weight*(pass_grid%nphi*pass_grid%dphi/(2*pi)) &
/(fast_ion%r*pass_grid%dv)/ngamma
else
states=rates*fast_ion%weight/(fast_ion%r*pass_grid%dv)/ngamma
endif
!! Attenuate states as the particle move through plasma
call attenuate(ri,rf,vi,states)
!! Store NPA Flux
flux = (dtheta/(2*pi))*sum(states)*(fast_ion%r*pass_grid%dv)
spread_loop: do it=1,25
theta = gyrange(1,ir) + (it-0.5)*dtheta/25
call gyro_trajectory(gs, theta, ri, vi)
call hit_npa_detector(ri, vi ,det, rf, det=ichan)
if(det.ne.ichan) then
if (inputs%verbose.ge.0)then
write(*,*) "PNPA_MC: Missed Detector ",ichan
endif
cycle spread_loop
endif
call store_npa(det,ri,rf,vi,flux/25,fast_ion%class, passive=.True.)
enddo spread_loop
enddo gyro_range_loop
enddo detector_loop
else !! Full Orbit
!! Convert to beam grid coordinates
call uvw_to_xyz(uvw, ri)
!! Calculate velocity vector
uvw_vi(1) = c*fast_ion%vr - s*fast_ion%vt
uvw_vi(2) = s*fast_ion%vr + c*fast_ion%vt
uvw_vi(3) = fast_ion%vz
vi = matmul(beam_grid%inv_basis,uvw_vi)
!! Check if particle hits a NPA detector
call hit_npa_detector(ri, vi ,det, rf)
if(det.eq.0) cycle gamma_loop
!! Calculate CX probability with beam and halo neutrals
call get_plasma(plasma, pos=ri)
call bt_cx_rates(plasma, plasma%denn, vi, beam_ion, rates)
if(sum(rates).le.0.) cycle gamma_loop
!! Weight CX rates by ion source density
if(particles%axisym) then
states=rates*fast_ion%weight*(pass_grid%nphi*pass_grid%dphi/(2*pi)) &
/(fast_ion%r*pass_grid%dv)/ngamma
else
states=rates*fast_ion%weight/(fast_ion%r*pass_grid%dv)/ngamma
endif
!! Attenuate states as the particle moves though plasma
call attenuate(ri,rf,vi,states)
!! Store NPA Flux
flux = sum(states)*(fast_ion%r*pass_grid%dv)
call store_npa(det,ri,rf,vi,flux,fast_ion%class,passive=.True.)
endif
enddo gamma_loop
enddo loop_over_fast_ions
!$OMP END PARALLEL DO
npart = pnpa%npart
#ifdef _MPI
call parallel_sum(npart)
call parallel_sum(pnpa%flux)
#endif
if(inputs%verbose.ge.1) then
write(*,'(T4,"Number of Passive NPA particles that hit a detector: ",i8)') npart
endif
end subroutine pnpa_mc