Calculates FIDA weights
subroutine fida_weights_mc
!+ Calculates FIDA weights
integer :: i,j,k,ic,ncell
integer(Int64) :: iion,ip
real(Float64), dimension(3) :: ri,rg !! start position
real(Float64), dimension(3) :: vi !! velocity of fast ions
integer, dimension(3) :: ind !! new actual cell
integer,dimension(5) :: neut_types=[1,2,3,4,5]
logical :: los_intersect
!! Determination of the CX rates
type(LocalProfiles) :: plasma
type(LocalEMFields) :: fields
real(Float64), dimension(nlevs) :: rates !! CX rates
!! Collisiional radiative model along track
integer :: ntrack
integer :: jj !! counter along track
type(ParticleTrack),dimension(beam_grid%ntrack) :: tracks
real(Float64) :: photons !! photon flux
real(Float64), dimension(nlevs) :: states !! Density of n-states
real(Float64), dimension(nlevs) :: denn
integer :: nwav
real(Float64) :: etov2, energy, pitch
real(Float64) :: dE, dP, dEdP
real(Float64), dimension(:), allocatable :: ebarr, ptcharr
integer, dimension(1) :: ienergy, ipitch
real(Float64), dimension(3) :: randomu3
!! Number of particles to launch
real(Float64) :: fbm_denf,phase_area
integer, dimension(beam_grid%ngrid) :: cell_ind
real(Float64), dimension(beam_grid%nx,beam_grid%ny,beam_grid%nz) :: papprox
integer(Int32), dimension(beam_grid%nx,beam_grid%ny,beam_grid%nz) :: nlaunch
nwav = inputs%nlambda_wght
!! define arrays
!! define energy - array
allocate(ebarr(inputs%ne_wght))
do i=1,inputs%ne_wght
ebarr(i)=real(i-0.5)*inputs%emax_wght/real(inputs%ne_wght)
enddo
dE = abs(ebarr(2)-ebarr(1))
!! define pitch - array
allocate(ptcharr(inputs%np_wght))
do i=1,inputs%np_wght
ptcharr(i)=real(i-0.5)*2./real(inputs%np_wght)-1.
enddo
dP = abs(ptcharr(2)-ptcharr(1))
dEdP = dE*dP
phase_area = dEdP*real(inputs%np_wght)*real(inputs%ne_wght)
!! allocate storage arrays
allocate(fweight%weight(nwav,inputs%ne_wght,inputs%np_wght,spec_chords%nchan))
allocate(fweight%mean_f(inputs%ne_wght,inputs%np_wght,spec_chords%nchan))
if(inputs%verbose.ge.1) then
write(*,'(T3,"Number of Channels: ",i5)') spec_chords%nchan
write(*,'(T3,"Nlambda: ",i4)') nwav
write(*,'(T3,"Nenergy: ",i3)') inputs%ne_wght
write(*,'(T3,"Maximum Energy: ",f7.2)') inputs%emax_wght
write(*,'(T3,"LOS averaged: ",a)') "False"
endif
!! zero out arrays
fweight%weight = 0.d0
fweight%mean_f = 0.d0
etov2 = 1.d0/(v2_to_E_per_amu*inputs%ab)
!! Estimate how many particles to launch in each cell
papprox=0.d0
do ic=1,beam_grid%ngrid
call ind2sub(beam_grid%dims,ic,ind)
i = ind(1) ; j = ind(2) ; k = ind(3)
call get_plasma(plasma,ind=ind)
if(.not.plasma%in_plasma) cycle
papprox(i,j,k)=(sum(neut%full(:,i,j,k)) + &
sum(neut%half(:,i,j,k)) + &
sum(neut%third(:,i,j,k)) + &
sum(neut%dcx(:,i,j,k)) + &
sum(neut%halo(:,i,j,k)))
enddo
ncell = 0
do ic=1,beam_grid%ngrid
call ind2sub(beam_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(10*inputs%n_fida,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(guided) private(ic,i,j,k,ind,iion,vi,ri,rg,ienergy,ipitch, &
!$OMP tracks,ntrack,jj,plasma,fields,rates,denn,states,photons,energy,pitch, &
!$OMP los_intersect,randomu3,fbm_denf)
loop_over_cells: do ic = istart, ncell, istep
call ind2sub(beam_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 uniformally
call randind(inputs%ne_wght, ienergy)
call randind(inputs%np_wght, ipitch)
call randu(randomu3)
energy = ebarr(ienergy(1)) + dE*(randomu3(1)-0.5)
pitch = ptcharr(ipitch(1)) + dP*(randomu3(2)-0.5)
if(energy.le.0) cycle loop_over_fast_ions
call randu(randomu3)
rg = [beam_grid%xc(i),beam_grid%yc(j),beam_grid%zc(k)] + beam_grid%dr*(randomu3-0.5)
!! Get velocity
call get_fields(fields,pos=rg)
if(.not.fields%in_plasma) cycle loop_over_fast_ions
call gyro_correction(fields,energy,pitch,ri,vi)
fbm_denf = 0.0
if (inputs%dist_type.eq.1) then
call get_ep_denf(energy,pitch,fbm_denf,coeffs=fields%b)
endif
!! Find the particles path through the beam grid
call track(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 beam and halo neutrals
call get_beam_cx_rate(tracks(1)%ind, ri, vi, beam_ion, neut_types, rates)
if(sum(rates).le.0.) cycle loop_over_fast_ions
states=rates*1.d20
!! 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)
call colrad(plasma,beam_ion, vi, tracks(jj)%time, states, denn, photons)
call store_fw_photons(ienergy(1), ipitch(1), &
tracks(jj)%pos, vi, fbm_denf, photons/nlaunch(i,j,k))
enddo loop_along_track
enddo loop_over_fast_ions
enddo loop_over_cells
!$OMP END PARALLEL DO
fweight%weight = ((1.d-20)*phase_area/dEdP)*fweight%weight
fweight%mean_f = ((1.d-20)*phase_area/dEdP)*fweight%mean_f
if(inputs%verbose.ge.1) then
write(*,*) 'write fida weights: ' , time(time_start)
endif
#ifdef _MPI
call parallel_sum(fweight%weight)
call parallel_sum(fweight%mean_f)
if(my_rank().eq.0) call write_fida_weights()
#else
call write_fida_weights()
#endif
end subroutine fida_weights_mc