Calculate Passive FIDA emission using a Monte Carlo Fast-ion distribution
subroutine pfida_mc
!+ Calculate Passive FIDA emission using a Monte Carlo Fast-ion distribution
integer :: iion,igamma,ngamma
type(FastIon) :: fast_ion
type(LocalEMFields) :: fields
type(LocalProfiles) :: plasma
real(Float64) :: phi
real(Float64), dimension(3) :: ri !! start position
real(Float64), dimension(3) :: vi !! velocity of fast ions
real(Float64), dimension(3) :: xyz_vi
!! Determination of the CX probability
real(Float64), dimension(nlevs) :: denn !! neutral dens (n=1-4)
real(Float64), dimension(nlevs) :: rates !! CX rates
!! Collisiional radiative model along track
real(Float64), dimension(nlevs) :: states ! Density of n-states
type(ParticleTrack), dimension(pass_grid%ntrack) :: tracks
integer :: ntrack
logical :: los_intersect
integer :: jj !! counter along track
real(Float64) :: photons !! photon flux
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_pfida)/particles%nparticle)
endif
if(inputs%verbose.ge.1) then
write(*,'(T6,"# of markers: ",i10)') int(particles%nparticle*ngamma,Int64)
endif
!$OMP PARALLEL DO schedule(dynamic,1) private(iion,igamma,fast_ion,vi,ri,phi,tracks,s,c,&
!$OMP& randomu,plasma,fields,ntrack,jj,rates,denn,los_intersect,states,photons,xyz_vi)
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
ri(1) = fast_ion%r*c
ri(2) = fast_ion%r*s
ri(3) = fast_ion%z
if(inputs%dist_type.eq.2) then
!! Get electomagnetic fields
call get_fields(fields, pos=ri, input_coords=1, output_coords=1)
!! Correct for gyro motion and get particle position and velocity
call gyro_correction(fields, fast_ion%energy, fast_ion%pitch, ri, vi)
else !! Full Orbit
!! Calculate velocity vector
vi(1) = c*fast_ion%vr - s*fast_ion%vt
vi(2) = s*fast_ion%vr + c*fast_ion%vt
vi(3) = fast_ion%vz
endif
xyz_vi = matmul(beam_grid%inv_basis,vi)
!! Track particle through grid
call track_cylindrical(ri, vi, tracks, ntrack, los_intersect)
if(.not.los_intersect) cycle gamma_loop
if(ntrack.eq.0) cycle gamma_loop
!! Calculate CX probability
call get_plasma(plasma, pos=ri, input_coords=1)
call bt_cx_rates(plasma, plasma%denn, xyz_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
!! Calculate the spectra produced in each cell along the path
loop_along_track: do jj=1,ntrack
call get_plasma(plasma,pos=tracks(jj)%pos,input_coords=1)
call colrad(plasma,beam_ion, xyz_vi, tracks(jj)%time, states, denn, photons)
call store_fida_photons(tracks(jj)%pos, xyz_vi, photons, fast_ion%class,passive=.True.)
enddo loop_along_track
enddo gamma_loop
enddo loop_over_fast_ions
!$OMP END PARALLEL DO
#ifdef _MPI
call parallel_sum(spec%pfida)
#endif
end subroutine pfida_mc