store_npa Subroutine

subroutine store_npa(det, ri, rf, vn, flux, orbit_class, passive)
    !+ Store NPA particles in [[libfida:npa]]
    integer, intent(in)                     :: det
        !+ Detector/Channel Number
    real(Float64), dimension(3), intent(in) :: ri
        !+ Birth position in beam coordinates [cm]
    real(Float64), dimension(3), intent(in) :: rf
        !+ Detector position in beam coordinates [cm]
    real(Float64), dimension(3), intent(in) :: vn
        !+ Particle velocity [cm/s]
    real(Float64), intent(in)               :: flux
        !+ Neutral flux [neutrals/s]
    integer, intent(in), optional           :: orbit_class
        !+ Orbit class ID
    logical, intent(in), optional           :: passive
        !+ Indicates whether npa particle is passive

    integer :: iclass, oclass
    type(LocalEMFields) :: fields
    real(Float64), dimension(3) :: uvw_ri, uvw_rf,vn_norm
    real(Float64) :: energy, pitch, dE
    integer(Int32), dimension(1) :: ienergy
    type(NPAParticle), dimension(:), allocatable :: parts
    logical :: pas = .False.

    if(present(orbit_class)) then
        iclass = min(orbit_class,particles%nclass)
        oclass = orbit_class
    else
        iclass = 1
        oclass = 1
    endif

    if(present(passive)) pas = passive

    ! Convert to machine coordinates
    call xyz_to_uvw(ri,uvw_ri)
    call xyz_to_uvw(rf,uvw_rf)

    ! Calculate energy
    energy = beam_mass*v2_to_E_per_amu*dot_product(vn,vn)
    if(pas) then
        dE = pnpa%energy(2)-pnpa%energy(1)
    else
        dE = npa%energy(2)-npa%energy(1)
    endif


    ! Calculate pitch if distribution actually uses pitch
    if(inputs%dist_type.le.2) then
        call get_fields(fields, pos = ri)
        vn_norm = vn/norm2(vn)
        pitch = dot_product(fields%b_norm,vn_norm)
    else
        pitch = 0.d0
    endif

    if(pas) then
        !$OMP CRITICAL(store_npa_1)
        pnpa%npart = pnpa%npart + 1
        if(pnpa%npart.gt.pnpa%nmax) then
            pnpa%nmax = int(pnpa%nmax*2)
            allocate(parts(pnpa%nmax))
            parts(1:(pnpa%npart-1)) = pnpa%part
            deallocate(pnpa%part)
            call move_alloc(parts, pnpa%part)
        endif
        pnpa%part(pnpa%npart)%detector = det
        pnpa%part(pnpa%npart)%class = oclass
        pnpa%part(pnpa%npart)%xi = uvw_ri(1)
        pnpa%part(pnpa%npart)%yi = uvw_ri(2)
        pnpa%part(pnpa%npart)%zi = uvw_ri(3)
        pnpa%part(pnpa%npart)%xf = uvw_rf(1)
        pnpa%part(pnpa%npart)%yf = uvw_rf(2)
        pnpa%part(pnpa%npart)%zf = uvw_rf(3)
        pnpa%part(pnpa%npart)%energy = energy
        pnpa%part(pnpa%npart)%pitch = pitch
        pnpa%part(pnpa%npart)%weight = flux
        ienergy = minloc(abs(pnpa%energy - energy))
        pnpa%flux(ienergy(1),det,iclass) = &
            pnpa%flux(ienergy(1),det,iclass) + flux/dE
        !$OMP END CRITICAL(store_npa_1)
    else
        !$OMP CRITICAL(store_npa_2)
        npa%npart = npa%npart + 1
        if(npa%npart.gt.npa%nmax) then
            npa%nmax = int(npa%nmax*2)
            allocate(parts(npa%nmax))
            parts(1:(npa%npart-1)) = npa%part
            deallocate(npa%part)
            call move_alloc(parts, npa%part)
        endif
        npa%part(npa%npart)%detector = det
        npa%part(npa%npart)%class = oclass
        npa%part(npa%npart)%xi = uvw_ri(1)
        npa%part(npa%npart)%yi = uvw_ri(2)
        npa%part(npa%npart)%zi = uvw_ri(3)
        npa%part(npa%npart)%xf = uvw_rf(1)
        npa%part(npa%npart)%yf = uvw_rf(2)
        npa%part(npa%npart)%zf = uvw_rf(3)
        npa%part(npa%npart)%energy = energy
        npa%part(npa%npart)%pitch = pitch
        npa%part(npa%npart)%weight = flux
        ienergy = minloc(abs(npa%energy - energy))
        npa%flux(ienergy(1),det,iclass) = &
            npa%flux(ienergy(1),det,iclass) + flux/dE
        !$OMP END CRITICAL(store_npa_2)
    endif

end subroutine store_npa