! Copyright 2024 - The Minton Group at Purdue University
! This file is part of Swiftest.
! Swiftest is free software: you can redistribute it and/or modify it under the terms of the GNU General Public License 
! as published by the Free Software Foundation, either version 3 of the License, or (at your option) any later version.
! Swiftest is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty 
! of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details.
! You should have received a copy of the GNU General Public License along with Swiftest. 
! If not, see: https://www.gnu.org/licenses. 

submodule (encounter) s_encounter_check
   use swiftest
contains

   module subroutine encounter_check_all_plpl(param, npl, r, v, renc, dt, nenc, index1, index2, lvdotr)
      !! author: David A. Minton
      !!
      !! Check for encounters between massive bodies. Choose between the standard triangular or the Sort & Sweep method based on 
      !! user inputs
      implicit none
      ! Arguments
      class(base_parameters),                  intent(inout) :: param  
         !! Current Swiftest run configuration parameter5s
      integer(I4B),                            intent(in)    :: npl    
         !! Total number of massive bodies
      real(DP),     dimension(:),              intent(in)    :: renc   
         !! Critical radii of massive bodies that defines an encounter 
      real(DP),     dimension(:,:),            intent(in)    :: r      
         !! Position vectors of massive bodies
      real(DP),     dimension(:,:),            intent(in)    :: v      
         !! Velocity vectors of massive bodies
      real(DP),                                intent(in)    :: dt     
         !! Step size
      integer(I8B),                            intent(out)   :: nenc   
         !! Total number of encounters
      integer(I4B), dimension(:), allocatable, intent(out)   :: index1 
         !! List of indices for body 1 in each encounter
      integer(I4B), dimension(:), allocatable, intent(out)   :: index2 
         !! List of indices for body 2 in each encounter
      logical,      dimension(:), allocatable, intent(out)   :: lvdotr 
         !! Logical flag indicating the sign of v .dot. x

      if (param%lencounter_sas_plpl) then
         call encounter_check_all_sort_and_sweep_plpl(npl, r, v, renc, dt, nenc, index1, index2, lvdotr)
      else
         call encounter_check_all_triangular_plpl(npl, r, v, renc, dt, nenc, index1, index2, lvdotr) 
      end if

      return
   end subroutine encounter_check_all_plpl


   module subroutine encounter_check_all_plplm(param, nplm, nplt, rplm, vplm, rplt, vplt, rencm, renct, dt, nenc, index1, index2, &
                                               lvdotr)
      !! author: David A. Minton
      !!
      !! Check for encounters between fully interacting massive bodies partially interacting massive bodies. Choose between the
      !! standard triangular or the Sort & Sweep method based on user inputs
      implicit none
      ! Arguments
      class(base_parameters),                  intent(inout) :: param  
         !! Current Swiftest run configuration parameter5s
      integer(I4B),                            intent(in)    :: nplm   
         !! Total number of fully interacting massive bodies 
      integer(I4B),                            intent(in)    :: nplt   
         !! Total number of partially interacting masive bodies (GM < GMTINY) 
      real(DP),     dimension(:,:),            intent(in)    :: rplm   
         !! Position vectors of fully interacting massive bodies
      real(DP),     dimension(:,:),            intent(in)    :: vplm   
         !! Velocity vectors of fully interacting massive bodies
      real(DP),     dimension(:,:),            intent(in)    :: rplt   
         !! Position vectors of partially interacting massive bodies
      real(DP),     dimension(:,:),            intent(in)    :: vplt   
         !! Velocity vectors of partially interacting massive bodies
      real(DP),     dimension(:),              intent(in)    :: rencm  
         !! Critical radii of fully interacting massive bodies that defines an encounter
      real(DP),     dimension(:),              intent(in)    :: renct  
         !! Critical radii of partially interacting massive bodies that defines an encounter
      real(DP),                                intent(in)    :: dt     
         !! Step size
      integer(I8B),                            intent(out)   :: nenc   
         !! Total number of encounters
      integer(I4B), dimension(:), allocatable, intent(out)   :: index1 
         !! List of indices for body 1 in each encounter
      integer(I4B), dimension(:), allocatable, intent(out)   :: index2 
         !! List of indices for body 2 in each encounter
      logical,      dimension(:), allocatable, intent(out)   :: lvdotr 
         !! Logical flag indicating the sign of v .dot. x
      ! Internals
      logical,      dimension(:), allocatable :: plmplt_lvdotr 
         !! Logical flag indicating the sign of v .dot. x in the plm-plt group
      integer(I4B), dimension(:), allocatable :: plmplt_index1 
         !! List of indices for body 1 in each encounter in the plm-plt group
      integer(I4B), dimension(:), allocatable :: plmplt_index2 
         !! List of indices for body 2 in each encounter in the plm-lt group
      integer(I8B)                            :: plmplt_nenc   
         !! Number of encounters of the plm-plt group
      class(base_parameters), allocatable :: tmp_param     
         !! Temporary parameter structure to turn off adaptive timer for the pl-pl phase if necessary
      integer(I8B), dimension(:), allocatable :: ind
      integer(I4B), dimension(:), allocatable :: itmp
      logical, dimension(:), allocatable :: ltmp

      allocate(tmp_param, source=param)

      ! Start with the pl-pl group
      call encounter_check_all_plpl(tmp_param, nplm, rplm, vplm, rencm, dt, nenc, index1, index2, lvdotr)

      if (param%lencounter_sas_plpl) then
         call encounter_check_all_sort_and_sweep_plplm(nplm, nplt, rplm, vplm, rplt, vplt, rencm, renct, dt, &
                                                       plmplt_nenc, plmplt_index1, plmplt_index2, plmplt_lvdotr)
      else
         call encounter_check_all_triangular_plplm(nplm, nplt, rplm, vplm, rplt, vplt, rencm, renct, dt, &
                                                   plmplt_nenc, plmplt_index1, plmplt_index2, plmplt_lvdotr) 
      end if

      if (plmplt_nenc > 0) then ! Consolidate the two lists
         allocate(itmp(nenc+plmplt_nenc))
         if (nenc > 0) itmp(1:nenc) = index1(1:nenc)
         itmp(nenc+1:nenc+plmplt_nenc) = plmplt_index1(1:plmplt_nenc)
         call move_alloc(itmp, index1)
         allocate(itmp(nenc+plmplt_nenc))
         if (nenc > 0) itmp(1:nenc) = index2(1:nenc)
         itmp(nenc+1:nenc+plmplt_nenc) = plmplt_index2(1:plmplt_nenc) + nplm ! Be sure to shift these indices back to their natural 
                                                                             ! range
         call move_alloc(itmp, index2)
         allocate(ltmp(nenc+plmplt_nenc))
         if (nenc > 0) ltmp(1:nenc) = lvdotr(1:nenc)
         ltmp(nenc+1:nenc+plmplt_nenc) = plmplt_lvdotr(1:plmplt_nenc)
         call move_alloc(ltmp, lvdotr)
         nenc = nenc + plmplt_nenc

         call util_sort(index1, ind)
         call util_sort_rearrange(index1, ind, nenc)
         call util_sort_rearrange(index2, ind, nenc)
         call util_sort_rearrange(lvdotr, ind, nenc)

      end if

      return
   end subroutine encounter_check_all_plplm


   module subroutine encounter_check_all_pltp(param, npl, ntp, rpl, vpl, rtp, vtp, renc, dt, nenc, index1, index2, lvdotr)
      !! author: David A. Minton
      !!
      !! Check for encounters between massive bodies and test particles. Choose between the standard triangular or the Sort & Sweep
      !! method based on user inputs
      !!
      implicit none
      ! Arguments
      class(base_parameters),                  intent(inout) :: param  
         !! Current Swiftest run configuration parameter5s
      integer(I4B),                            intent(in)    :: npl    
         !! Total number of massive bodies 
      integer(I4B),                            intent(in)    :: ntp    
         !! Total number of test particles 
      real(DP),     dimension(:,:),            intent(in)    :: rpl    
         !! Position vectors of massive bodies
      real(DP),     dimension(:,:),            intent(in)    :: vpl    
         !! Velocity vectors of massive bodies
      real(DP),     dimension(:,:),            intent(in)    :: rtp    
         !! Position vectors of test particlse
      real(DP),     dimension(:,:),            intent(in)    :: vtp    
         !! Velocity vectors of test particles
      real(DP),     dimension(:),              intent(in)    :: renc   
         !! Critical radii of massive bodies that defines an encounter
      real(DP),                                intent(in)    :: dt     
         !! Step size
      integer(I8B),                            intent(out)   :: nenc   
         !! Total number of encounters
      integer(I4B), dimension(:), allocatable, intent(out)   :: index1 
         !! List of indices for body 1 in each encounter
      integer(I4B), dimension(:), allocatable, intent(out)   :: index2 
         !! List of indices for body 2 in each encounter
      logical,      dimension(:), allocatable, intent(out)   :: lvdotr 
         !! Logical flag indicating the sign of v .dot. x

      if (param%lencounter_sas_pltp) then
         call encounter_check_all_sort_and_sweep_pltp(npl, ntp, rpl, vpl, rtp, vtp, renc, dt, nenc, index1, index2, lvdotr)
      else
         call encounter_check_all_triangular_pltp(npl, ntp, rpl, vpl, rtp, vtp, renc, dt, nenc, index1, index2, lvdotr) 
      end if

      return
   end subroutine encounter_check_all_pltp


   subroutine encounter_check_all_sort_and_sweep_plpl(npl, r, v, renc, dt, nenc, index1, index2, lvdotr)
      !! author: David A. Minton
      !!
      !! Check for encounters between massive bodies. 
      !! This is the sort and sweep version
      !! References: Adapted from Christer Ericson's _Real Time Collision Detection_
      !!
      implicit none
      ! Arguments
      integer(I4B),                            intent(in)  :: npl    
         !! Total number of massive bodies
      real(DP),     dimension(:,:),            intent(in)  :: r      
         !! Position vectors of massive bodies
      real(DP),     dimension(:,:),            intent(in)  :: v      
         !! Velocity vectors of massive bodies
      real(DP),     dimension(:),              intent(in)  :: renc   
         !! Critical radii of massive bodies that defines an encounter 
      real(DP),                                intent(in)  :: dt     
         !! Step size
      integer(I8B),                            intent(out) :: nenc   
         !! Total number of encounters
      integer(I4B), dimension(:), allocatable, intent(out) :: index1 
         !! List of indices for body 1 in each encounter
      integer(I4B), dimension(:), allocatable, intent(out) :: index2 
         !! List of indices for body 2 in each encounter
      logical,      dimension(:), allocatable, intent(out) :: lvdotr 
         !! Logical flag indicating the sign of v .dot. x
      ! Internals
      integer(I4B) :: i, n
      integer(I4B), save :: npl_last = 0
      type(encounter_bounding_box), save :: boundingbox
      real(DP), dimension(npl) :: rmin,rmax
      real(DP) :: rmag

      if (npl == 0) return

      ! If this is the first time through, build the index lists
      n = 2 * npl
      if (npl_last /= npl) then
         call boundingbox%setup(npl, npl_last)
         npl_last = npl
      end if

#ifdef DOCONLOC
      do concurrent (i = 1:npl) shared(r,renc,rmin,rmax) local(rmag)
#else
      do concurrent (i = 1:npl)
#endif
         rmag = norm2(r(:,i))
         rmax(i) = rmag + RSWEEP_FACTOR * renc(i)
         rmin(i) = rmag - RSWEEP_FACTOR * renc(i)
      end do

      call boundingbox%aabb%sort(npl, [rmin,rmax]) 

      call boundingbox%sweep(npl, r, v, renc, dt, nenc, index1, index2, lvdotr) 

      return
   end subroutine encounter_check_all_sort_and_sweep_plpl


   subroutine encounter_check_all_sort_and_sweep_plplm(nplm, nplt, rplm, vplm, rplt, vplt, rencm, renct, dt, nenc, index1, index2, &
                                                      lvdotr)
      !! author: David A. Minton
      !!
      !! Check for encounters between massive bodies and test particles. 
      !! This is the sort and sweep version
      !! References: Adapted from Christer Ericson's _Real Time Collision Detection_
      !!
      implicit none
      ! Arguments
      integer(I4B),                            intent(in)  :: nplm   
         !! Total number of fully interacting massive bodies 
      integer(I4B),                            intent(in)  :: nplt   
         !! Total number of partially interacting masive bodies (GM < GMTINY) 
      real(DP),     dimension(:,:),            intent(in)  :: rplm   
         !! Position vectors of fully interacting massive bodies
      real(DP),     dimension(:,:),            intent(in)  :: vplm   
         !! Velocity vectors of fully interacting massive bodies
      real(DP),     dimension(:,:),            intent(in)  :: rplt   
         !! Position vectors of partially interacting massive bodies
      real(DP),     dimension(:,:),            intent(in)  :: vplt   
         !! Velocity vectors of partially interacting massive bodies
      real(DP),     dimension(:),              intent(in)  :: rencm  
         !! Critical radii of fully interacting massive bodies that defines an encounter
      real(DP),     dimension(:),              intent(in)  :: renct  
         !! Critical radii of partially interacting massive bodies that defines an encounter
      real(DP),                                intent(in)  :: dt     
         !! Step size
      integer(I8B),                            intent(out) :: nenc   
         !! Total number of encounter
      integer(I4B), dimension(:), allocatable, intent(out) :: index1 
         !! List of indices for body 1 in each encounter
      integer(I4B), dimension(:), allocatable, intent(out) :: index2 
         !! List of indices for body 2 in each encounter
      logical,      dimension(:), allocatable, intent(out) :: lvdotr 
         !! Logical flag indicating the sign of v .dot. x
      ! Internals
      type(encounter_bounding_box), save :: boundingbox
      integer(I4B) :: i, n, ntot
      integer(I4B), save :: ntot_last = 0
      real(DP), dimension(nplm+nplt) :: rmin,rmax
      real(DP) :: rmag

      ! If this is the first time through, build the index lists
      if ((nplm == 0) .or. (nplt == 0)) return

      ntot = nplm + nplt
      n = 2 * ntot
      if (ntot /= ntot_last) then
         call boundingbox%setup(ntot, ntot_last)
         ntot_last = ntot
      end if

#ifdef DOCONLOC
      do concurrent (i = 1:nplm) shared(rmin,rmax,rplm,rencm) local(rmag)
#else
      do concurrent (i = 1:nplm)
#endif
         rmag = norm2(rplm(:,i))
         rmax(i) = rmag + RSWEEP_FACTOR * rencm(i)
         rmin(i) = rmag - RSWEEP_FACTOR * rencm(i)
      end do
#ifdef DOCONLOC
      do concurrent (i = 1:nplt) shared(rmin,rmax,rplt,renct) local(rmag)
#else
      do concurrent (i = 1:nplt)
#endif
         rmag = norm2(rplt(:,i))
         rmax(nplm+i) = rmag + RSWEEP_FACTOR * renct(i)
         rmin(nplm+i) = rmag - RSWEEP_FACTOR * renct(i)
      end do

      call boundingbox%aabb%sort(ntot, [rmin, rmax])

      call boundingbox%sweep(nplm, nplt, rplm, vplm, rplt, vplt, rencm, renct, dt, nenc, index1, index2, lvdotr) 

      return
   end subroutine encounter_check_all_sort_and_sweep_plplm


   subroutine encounter_check_all_sort_and_sweep_pltp(npl, ntp, rpl, vpl, rtp, vtp, rencpl, dt, nenc, index1, index2, lvdotr)
      !! author: David A. Minton
      !!
      !! Check for encounters between massive bodies and test particles. 
      !! This is the sort and sweep version
      !! References: Adapted from Christer Ericson's _Real Time Collision Detection_
      !!
      implicit none
      ! Arguments
      integer(I4B),                            intent(in)  :: npl    
         !! Total number of massive bodies 
      integer(I4B),                            intent(in)  :: ntp    
         !! Total number of test particles 
      real(DP),     dimension(:,:),            intent(in)  :: rpl    
         !! Position vectors of massive bodies
      real(DP),     dimension(:,:),            intent(in)  :: vpl    
         !! Velocity vectors of massive bodies
      real(DP),     dimension(:,:),            intent(in)  :: rtp    
         !! Position vectors of massive bodies
      real(DP),     dimension(:,:),            intent(in)  :: vtp    
         !! Velocity vectors of massive bodies
      real(DP),     dimension(:),              intent(in)  :: rencpl 
         !! Critical radii of massive bodies that defines an encounter
      real(DP),                                intent(in)  :: dt     
         !! Step size
      integer(I8B),                            intent(out) :: nenc   
         !! Total number of encounter
      integer(I4B), dimension(:), allocatable, intent(out) :: index1 
         !! List of indices for body 1 in each encounter
      integer(I4B), dimension(:), allocatable, intent(out) :: index2 
         !! List of indices for body 2 in each encounter
      logical,      dimension(:), allocatable, intent(out) :: lvdotr 
         !! Logical flag indicating the sign of v .dot. x
      ! Internals
      type(encounter_bounding_box), save :: boundingbox
      integer(I4B) :: i, n, ntot
      integer(I4B), save :: ntot_last = 0
      real(DP), dimension(npl+ntp) :: rmin,rmax
      real(DP), dimension(ntp) :: renctp
      real(DP) :: rmag

      ! If this is the first time through, build the index lists
      if ((ntp == 0) .or. (npl == 0)) return

      ntot = npl + ntp
      n = 2 * ntot
      if (ntot /= ntot_last) then
         call boundingbox%setup(ntot, ntot_last)
         ntot_last = ntot
      end if

      renctp(:) = 0.0_DP

#ifdef DOCONLOC
      do concurrent (i = 1:npl) shared(rmin,rmax,rpl,rencpl) local(rmag)
#else
      do concurrent (i = 1:npl)
#endif
         rmag = norm2(rpl(:,i))
         rmax(i) = rmag + RSWEEP_FACTOR * rencpl(i)
         rmin(i) = rmag - RSWEEP_FACTOR * rencpl(i)
      end do
#ifdef DOCONLOC
      do concurrent (i = 1:ntp) shared(rmin,rmax,rtp,renctp) local(rmag)
#else
      do concurrent (i = 1:ntp) 
#endif
         rmag = norm2(rtp(:,i))
         rmax(npl+i) = rmag + RSWEEP_FACTOR * renctp(i)
         rmin(npl+i) = rmag - RSWEEP_FACTOR * renctp(i)
      end do

      call boundingbox%aabb%sort(ntot, [rmin, rmax])

      call boundingbox%sweep(npl, ntp, rpl, vpl, rtp, vtp, rencpl, renctp, dt, nenc, index1, index2, lvdotr) 

      return
   end subroutine encounter_check_all_sort_and_sweep_pltp


   pure subroutine encounter_check_all_sweep_one(i, n, xi, yi, zi, vxi, vyi, vzi, x, y, z, vx, vy, vz, renci, renc, dt, &
                                                 ind_arr, lgood, nenci, index1, index2, lvdotr)
      !! author: David A. Minton
      !!
      !! Check for encounters between the ith body and all other bodies.
      !! This is used in the narrow phase of the sort & sweep algorithm
      implicit none
      ! Arguments
      integer(I4B),                            intent(in)    :: i             
         !! Index of the ith body that is being checked
      integer(I4B),                            intent(in)    :: n             
         !! Total number of bodies being checked
      real(DP),                                intent(in)    :: xi, yi, zi    
         !! Position vector components of the ith body
      real(DP),                                intent(in)    :: vxi, vyi, vzi 
         !! Velocity vector components of the ith body
      real(DP),     dimension(:),              intent(in)    :: x, y, z       
         !! Arrays of position vector components of all bodies
      real(DP),     dimension(:),              intent(in)    :: vx, vy, vz    
         !! Arrays of velocity vector components of all bodies
      real(DP),                                intent(in)    :: renci         
         !! Encounter radius of the ith body
      real(DP),     dimension(:),              intent(in)    :: renc          
         !! Array of encounter radii of all bodies
      real(DP),                                intent(in)    :: dt            
         !! Step size
      integer(I4B), dimension(:),              intent(in)    :: ind_arr       
         !! Index array [1, 2, ..., n]
      logical,      dimension(:),              intent(in)    :: lgood         
         !! Logical array mask where true values correspond to bodies selected in the broad phase
      integer(I8B),                            intent(out)   :: nenci         
         !! Total number of encountering bodies
      integer(I4B), dimension(:), allocatable, intent(inout) :: index1        
         !! Array of indices of the ith body of size nenci [i, i, ..., i]
      integer(I4B), dimension(:), allocatable, intent(inout) :: index2        
         !! Array of indices of the encountering bodies of size nenci 
      logical,      dimension(:), allocatable, intent(inout) :: lvdotr        
         !! v.dot.r direction array
      ! Internals
      integer(I4B) :: j
      real(DP) :: xr, yr, zr, vxr, vyr, vzr, renc12
      logical, dimension(n) :: lencounteri, lvdotri

      lencounteri(:) = .false.
#ifdef DOCONLOC
      do concurrent(j = 1:n, lgood(j)) shared(lgood,lencounteri,lvdotri,x,y,z,vx,vy,vz,renci,renc) &
                                       local(xr,yr,zr,vxr,vyr,vzr,renc12)
#else
      do concurrent(j = 1:n, lgood(j))
#endif
         xr = x(j) - xi
         yr = y(j) - yi
         zr = z(j) - zi
         vxr = vx(j) - vxi
         vyr = vy(j) - vyi
         vzr = vz(j) - vzi
         renc12 = renci + renc(j)
         call encounter_check_one(xr, yr, zr, vxr, vyr, vzr, renc12, dt, lencounteri(j), lvdotri(j))
      end do
      if (any(lencounteri(:))) then
         nenci = count(lencounteri(:))
         allocate(lvdotr(nenci), index1(nenci), index2(nenci))
         index1(:) = i
         index2(:) = pack(ind_arr(1:n), lencounteri(1:n)) 
         lvdotr(:) = pack(lvdotri(1:n), lencounteri(1:n)) 
      end if

      return
   end subroutine encounter_check_all_sweep_one


   subroutine encounter_check_all_triangular_one(i, n, xi, yi, zi, vxi, vyi, vzi, x, y, z, vx, vy, vz, renci, renc, &
                                                      dt, ind_arr, lenci)
      !! author: David A. Minton
      !!
      !! Check for encounters between the ith body and all other bodies.
      !! This is the upper triangular (double loop) version.
      implicit none
      ! Arguments
      integer(I4B),                       intent(in)  :: i             
         !! Index of the ith body that is being checked
      integer(I4B),                       intent(in)  :: n             
         !! Total number of bodies being checked
      real(DP),                           intent(in)  :: xi, yi, zi    
         !! Position vector components of the ith body
      real(DP),                           intent(in)  :: vxi, vyi, vzi 
         !! Velocity vector components of the ith body
      real(DP),             dimension(:), intent(in)  :: x, y, z       
         !! Arrays of position vector components of all bodies
      real(DP),             dimension(:), intent(in)  :: vx, vy, vz    
         !! Arrays of velocity vector components of all bodies
      real(DP),                           intent(in)  :: renci         
         !! Encounter radius of the ith body
      real(DP),             dimension(:), intent(in)  :: renc          
         !! Array of encounter radii of all bodies
      real(DP),                           intent(in)  :: dt            
         !! Step size
      integer(I4B),         dimension(:), intent(in)  :: ind_arr       
         !! Index array [1, 2, ..., n]
      class(encounter_list),              intent(out) :: lenci         
         !! Output encounter lists containing number of encounters, the v.dot.r direction array, and the index list of encountering
         !! bodies 
      ! Internals
      integer(I4B) :: j
      integer(I8B) :: nenci
      real(DP) :: xr, yr, zr, vxr, vyr, vzr, renc12
      logical, dimension(n) :: lencounteri, lvdotri

#ifdef DOCONLOC
      do concurrent(j = i+1:n) shared(lencounteri, lvdotri, renci, renc) local(xr,yr,zr,vxr,vyr,vzr,renc12)
#else
      do concurrent(j = i+1:n)
#endif
         xr = x(j) - xi
         yr = y(j) - yi
         zr = z(j) - zi
         vxr = vx(j) - vxi
         vyr = vy(j) - vyi
         vzr = vz(j) - vzi
         renc12 = renci + renc(j)
         call encounter_check_one(xr, yr, zr, vxr, vyr, vzr, renc12, dt, lencounteri(j), lvdotri(j))
      end do
      nenci = count(lencounteri(i+1:n))
      if (nenci > 0_I8B) then
         allocate(lenci%lvdotr(nenci), lenci%index1(nenci), lenci%index2(nenci))
         lenci%nenc = nenci
         lenci%index1(:) = i
         lenci%index2(:) = pack(ind_arr(i+1:n), lencounteri(i+1:n)) 
         lenci%lvdotr(:) = pack(lvdotri(i+1:n), lencounteri(i+1:n)) 
      end if

      return
   end subroutine encounter_check_all_triangular_one


   subroutine encounter_check_all_triangular_plpl(npl, r, v, renc, dt, nenc, index1, index2, lvdotr)
      !! author: David A. Minton
      !!
      !! Check for encounters between massive bodies. Split off from the main subroutine for performance
      !! This is the upper triangular (double loop) version.
      implicit none
      ! Arguments
      integer(I4B),                            intent(in)  :: npl    
         !! Total number of massive bodies
      real(DP),     dimension(:,:),            intent(in)  :: r      
         !! Position vectors of massive bodies
      real(DP),     dimension(:,:),            intent(in)  :: v      
         !! Velocity vectors of massive bodies
      real(DP),     dimension(:),              intent(in)  :: renc   
         !! Critical radii of massive bodies that defines an encounter 
      real(DP),                                intent(in)  :: dt     
         !! Step size
      integer(I8B),                            intent(out) :: nenc   
         !! Total number of encounters
      integer(I4B), dimension(:), allocatable, intent(out) :: index1 
         !! List of indices for body 1 in each encounter
      integer(I4B), dimension(:), allocatable, intent(out) :: index2 
         !! List of indices for body 2 in each encounter
      logical,      dimension(:), allocatable, intent(out) :: lvdotr 
         !! Logical flag indicating the sign of v .dot. x
      ! Internals
      integer(I4B) :: i
      integer(I4B), dimension(:), allocatable, save :: ind_arr
      type(collision_list_plpl), dimension(npl) :: lenc

      call swiftest_util_index_array(ind_arr, npl) 

      !$omp parallel do default(private) schedule(static)&
      !$omp shared(r, v, renc, lenc, ind_arr) &
      !$omp firstprivate(npl, dt)
      do i = 1,npl
         call encounter_check_all_triangular_one(i, npl, r(1,i), r(2,i), r(3,i), &
                                                         v(1,i), v(2,i), v(3,i), &
                                                         r(1,:), r(2,:), r(3,:), &
                                                         v(1,:), v(2,:), v(3,:), &
                                                         renc(i), renc(:), dt, ind_arr(:), lenc(i))
         if (lenc(i)%nenc > 0) lenc(i)%index1(:) = i
      end do
      !$omp end parallel do

      call encounter_check_collapse_ragged_list(lenc, npl, nenc, index1, index2, lvdotr)

      return
   end subroutine encounter_check_all_triangular_plpl


   subroutine encounter_check_all_triangular_plplm(nplm, nplt, rplm, vplm, rplt, vplt, rencm, renct, dt, &
                                                   nenc, index1, index2, lvdotr)
      !! author: David A. Minton
      !!
      !! Check for encounters between massive bodies and test particles. Split off from the main subroutine for performance
      !! This is the upper triangular (double loop) version.
      implicit none
      ! Arguments
      integer(I4B),                            intent(in)  :: nplm   
         !! Total number of fully interacting massive bodies 
      integer(I4B),                            intent(in)  :: nplt   
         !! Total number of partially interacting masive bodies (GM < GMTINY) 
      real(DP),     dimension(:,:),            intent(in)  :: rplm   
         !! Position vectors of fully interacting massive bodies
      real(DP),     dimension(:,:),            intent(in)  :: vplm   
         !! Velocity vectors of fully interacting massive bodies
      real(DP),     dimension(:,:),            intent(in)  :: rplt   
         !! Position vectors of partially interacting massive bodies
      real(DP),     dimension(:,:),            intent(in)  :: vplt   
         !! Velocity vectors of partially interacting massive bodies
      real(DP),     dimension(:),              intent(in)  :: rencm  
         !! Critical radii of fully interacting massive bodies that defines an encounter
      real(DP),     dimension(:),              intent(in)  :: renct  
         !! Critical radii of partially interacting massive bodies that defines an encounter
      real(DP),                                intent(in)  :: dt     
         !! Step size
      integer(I8B),                            intent(out) :: nenc   
         !! Total number of encounters
      integer(I4B), dimension(:), allocatable, intent(out) :: index1 
         !! List of indices for body 1 in each encounter
      integer(I4B), dimension(:), allocatable, intent(out) :: index2 
         !! List of indices for body 2 in each encounter
      logical,      dimension(:), allocatable, intent(out) :: lvdotr 
         !! Logical flag indicating the sign of v .dot. x
      ! Internals
      integer(I4B) :: i
      integer(I4B), dimension(:), allocatable, save :: ind_arr
      type(collision_list_plpl), dimension(nplm) :: lenc

      call swiftest_util_index_array(ind_arr, nplt)

      !$omp parallel do default(private) schedule(dynamic)&
      !$omp shared(rplm, vplm, rplt, vplt, rencm, renct, lenc, ind_arr) &
      !$omp firstprivate(nplm, nplt, dt)
      do i = 1, nplm
         call encounter_check_all_triangular_one(0, nplt, rplm(1,i), rplm(2,i), rplm(3,i), &
                                                          vplm(1,i), vplm(2,i), vplm(3,i), &
                                                          rplt(1,:), rplt(2,:), rplt(3,:), &
                                                          vplt(1,:), vplt(2,:), vplt(3,:), &
                                                          rencm(i), renct(:), dt, ind_arr(:), lenc(i))
         if (lenc(i)%nenc > 0) lenc(i)%index1(:) = i
      end do
      !$omp end parallel do

      call encounter_check_collapse_ragged_list(lenc, nplm, nenc, index1, index2, lvdotr)

      return
   end subroutine encounter_check_all_triangular_plplm


   subroutine encounter_check_all_triangular_pltp(npl, ntp, rpl, vpl, rtp, vtp, renc, dt, &
                                                  nenc, index1, index2, lvdotr)
      !! author: David A. Minton
      !!
      !! Check for encounters between massive bodies and test particles. Split off from the main subroutine for performance
      !! This is the upper triangular (double loop) version.
      implicit none
      ! Arguments
      integer(I4B),                            intent(in)  :: npl    
         !! Total number of massive bodies 
      integer(I4B),                            intent(in)  :: ntp    
         !! Total number of test particles 
      real(DP),     dimension(:,:),            intent(in)  :: rpl    
         !! Position vectors of massive bodies
      real(DP),     dimension(:,:),            intent(in)  :: vpl    
         !! Velocity vectors of massive bodies
      real(DP),     dimension(:,:),            intent(in)  :: rtp    
         !! Position vectors of massive bodies
      real(DP),     dimension(:,:),            intent(in)  :: vtp    
         !! Velocity vectors of massive bodies
      real(DP),     dimension(:),              intent(in)  :: renc   
         !! Critical radii of massive bodies that defines an encounter
      real(DP),                                intent(in)  :: dt     
         !! Step size
      integer(I8B),                            intent(out) :: nenc   
         !! Total number of encounters
      integer(I4B), dimension(:), allocatable, intent(out) :: index1 
         !! List of indices for body 1 in each encounter
      integer(I4B), dimension(:), allocatable, intent(out) :: index2 
         !! List of indices for body 2 in each encounter
      logical,      dimension(:), allocatable, intent(out) :: lvdotr 
         !! Logical flag indicating the sign of v .dot. x
      ! Internals
      integer(I4B) :: i
      integer(I4B), dimension(:), allocatable, save :: ind_arr
      type(collision_list_pltp), dimension(npl) :: lenc
      real(DP), dimension(ntp) :: renct

      call swiftest_util_index_array(ind_arr, ntp)
      renct(:) = 0.0_DP

      !$omp parallel do default(private) schedule(dynamic)&
      !$omp shared(rpl, vpl, rtp, vtp, renc, renct, lenc, ind_arr) &
      !$omp firstprivate(npl, ntp, dt)
      do i = 1, npl
         call encounter_check_all_triangular_one(0, ntp, rpl(1,i), rpl(2,i), rpl(3,i), &
                                                         vpl(1,i), vpl(2,i), vpl(3,i), &
                                                         rtp(1,:), rtp(2,:), rtp(3,:), &
                                                         vtp(1,:), vtp(2,:), vtp(3,:), &
                                                         renc(i), renct(:), dt, ind_arr(:), lenc(i))
         if (lenc(i)%nenc > 0) lenc(i)%index1(:) = i
      end do
      !$omp end parallel do

      call encounter_check_collapse_ragged_list(lenc, npl, nenc, index1, index2, lvdotr)

      return
   end subroutine encounter_check_all_triangular_pltp


   elemental module subroutine encounter_check_one(xr, yr, zr, vxr, vyr, vzr, renc, dt, lencounter, lvdotr)
      !! author: David A. Minton
      !!
      !! Determine whether a test particle and planet are having or will have an encounter within the next time step
      !!
      !! Adapted from David E. Kaufmann's Swifter routine: rmvs_chk_ind.f90
      !! Adapted from Hal Levison's Swift routine rmvs_chk_ind.f
      implicit none
      ! Arguments
      real(DP), intent(in)  :: xr, yr, zr    
         !! Relative distance vector components
      real(DP), intent(in)  :: vxr, vyr, vzr 
         !! Relative velocity vector components
      real(DP), intent(in)  :: renc        
         !! Square of the critical encounter distance
      real(DP), intent(in)  :: dt            
         !! Step size
      logical,  intent(out) :: lencounter    
         !! Flag indicating that an encounter has occurred
      logical,  intent(out) :: lvdotr        
         !! Logical flag indicating the direction of the v .dot. r vector
      ! Internals
      real(DP) :: r2crit, r2min, r2, v2, vdotr, tmin

      r2 = xr**2 + yr**2 + zr**2
      r2crit = renc**2

      if (r2 > r2crit) then 
         vdotr = vxr * xr + vyr * yr + vzr * zr
         if (vdotr > 0.0_DP) then
            r2min = r2
         else
            v2 = vxr**2 + vyr**2 + vzr**2
            if (v2 <= VSMALL) then
               r2min = r2
            else
               tmin = -vdotr / v2
         
               if (tmin < dt) then
                  r2min = r2 - vdotr**2 / v2
               else
                  r2min = r2 + 2 * vdotr * dt + v2 * dt**2
               end if
            end if
         end if
      else
         vdotr = -1.0_DP
         r2min = r2
      end if

      lvdotr = (vdotr < 0.0_DP)
      lencounter = lvdotr .and. (r2min <= r2crit) 

      return
   end subroutine encounter_check_one


   module subroutine encounter_check_collapse_ragged_list(ragged_list, n1, nenc, index1, index2, lvdotr)
      !! author: David A. Minton
      !!    
      !! Collapses a ragged index list (one encounter list per body) into a pair of index arrays and a vdotr logical array (optional)
      implicit none
      ! Arguments
      class(encounter_list), dimension(:),              intent(in)            :: ragged_list 
         !! The ragged encounter list
      integer(I4B),                                     intent(in)            :: n1          
         !! Number of bodies 1
      integer(I8B),                                     intent(out)           :: nenc        
         !! Total number of encountersj 
      integer(I4B),          dimension(:), allocatable, intent(out)           :: index1      
         !! Array of indices for body 1
      integer(I4B),          dimension(:), allocatable, intent(out)           :: index2      
         !! Array of indices for body 1
      logical,               dimension(:), allocatable, intent(out), optional :: lvdotr      
         !! Array indicating which bodies are approaching
      ! Internals
      integer(I4B) :: i
      integer(I8B) :: j1, j0, nenci
      integer(I8B), dimension(n1) :: ibeg

      associate(nenc_arr => ragged_list(:)%nenc)
         nenc = sum(nenc_arr(:))
      end associate
      if (nenc == 0) return

      allocate(index1(nenc))
      allocate(index2(nenc))
      if (present(lvdotr)) allocate(lvdotr(nenc))

      j0 = 1
      do i = 1, n1
         nenci = ragged_list(i)%nenc
         if (nenci == 0) cycle
         ibeg(i) = j0
         j0 = j0 + nenci
      end do

      !$omp parallel do simd default(private) schedule(simd: static)&
      !$omp shared(ragged_list, index1, index2, ibeg, lvdotr) &
      !$omp firstprivate(n1)
      do i = 1,n1
         if (ragged_list(i)%nenc == 0_I8B) cycle
         nenci = ragged_list(i)%nenc
         j0 = ibeg(i)
         j1 = j0 + nenci - 1_I8B
         index1(j0:j1) = ragged_list(i)%index1(:)
         index2(j0:j1) = ragged_list(i)%index2(:)
         if (present(lvdotr)) lvdotr(j0:j1) = ragged_list(i)%lvdotr(:)
      end do
      !$omp end parallel do simd

      return
   end subroutine encounter_check_collapse_ragged_list


   subroutine encounter_check_remove_duplicates(n, nenc, index1, index2, lvdotr)
      !! author: David A. Minton
      !!
      !! Takes the candidate encounter lists that came out of the sort & sweep method and remove any duplicates.
      implicit none
      ! Arguments
      integer(I4B),                            intent(in)    :: n          
         !! Number of bodies 
      integer(I8B),                            intent(inout) :: nenc       
         !! Number of encountering bodies (input is the broad phase value, output is the final narrow phase value)
      integer(I4B), dimension(:), allocatable, intent(inout) :: index1     
         !! List of indices for body 1 in each encounter
      integer(I4B), dimension(:), allocatable, intent(inout) :: index2     
         !! List of indices for body 2 in each encounter
      logical,      dimension(:), allocatable, intent(inout) :: lvdotr     
         !! Logical flag indicating the sign of v .dot. x
      ! Internals
      integer(I4B) :: i, i0
      integer(I8B) :: j, k, klo, khi, nenci
      integer(I4B), dimension(:), allocatable :: itmp
      integer(I8B), dimension(:), allocatable :: ind
      integer(I8B), dimension(n) :: ibeg, iend
      logical, dimension(:), allocatable :: ltmp
      logical, dimension(nenc) :: lencounter

      if (nenc == 0) then
         if (allocated(index1)) deallocate(index1)
         if (allocated(index2)) deallocate(index2)
         if (allocated(lvdotr)) deallocate(lvdotr)
         return
      end if

      call util_sort(index1, ind)
      call util_sort_rearrange(index1, ind, nenc)
      call util_sort_rearrange(index2, ind, nenc)
      call util_sort_rearrange(lvdotr, ind, nenc)

      ! Get the bounds on the bodies in the first index
      ibeg(:) = n
      iend(:) = 1_I8B
      i0 = index1(1) 
      ibeg(i0) = 1_I8B
      do k = 2_I8B, nenc 
         i = index1(k)
         if (i /= i0) then
            iend(i0) = k - 1_I8B
            ibeg(i) = k
            i0 = i
         end if
         if (k == nenc) iend(i) = k
      end do

      lencounter(:) = .true.
      ! Sort on the second index and remove duplicates 
      if (allocated(itmp)) deallocate(itmp)
      allocate(itmp, source=index2)
#ifdef DOCONLOC
      do concurrent(i = 1:n, iend(i) - ibeg(i) > 0_I8B) shared(iend,ibeg,index2,lencounter,itmp) local(klo,khi,nenci,j)
#else
      do concurrent(i = 1:n, iend(i) - ibeg(i) > 0_I8B)
#endif
         klo = ibeg(i)
         khi = iend(i)
         nenci = khi - klo + 1_I8B
         if (allocated(ind)) deallocate(ind)
         call util_sort(index2(klo:khi), ind)
         index2(klo:khi) = itmp(klo - 1_I8B + ind(:))
         do j = klo + 1_I8B, khi
            if (index2(j) == index2(j - 1_I8B)) lencounter(j) = .false. 
         end do
      end do

      if (count(lencounter(:)) == nenc) return
      nenc = count(lencounter(:)) ! Count the true number of encounters

      if (allocated(itmp)) deallocate(itmp)
      allocate(itmp(nenc))
      itmp(:) = pack(index1(:), lencounter(:))
      call move_alloc(itmp, index1)

      allocate(itmp(nenc))
      itmp(:) = pack(index2(:), lencounter(:))
      call move_alloc(itmp, index2)

      allocate(ltmp(nenc))
      ltmp(:) = pack(lvdotr(:), lencounter(:))
      call move_alloc(ltmp, lvdotr)

      return
   end subroutine encounter_check_remove_duplicates


   pure module subroutine encounter_check_sort_aabb_1D(self, n, extent_arr)
      !! author: David A. Minton
      !!
      !! Sorts the bounding box extents along a single dimension prior to the sweep phase. 
      !! This subroutine sets the sorted index array (ind) and the beginning/ending index list (beg & end)
      implicit none
      ! Arguments
      class(encounter_bounding_box_1D), intent(inout) :: self       
         !! Bounding box structure along a single dimension
      integer(I4B),                     intent(in)    :: n          
         !! Number of bodies with extents
      real(DP), dimension(:),           intent(in)    :: extent_arr 
         !! Array of extents of size 2*n
      ! Internals
      integer(I8B) :: i, k

      call util_sort(extent_arr, self%ind)

#ifdef DOCONLOC
      do concurrent(k = 1_I8B:2_I8B * n) shared(self,n) local(i)
#else
      do concurrent(k = 1_I8B:2_I8B * n)
#endif
         i = self%ind(k)
         if (i <= n) then
            self%ibeg(i) = k
         else
            self%iend(i - n) = k
         end if
      end do

      return
   end subroutine encounter_check_sort_aabb_1D 


   module subroutine encounter_check_sweep_aabb_double_list(self, n1, n2, r1, v1, r2, v2, renc1, renc2, dt, &
                                                            nenc, index1, index2, lvdotr)
      !! author: David A. Minton
      !!
      !! Sweeps the sorted bounding box extents and returns the true encounters (combines broad and narrow phases)
      !! Double list version (e.g. pl-tp or plm-plt)
      implicit none
      ! Arguments
      class(encounter_bounding_box),           intent(inout) :: self       
         !! Multi-dimensional bounding box structure
      integer(I4B),                            intent(in)    :: n1         
         !! Number of bodies 1
      integer(I4B),                            intent(in)    :: n2         
         !! Number of bodies 2
      real(DP),     dimension(:,:),            intent(in)    :: r1, v1     
         !! Array of position and velocity vectorrs for bodies 1
      real(DP),     dimension(:,:),            intent(in)    :: r2, v2     
         !! Array of position and velocity vectorrs for bodies 2
      real(DP),     dimension(:),              intent(in)    :: renc1      
         !! Radius of encounter regions of bodies 1
      real(DP),     dimension(:),              intent(in)    :: renc2      
         !! Radius of encounter regions of bodies 2
      real(DP),                                intent(in)    :: dt         
         !! Step size
      integer(I8B),                            intent(out)   :: nenc       
         !! Total number of encounter candidates
      integer(I4B), dimension(:), allocatable, intent(out)   :: index1     
         !! List of indices for body 1 in each encounter candidate pair
      integer(I4B), dimension(:), allocatable, intent(out)   :: index2     
         !! List of indices for body 2 in each encounter candidate pair
      logical,      dimension(:), allocatable, intent(out)   :: lvdotr     
         !! Logical array indicating which pairs are approaching
      ! Internals
      integer(I4B) :: ii, i, ntot, nbox, dim
      logical, dimension(n1+n2) :: loverlap
      logical, dimension(2*(n1+n2)) :: llist1
      integer(I4B), dimension(2*(n1+n2)) :: ext_ind
      type(collision_list_pltp), dimension(n1+n2) :: lenc         
         !! Array of encounter lists (one encounter list per body)
      integer(I4B), dimension(:), allocatable, save :: ind_arr
      integer(I8B) :: ibeg, iend
      real(DP), dimension(2*(n1+n2)) :: xind, yind, zind, vxind, vyind, vzind, rencind

      ntot = n1 + n2
      call swiftest_util_index_array(ind_arr, ntot)

      loverlap(:) = (self%aabb%ibeg(:) + 1_I8B) < (self%aabb%iend(:) - 1_I8B)
      where(self%aabb%ind(:) > ntot)
         ext_ind(:) = self%aabb%ind(:) - ntot
      elsewhere
         ext_ind(:) = self%aabb%ind(:)
      endwhere
      llist1(:) = ext_ind(:) <= n1
      where(.not.llist1(:)) ext_ind(:) = ext_ind(:) - n1


      dim = 1
      where(llist1(:))
         xind(:) = r1(1,ext_ind(:))
         yind(:) = r1(2,ext_ind(:))
         zind(:) = r1(3,ext_ind(:))
         vxind(:) = v1(1,ext_ind(:))
         vyind(:) = v1(2,ext_ind(:))
         vzind(:) = v1(3,ext_ind(:))
         rencind(:) = renc1(ext_ind(:))
      elsewhere
         xind(:) = r2(1,ext_ind(:))
         yind(:) = r2(2,ext_ind(:))
         zind(:) = r2(3,ext_ind(:))
         vxind(:) = v2(1,ext_ind(:))
         vyind(:) = v2(2,ext_ind(:))
         vzind(:) = v2(3,ext_ind(:))
         rencind(:) = renc2(ext_ind(:))
      endwhere

      where(.not.loverlap(:)) lenc(:)%nenc = 0
      !$omp parallel default(private) &
      !$omp shared(self, ext_ind, lenc, loverlap, r1, v1, r2, v2, renc1, renc2, xind, yind, zind, vxind, vyind, vzind, rencind, &
      !$omp        llist1) &
      !$omp firstprivate(ntot, n1, n2, dt, dim) 
     
      ! Do the first group of bodies (i is in list 1, all the others are from list 2)
      !$omp do schedule(static)
      do i = 1, n1
         if (loverlap(i)) then
            ibeg =  self%aabb%ibeg(i) + 1_I8B
            iend =  self%aabb%iend(i) - 1_I8B
            nbox = int(iend - ibeg, kind=I4B) + 1
            call encounter_check_all_sweep_one(i, nbox, r1(1,i), r1(2,i), r1(3,i), v1(1,i), v1(2,i), v1(3,i), &
                                                         xind(ibeg:iend), yind(ibeg:iend), zind(ibeg:iend),&
                                                         vxind(ibeg:iend), vyind(ibeg:iend), vzind(ibeg:iend), &
                                                         renc1(i), rencind(ibeg:iend), dt, ext_ind(ibeg:iend), &
                                                         .not.llist1(ibeg:iend), lenc(i)%nenc, lenc(i)%index1, lenc(i)%index2, &
                                                         lenc(i)%lvdotr)
         end if
      end do
      !$omp end do nowait

      ! Do the second group of bodies (i is in list 2, all the others are in list 1)
      !$omp do schedule(static)
      do i = n1+1, ntot
         if (loverlap(i)) then
            ibeg =  self%aabb%ibeg(i) + 1_I8B
            iend =  self%aabb%iend(i) - 1_I8B
            nbox = int(iend - ibeg, kind=I4B) + 1
            ii = i - n1
            call encounter_check_all_sweep_one(ii, nbox, r2(1,ii), r2(2,ii), r2(3,ii), v2(1,ii), v2(2,ii), v2(3,ii), &
                                                          xind(ibeg:iend), yind(ibeg:iend), zind(ibeg:iend),&
                                                          vxind(ibeg:iend), vyind(ibeg:iend), vzind(ibeg:iend), &
                                                          renc2(ii), rencind(ibeg:iend), dt, ext_ind(ibeg:iend), &
                                                          llist1(ibeg:iend), lenc(i)%nenc, lenc(i)%index2, lenc(i)%index1, &
                                                          lenc(i)%lvdotr)
         end if
      end do
      !$omp end do nowait

      !$omp end parallel

      call encounter_check_collapse_ragged_list(lenc, ntot, nenc, index1, index2, lvdotr)

      call encounter_check_remove_duplicates(ntot, nenc, index1, index2, lvdotr)

      return
   end subroutine encounter_check_sweep_aabb_double_list


   module subroutine encounter_check_sweep_aabb_single_list(self, n, r, v, renc, dt, nenc, index1, index2, lvdotr)
      !! author: David A. Minton
      !!
      !! Sweeps the sorted bounding box extents and returns the true encounters (combines broad and narrow phases)
      !! Single list version (e.g. pl-pl)
      implicit none
      ! Arguments
      class(encounter_bounding_box),           intent(inout) :: self       
         !! Multi-dimensional bounding box structure
      integer(I4B),                            intent(in)    :: n          
         !! Number of bodies
      real(DP),     dimension(:,:),            intent(in)    :: r, v       
         !! Array of position and velocity vectors 
      real(DP),     dimension(:),              intent(in)    :: renc       
         !! Radius of encounter regions of bodies 1
      real(DP),                                intent(in)    :: dt         
         !! Step size
      integer(I8B),                            intent(out)   :: nenc       
         !! Total number of encounter candidates
      integer(I4B), dimension(:), allocatable, intent(out)   :: index1     
         !! List of indices for one body in each encounter candidate pair
      integer(I4B), dimension(:), allocatable, intent(out)   :: index2     
         !! List of indices for the other body in each encounter candidate pair
      logical,      dimension(:), allocatable, intent(out)   :: lvdotr     
         !! Logical array indicating which pairs are approaching
      ! Internals
      integer(I4B) :: i, dim, itmp, nbox
      integer(I8B) :: k
      logical, dimension(n) :: loverlap
      logical, dimension(2*n) :: lencounteri
      real(DP), dimension(2*n) :: xind, yind, zind, vxind, vyind, vzind, rencind
      integer(I4B), dimension(2*n) :: ext_ind
      type(collision_list_plpl), dimension(n) :: lenc         
         !! Array of encounter lists (one encounter list per body)
      integer(I4B), dimension(:), allocatable, save :: ind_arr
      integer(I8B) :: ibeg, iend

      call swiftest_util_index_array(ind_arr, n)
      dim = 1

      ! Sweep the intervals for each of the massive bodies along one dimension
      ! This will build a ragged pair of index lists inside of the lenc data structure
      where(self%aabb%ind(:) > n)
         ext_ind(:) = self%aabb%ind(:) - n
      elsewhere
         ext_ind(:) = self%aabb%ind(:)
      endwhere

      xind(:) = r(1,ext_ind(:))
      yind(:) = r(2,ext_ind(:))
      zind(:) = r(3,ext_ind(:))
      vxind(:) = v(1,ext_ind(:))
      vyind(:) = v(2,ext_ind(:))
      vzind(:) = v(3,ext_ind(:))
      rencind(:) = renc(ext_ind(:))

      loverlap(:) = (self%aabb%ibeg(:) + 1_I8B) < (self%aabb%iend(:) - 1_I8B)
      where(.not.loverlap(:)) lenc(:)%nenc = 0

      !$omp parallel do default(private) schedule(static)&
      !$omp shared(self, ext_ind, lenc, loverlap, r, v, renc, xind, yind, zind, vxind, vyind, vzind, rencind) &
      !$omp firstprivate(n, dt, dim) 
      do i = 1, n
         if (loverlap(i)) then
            ibeg =  self%aabb%ibeg(i) + 1_I8B
            iend =  self%aabb%iend(i) - 1_I8B
            nbox = int(iend - ibeg, kind=I4B) + 1
            lencounteri(ibeg:iend) = .true.
            call encounter_check_all_sweep_one(i, nbox, r(1,i), r(2,i), r(3,i), v(1,i), v(2,i), v(3,i), &
                                                      xind(ibeg:iend), yind(ibeg:iend), zind(ibeg:iend),&
                                                      vxind(ibeg:iend), vyind(ibeg:iend), vzind(ibeg:iend), &
                                                      renc(i), rencind(ibeg:iend), dt, ext_ind(ibeg:iend), &
                                                      lencounteri(ibeg:iend), lenc(i)%nenc, lenc(i)%index1, lenc(i)%index2, &
                                                      lenc(i)%lvdotr)
            end if
      end do
      !$omp end parallel do

      call encounter_check_collapse_ragged_list(lenc, n, nenc, index1, index2, lvdotr)

      ! By convention, we always assume that index1 < index2, and so we must swap any that are out of order
#ifdef DOCONLOC
      do concurrent(k = 1_I8B:nenc, index1(k) > index2(k)) shared(index1,index2) local(itmp)
#else
      do concurrent(k = 1_I8B:nenc, index1(k) > index2(k))
#endif
         itmp = index1(k)
         index1(k) = index2(k)
         index2(k) = itmp
      end do

      call encounter_check_remove_duplicates(n, nenc, index1, index2, lvdotr)

      return
   end subroutine encounter_check_sweep_aabb_single_list

end submodule s_encounter_check
