Calculates LOS averaged FIDA weights
subroutine fida_weights_los
!+ Calculates LOS averaged FIDA weights
type(LocalProfiles) :: plasma, plasma_cell
type(LocalEMFields) :: fields, fields_cell
real(Float64) :: denf
real(Float64) :: wght, wght_tot
real(Float64) :: photons !! photon flux
real(Float64) :: length
type(ParticleTrack), dimension(beam_grid%ntrack) :: tracks
integer :: nwav
integer(Int32) :: i, j, k, ienergy, cid, cind
integer(Int32) :: ipitch, igyro, icell, ichan
real(Float64), dimension(:), allocatable :: ebarr,ptcharr,phiarr
real(Float64), dimension(:,:), allocatable :: mean_f
real(Float64), dimension(3) :: vi, vi_norm, vp
real(Float64), dimension(3) :: vnbi_f, vnbi_h, vnbi_t, vhalo
real(Float64), dimension(3) :: r_enter, r_exit
real(Float64) :: vabs, dE, dP
!! Determination of the CX probability
real(Float64), dimension(nlevs) :: fdens,hdens,tdens,dcxdens,halodens
real(Float64), dimension(nlevs) :: rates
real(Float64), dimension(nlevs) :: states ! Density of n-states
real(Float64), dimension(nlevs) :: denn ! Density of n-states
!! COLRAD
real(Float64) :: dt, max_dens, dlength, sigma_pi
type(LOSInters) :: inter
real(Float64) :: eb, ptch, phi
!! Solution of differential equation
integer, dimension(3) :: ind !!actual cell
real(Float64), dimension(3) :: ri
integer(Int32) :: ntrack
logical :: inp
real(Float64):: etov2, dEdP
nwav = inputs%nlambda_wght
!! 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
!! define gyro - array
allocate(phiarr(inputs%nphi_wght))
do i=1,inputs%nphi_wght
phiarr(i)=real(i-0.5)*2.d0*pi/real(inputs%nphi_wght)
enddo
!! allocate storage arrays
allocate(fweight%mean_f(inputs%ne_wght,inputs%np_wght,spec_chords%nchan))
allocate(fweight%weight(nwav,inputs%ne_wght,inputs%np_wght,spec_chords%nchan))
allocate(mean_f(inputs%ne_wght,inputs%np_wght))
!! zero out arrays
fweight%weight = 0.d0
fweight%mean_f = 0.d0
mean_f = 0.d0
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,"Npitch: ",i3)') inputs%np_wght
write(*,'(T3,"Ngyro: ", i3)') inputs%nphi_wght
write(*,'(T3,"Maximum Energy: ",f7.2)') inputs%emax_wght
write(*,'(T3,"LOS averaged: ",a)') "True"
write(*,*) ''
endif
etov2 = 1.0/(v2_to_E_per_amu*inputs%ab)
chan_loop: do ichan=1,spec_chords%nchan
fdens = 0.d0 ; hdens = 0.d0 ; tdens = 0.d0
halodens = 0.d0 ; dcxdens = 0.d0
plasma = plasma*0.d0
fields = fields*0.d0
wght_tot = 0.d0
mean_f = 0.d0
do k=1,beam_grid%nz
do j=1,beam_grid%ny
x_loop: do i=1,beam_grid%nx
inter = spec_chords%inter(i,j,k)
cid = 0
cind = 0
do while (cid.ne.ichan.and.cind.lt.inter%nchan)
cind = cind + 1
cid = inter%los_elem(cind)%id
enddo
if(cid.eq.ichan) then
ind = [i,j,k]
ri = [beam_grid%xc(i), beam_grid%yc(j), beam_grid%zc(k)]
call in_plasma(ri,inp)
if(.not.inp) cycle x_loop
dlength = inter%los_elem(cind)%length
fdens = fdens + neut%full(:,i,j,k)*dlength
hdens = hdens + neut%half(:,i,j,k)*dlength
tdens = tdens + neut%third(:,i,j,k)*dlength
dcxdens = dcxdens + neut%dcx(:,i,j,k)*dlength
halodens = halodens + neut%halo(:,i,j,k)*dlength
wght = (neut%full(3,i,j,k) + neut%half(3,i,j,k) + &
neut%third(3,i,j,k) + neut%dcx(3,i,j,k) + &
neut%halo(3,i,j,k))*dlength
call get_plasma(plasma_cell,pos=ri)
call get_fields(fields_cell,pos=ri)
plasma = plasma + wght*plasma_cell
fields = fields + wght*fields_cell
if (inputs%dist_type.eq.1) then
do ipitch=1,inputs%np_wght
do ienergy=1,inputs%ne_wght
call get_ep_denf(ebarr(ienergy),ptcharr(ipitch),denf,coeffs=fields_cell%b)
mean_f(ienergy,ipitch) = mean_f(ienergy,ipitch) + wght*denf
enddo
enddo
endif
wght_tot = wght_tot + wght
endif
enddo x_loop
enddo
enddo
if(wght_tot.le.0) then
if(inputs%verbose.ge.1) then
write(*,'(T4,"Skipping channel ",i5,": Neutral density is zero")') ichan
endif
cycle chan_loop
else
plasma = plasma/wght_tot
plasma%in_plasma = .True.
fields = fields/wght_tot
fields%in_plasma= .True.
mean_f = mean_f/wght_tot
if(inputs%verbose.ge.1) then
write(*,'(T4,"Channel: ",i5)') ichan
write(*,'(T4,"Radius: ",f7.2)') spec_chords%radius(ichan)
write(*,'(T4,"Mean Fast-ion Density: ",ES14.5)') sum(mean_f)*dEdP
write(*,*)''
endif
endif
ri = plasma%pos
vp = ri - spec_chords%los(ichan)%lens
vnbi_f = ri - nbi%src
vnbi_f = vnbi_f/norm2(vnbi_f)*nbi%vinj
vnbi_h = vnbi_f/sqrt(2.d0)
vnbi_t = vnbi_f/sqrt(3.d0)
sigma_pi = spec_chords%los(ichan)%sigma_pi
dlength = 1.d0
!$OMP PARALLEL DO schedule(guided) collapse(3) private(eb,vabs,ptch,phi,vi,vi_norm, &
!$OMP& r_enter,r_exit,length,max_dens,ind,tracks,ntrack,dt,icell,states,rates, &
!$OMP& vhalo,denn,denf,photons,ienergy,ipitch,igyro)
do ienergy=istart,inputs%ne_wght,istep
do ipitch=1,inputs%np_wght
do igyro=1,inputs%nphi_wght
eb = ebarr(ienergy)
vabs = sqrt(eb*etov2)
ptch = ptcharr(ipitch)
phi = phiarr(igyro)
call pitch_to_vec(ptch,phi,fields,vi_norm)
vi = vabs*vi_norm
call grid_intersect(ri,vi,length,r_enter,r_exit)
call track(r_enter, vi, tracks, ntrack)
max_dens = 0.d0
do icell=1,ntrack
ind = tracks(icell)%ind
tracks(icell)%flux = neut%full(3,ind(1),ind(2),ind(3)) + &
neut%half(3,ind(1),ind(2),ind(3)) + &
neut%third(3,ind(1),ind(2),ind(3)) + &
neut%dcx(3,ind(1),ind(2),ind(3)) + &
neut%halo(3,ind(1),ind(2),ind(3))
if(tracks(icell)%flux.gt.max_dens) max_dens=tracks(icell)%flux
enddo
dt = 0.d0
do icell=1,ntrack
if(tracks(icell)%flux.gt.(0.5*max_dens)) then
dt = dt + tracks(icell)%time
endif
enddo
states=0.d0
call bb_cx_rates(fdens,vi,vnbi_f,rates)
states = states + rates
call bb_cx_rates(hdens,vi,vnbi_h,rates)
states = states + rates
call bb_cx_rates(tdens,vi,vnbi_t,rates)
states = states + rates
call bt_cx_rates(plasma, dcxdens + halodens, vi, beam_ion, rates)
states = states + rates
call colrad(plasma,beam_ion,vi,dt,states,denn,photons)
denf = mean_f(ienergy,ipitch)*dEdP
photons = photons/real(inputs%nphi_wght)
call store_fw_photons_at_chan(ichan, ienergy, ipitch, &
vp, vi, fields, dlength, sigma_pi, denf, photons)
enddo
enddo
enddo
!$OMP END PARALLEL DO
enddo chan_loop
fweight%mean_f = fweight%mean_f/(dEdP)
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_los