Calculates approximate neutral beam emission (full, half, third) from user supplied neutrals file
subroutine nbi_spec
!+ Calculates approximate neutral beam emission (full, half, third)
!+ from user supplied neutrals file
integer :: ic, i, j, k, it
real(Float64), dimension(3) :: ri, vnbi, random3, rc
integer,dimension(3) :: ind
!! Determination of the CX probability
real(Float64) :: nbif_photons, nbih_photons, nbit_photons
real(Float64) :: f_wght, h_wght, t_wght
real(Float64) :: f_tot, h_tot, t_tot
real(Float64), dimension(n_stark,inputs%nlambda,spec_chords%nchan) :: full, half, third
real(Float64), dimension(:,:,:,:), allocatable :: fullstokes, halfstokes, thirdstokes
logical :: inp
integer :: n = 10000
!$OMP PARALLEL DO schedule(dynamic,1) private(i,j,k,ic,ind, &
!$OMP& nbif_photons, nbih_photons, nbit_photons, rc, ri,inp, vnbi,&
!$OMP& random3,f_tot,h_tot,t_tot,f_wght,h_wght,t_wght,&
!$OMP& full,half,third, &
!$OMP& fullstokes, halfstokes, thirdstokes)
loop_over_cells: do ic = istart, spec_chords%ncell, istep
call ind2sub(beam_grid%dims,spec_chords%cell(ic),ind)
i = ind(1) ; j = ind(2) ; k = ind(3)
nbif_photons = neut%full%dens(3,i,j,k)*tables%einstein(2,3)
nbih_photons = neut%half%dens(3,i,j,k)*tables%einstein(2,3)
nbit_photons = neut%third%dens(3,i,j,k)*tables%einstein(2,3)
if((nbif_photons + nbih_photons + nbit_photons).le.0.0) then
cycle loop_over_cells
endif
rc = [beam_grid%xc(i), beam_grid%yc(j), beam_grid%zc(k)]
!Find a point in cell that is also in the plasma
ri = rc
call in_plasma(ri, inp)
do while (.not.inp)
call randu(random3)
ri = rc + beam_grid%dr*(random3 - 0.5)
call in_plasma(ri,inp)
enddo
allocate(fullstokes(n_stark,4,inputs%nlambda,spec_chords%nchan))
allocate(halfstokes(n_stark,4,inputs%nlambda,spec_chords%nchan))
allocate(thirdstokes(n_stark,4,inputs%nlambda,spec_chords%nchan))
f_tot = 0.0 ; h_tot = 0.0 ; t_tot = 0.0
full = 0.0 ; half = 0.0 ; third = 0.0
fullstokes = 0.0 ; halfstokes = 0.0 ; thirdstokes = 0.0
do it=1, n
!! Full Spectra
call mc_nbi_cell(ind, nbif_type, vnbi, f_wght)
f_tot = f_tot + f_wght
call store_photons(ri, vnbi, beam_lambda0, f_wght*nbif_photons, full, fullstokes)
!! Half Spectra
call mc_nbi_cell(ind, nbih_type, vnbi, h_wght)
h_tot = h_tot + h_wght
call store_photons(ri, vnbi, beam_lambda0, h_wght*nbih_photons, half, halfstokes)
!! Third Spectra
call mc_nbi_cell(ind, nbit_type, vnbi, t_wght)
t_tot = t_tot + t_wght
call store_photons(ri, vnbi, beam_lambda0, t_wght*nbit_photons, third, thirdstokes)
enddo
!$OMP CRITICAL(nbi_spec_1)
spec%full = spec%full + full/f_tot
spec%half = spec%half + half/h_tot
spec%third = spec%third + third/t_tot
spec%fullstokes = spec%fullstokes + fullstokes/f_tot
spec%halfstokes = spec%halfstokes + halfstokes/h_tot
spec%thirdstokes = spec%thirdstokes + thirdstokes/t_tot
!$OMP END CRITICAL(nbi_spec_1)
deallocate(fullstokes)
deallocate(halfstokes)
deallocate(thirdstokes)
enddo loop_over_cells
!$OMP END PARALLEL DO
#ifdef _MPI
!! Combine Spectra
call parallel_sum(spec%full)
call parallel_sum(spec%half)
call parallel_sum(spec%third)
#endif
end subroutine nbi_spec