fida_mc Subroutine

subroutine fida_mc
    !+ Calculate Active 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
    !! 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
    integer :: ntrack
    type(ParticleTrack), dimension(beam_grid%ntrack) :: tracks
    logical :: los_intersect
    integer :: jj      !! counter along track
    real(Float64) :: photons !! photon flux
    integer, dimension(5) :: neut_types=[1,2,3,4,5]
    real(Float64), dimension(3) :: uvw, uvw_vi
    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_fida)/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,uvw,uvw_vi,ntrack,jj,rates,denn,los_intersect,states,photons)
    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
        if(.not.fast_ion%beam_grid_cross_grid) cycle loop_over_fast_ions
        gamma_loop: do igamma=1,ngamma
            if(particles%axisym) then
                !! Pick random toroidal angle
                call randu(randomu)
                phi = fast_ion%beam_grid_phi_enter + fast_ion%delta_phi*randomu(1)
            else
                phi = fast_ion%phi
            endif
            s = sin(phi)
            c = cos(phi)

            !! Calculate position in machine coordinates
            uvw(1) = fast_ion%r*c
            uvw(2) = fast_ion%r*s
            uvw(3) = fast_ion%z

            !! Convert to beam grid coordinates
            call uvw_to_xyz(uvw, ri)

            if(inputs%dist_type.eq.2) then
                !! Get electomagnetic fields
                call get_fields(fields, pos=ri)

                !! 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
                uvw_vi(1) = c*fast_ion%vr - s*fast_ion%vt
                uvw_vi(2) = s*fast_ion%vr + c*fast_ion%vt
                uvw_vi(3) = fast_ion%vz
                vi = matmul(beam_grid%inv_basis,uvw_vi)
            endif

            !! Track particle through grid
            call track(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_beam_cx_rate(tracks(1)%ind,ri,vi,beam_ion,neut_types,rates)
            if(sum(rates).le.0.)cycle gamma_loop

            !! Weight CX rates by ion source density
            states=rates*fast_ion%weight*(fast_ion%delta_phi/(2*pi))/beam_grid%dv/ngamma

            !! 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_fida_photons(tracks(jj)%pos, vi, photons, fast_ion%class)
            enddo loop_along_track
        enddo gamma_loop
    enddo loop_over_fast_ions
    !$OMP END PARALLEL DO

#ifdef _MPI
    call parallel_sum(spec%fida)
#endif

end subroutine fida_mc