Calculate Passive FIDA emission using a Fast-ion distribution function F(E,p,r,z)
subroutine pfida_f
!+ Calculate Passive FIDA emission using a Fast-ion distribution function F(E,p,r,z)
integer :: i,j,k,ic,ncell,is
integer(Int64) :: iion
real(Float64), dimension(3) :: ri !! start position
real(Float64), dimension(3) :: vi !! velocity of fast ions
real(Float64), dimension(3) :: xyz_vi
real(Float64) :: denf !! fast-ion density
integer, dimension(3) :: ind !! new actual cell
logical :: los_intersect
!! Determination of the CX probability
type(LocalEMFields) :: fields
type(LocalProfiles) :: plasma
real(Float64), dimension(nlevs) :: rates, rates_is !! CX rates
!! Collisiional radiative model along track
integer :: ntrack
integer :: jj !! counter along track
type(ParticleTrack),dimension(pass_grid%ntrack) :: tracks
real(Float64) :: photons !! photon flux
real(Float64), dimension(nlevs) :: states !! Density of n-states
real(Float64), dimension(nlevs) :: denn
!! Number of particles to launch
real(Float64) :: max_papprox, eb, ptch
integer, dimension(pass_grid%ngrid) :: cell_ind
real(Float64), dimension(pass_grid%nr,pass_grid%nz,pass_grid%nphi) :: papprox
integer(Int32), dimension(pass_grid%nr,pass_grid%nz,pass_grid%nphi) :: nlaunch
!! Estimate how many particles to launch in each cell
papprox=0.d0
do ic=1,pass_grid%ngrid
call ind2sub(pass_grid%dims,ic,ind)
i = ind(1) ; j = ind(2) ; k = ind(3)
call get_plasma(plasma,ind=ind,input_coords=2)
if(.not.plasma%in_plasma) cycle
papprox(i,j,k) = sum(plasma%denn)*plasma%denf
enddo
max_papprox = maxval(papprox)
where (papprox.lt.(max_papprox*1.d-6))
papprox = 0.0
endwhere
ncell = 0
do ic=1,pass_grid%ngrid
call ind2sub(pass_grid%dims,ic,ind)
i = ind(1) ; j = ind(2) ; k = ind(3)
if(papprox(i,j,k).gt.0.0) then
ncell = ncell + 1
cell_ind(ncell) = ic
endif
enddo
call get_nlaunch_pass_grid(inputs%n_pfida, papprox, nlaunch)
if(inputs%verbose.ge.1) then
write(*,'(T6,"# of markers: ",i10)') sum(nlaunch)
endif
!! Loop over all cells that have neutrals
!$OMP PARALLEL DO schedule(dynamic,1) private(ic,i,j,k,ind,iion,vi,xyz_vi,ri,fields, &
!$OMP tracks,ntrack,jj,plasma,rates,rates_is,is,denn,states,photons,denf,eb,ptch,los_intersect)
loop_over_cells: do ic = istart, ncell, istep
call ind2sub(pass_grid%dims,cell_ind(ic),ind)
i = ind(1) ; j = ind(2) ; k = ind(3)
loop_over_fast_ions: do iion=1, nlaunch(i, j, k)
!! Sample fast ion distribution for velocity and position
call mc_fastion_pass_grid(ind, fields, eb, ptch, denf, output_coords=1)
if(denf.le.0.0) cycle loop_over_fast_ions
!! Correct for gyro motion and get particle position and velocity
call gyro_correction(fields, eb, ptch, fbm%A, ri, vi)
xyz_vi = matmul(beam_grid%inv_basis,vi)
!! Find the particles path through the interpolation grid
call track_cylindrical(ri, vi, tracks, ntrack,los_intersect)
if(.not.los_intersect) cycle loop_over_fast_ions
if(ntrack.eq.0) cycle loop_over_fast_ions
!! Calculate CX probability with edge neutrals
call get_plasma(plasma, pos=ri, input_coords=1)
rates = 0.0
do is=1,n_thermal
call bt_cx_rates(plasma, plasma%denn(:,is), thermal_mass(is), xyz_vi, rates_is)
rates = rates + rates_is
enddo
if(sum(rates).le.0.) cycle loop_over_fast_ions
!! Weight CX rates by ion source density
states=rates*denf
!! 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, fbm%A, xyz_vi, tracks(jj)%time, states, denn, photons)
call store_fida_photons(tracks(jj)%pos, xyz_vi, beam_lambda0, photons/nlaunch(i,j,k), passive=.True.)
if(inputs%calc_res.ge.1) then
call store_photon_birth(tracks(1)%pos, photons/nlaunch(i,j,k), spatres%pfida, passive=.True.)
endif
enddo loop_along_track
enddo loop_over_fast_ions
enddo loop_over_cells
!$OMP END PARALLEL DO
#ifdef _MPI
call parallel_sum(spec%pfida)
if(inputs%calc_res.ge.1) then
do jj=1,spec_chords%nchan
call parallel_merge_reservoirs(spatres%pfida(jj))
enddo
endif
#endif
end subroutine pfida_f