diff --git a/Makefile b/Makefile
index 34df2b87a..463b5883f 100644
--- a/Makefile
+++ b/Makefile
@@ -98,6 +98,11 @@ fast:
 	  ln -s $(SWIFTEST_HOME)/Makefile.Defines .; \
 	  ln -s $(SWIFTEST_HOME)/Makefile .; \
 	  make fastdir
+	cd $(SWIFTEST_HOME)/src/encounter; \
+	  rm -f Makefile.Defines Makefile; \
+	  ln -s $(SWIFTEST_HOME)/Makefile.Defines .; \
+	  ln -s $(SWIFTEST_HOME)/Makefile .; \
+	  make fastdir
 	cd $(SWIFTEST_HOME)/src/fraggle; \
 	  rm -f Makefile.Defines Makefile; \
 	  ln -s $(SWIFTEST_HOME)/Makefile.Defines .; \
@@ -232,6 +237,7 @@ clean:
 	cd $(SWIFTEST_HOME)/src/modules;  rm -f Makefile.Defines Makefile *.gc*
 	cd $(SWIFTEST_HOME)/src/discard; rm -f Makefile.Defines Makefile *.gc*
 	cd $(SWIFTEST_HOME)/src/drift;   rm -f Makefile.Defines Makefile *.gc*
+	cd $(SWIFTEST_HOME)/src/encounter;   rm -f Makefile.Defines Makefile *.gc*
 	cd $(SWIFTEST_HOME)/src/fraggle;   rm -f Makefile.Defines Makefile *.gc*
 	cd $(SWIFTEST_HOME)/src/gr;          rm -f Makefile.Defines Makefile *.gc*
 	cd $(SWIFTEST_HOME)/src/helio;    rm -f Makefile.Defines Makefile *.gc*
@@ -257,29 +263,4 @@ clean:
 	cd $(COLLRESOLVE_HOME); rm -rf autom4te.cache aux Makefile stamp-h1 configure config.status config.h config.log aclocal.m4 lib* *.in *.o *.lo cambioni2019/*.o cambioni2019/*.lo
 
 
-force:
-
-#******************************************************************************
-#
-#  Author(s)   : David E. Kaufmann
-#
-#  Revision Control System (RCS) Information
-#
-#  Source File : $RCSfile: Makefile,v $
-#  Full Path   : $Source: /d1/kaufmann/development/RCS/Makefile,v $
-#  Revision    : $Revision: 0.1 $
-#  Date        : $Date: 2003/04/15 22:56:34 $
-#  Programmer  : $Author: kaufmann $
-#  Locked By   : $Locker: kaufmann $
-#  State       : $State: Exp $
-#
-#  Modification History:
-#
-#  $Log: Makefile,v $
-#  Revision 0.1  2003/04/15 22:56:34  kaufmann
-#  Initial implementation
-#
-#
-#******************************************************************************
-
-
+force:
\ No newline at end of file
diff --git a/src/encounter/encounter.f90 b/src/encounter/encounter.f90
new file mode 100644
index 000000000..302f7b1b5
--- /dev/null
+++ b/src/encounter/encounter.f90
@@ -0,0 +1,243 @@
+submodule (swiftest_classes) s_encounter
+   use swiftest
+contains
+
+
+   module subroutine encounter_check_all_flat_plpl(nplplm, k_plpl, x, v, renc, dt, lencounter, loc_lvdotr)
+      !! author: David A. Minton
+      !!
+      !! Check for encounters between massive bodies. Split off from the main subroutine for performance.
+      !! This is the flat (single loop) version.
+      implicit none
+      ! Arguments
+      integer(I8B),                 intent(in)  :: nplplm     !! Total number of plm-pl encounters to test
+      integer(I4B), dimension(:,:), intent(in)  :: k_plpl     !! List of all pl-pl encounters
+      real(DP),     dimension(:,:), intent(in)  :: x          !! Position vectors of massive bodies
+      real(DP),     dimension(:,:), intent(in)  :: v          !! Velocity vectors of massive bodies
+      real(DP),     dimension(:),   intent(in)  :: renc      !! Hill's radii of massive bodies
+      real(DP),                     intent(in)  :: dt         !! Step size
+      logical,      dimension(:),   intent(out) :: lencounter !! Logical array indicating which pair is in an encounter state
+      logical,      dimension(:),   intent(out) :: loc_lvdotr !! Logical array indicating the sign of v .dot. x for each encounter
+      ! Internals
+      integer(I8B) :: k
+      integer(I4B) :: i, j
+      real(DP) :: xr, yr, zr, vxr, vyr, vzr, renc12
+      
+      !$omp parallel do simd default(private) schedule(static)&
+      !$omp shared(nplplm, k_plpl, x, v, renc, dt, lencounter, loc_lvdotr) &
+      !$omp lastprivate(xr, yr, zr, vxr, vyr, vzr, renc1, renc2)
+      do k = 1_I8B, nplplm
+         i = k_plpl(1, k)
+         j = k_plpl(2, k)
+         xr = x(1, j) - x(1, i)
+         yr = x(2, j) - x(2, i)
+         zr = x(3, j) - x(3, i)
+         vxr = v(1, j) - v(1, i)
+         vyr = v(2, j) - v(2, i)
+         vzr = v(3, j) - v(3, i)
+         renc12 = renc(i) + renc(j)
+         call encounter_check_one(xr, yr, zr, vxr, vyr, vzr, renc12, dt, lencounter(k), loc_lvdotr(k))
+      end do
+      !$omp end parallel do simd
+
+      return
+   end subroutine encounter_check_all_flat_plpl
+
+
+   module subroutine encounter_check_all_triangular_plpl(npl, nplm, x, v, renc, dt, lvdotr, index1, index2, nenc)
+      !! 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
+      integer(I4B),                            intent(in)  :: nplm   !! Number of fully interacting massive bodies
+      real(DP),     dimension(:,:),            intent(in)  :: x      !! 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
+      logical,      dimension(:), allocatable, intent(out) :: lvdotr !! Logical flag indicating the sign of v .dot. x
+      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
+      integer(I4B),                            intent(out) :: nenc   !! Total number of encounters
+      ! Internals
+      integer(I4B) :: i, j, nenci, j0, j1
+      real(DP) :: xr, yr, zr, vxr, vyr, vzr, renc12
+      logical, dimension(npl) :: lencounteri, lvdotri
+      integer(I4B), dimension(npl) :: ind_arr
+      type lenctype
+         logical, dimension(:), allocatable :: lvdotr
+         integer(I4B), dimension(:), allocatable :: index2
+         integer(I4B) :: nenc
+      end type
+      type(lenctype), dimension(nplm) :: lenc
+   
+      ind_arr(:) = [(i, i = 1, npl)]
+      !$omp parallel do default(private) schedule(static)&
+      !$omp shared(npl, nplm, x, v, renc, dt, lenc, ind_arr) &
+      !$omp lastprivate(xr, yr, zr, vxr, vyr, vzr, renc12, lencounteri, lvdotri)
+      do i = 1, nplm
+         do concurrent(j = i+1:npl)
+            xr = x(1, j) - x(1, i)
+            yr = x(2, j) - x(2, i)
+            zr = x(3, j) - x(3, i)
+            vxr = v(1, j) - v(1, i)
+            vyr = v(2, j) - v(2, i)
+            vzr = v(3, j) - v(3, i)
+            renc12 = renc(i) + 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:npl))
+         if (nenci > 0) then
+            allocate(lenc(i)%lvdotr(nenci), lenc(i)%index2(nenci))
+            lenc(i)%nenc = nenci
+            lenc(i)%lvdotr(:) = pack(lvdotri(i+1:npl), lencounteri(i+1:npl)) 
+            lenc(i)%index2(:) = pack(ind_arr(i+1:npl), lencounteri(i+1:npl)) 
+         end if
+      end do
+      !$omp end parallel do
+
+      associate(nenc_arr => lenc(:)%nenc)
+         nenc = sum(nenc_arr(1:nplm))
+      end associate
+      if (nenc > 0) then
+         allocate(lvdotr(nenc))
+         allocate(index1(nenc))
+         allocate(index2(nenc))
+         j0 = 1
+         do i = 1, nplm
+            if (lenc(i)%nenc > 0) then
+               j1 = j0 + lenc(i)%nenc - 1
+               lvdotr(j0:j1) = lenc(i)%lvdotr(:)
+               index1(j0:j1) = i
+               index2(j0:j1) = lenc(i)%index2(:)
+               j0 = j1 + 1
+            end if
+         end do
+      end if
+
+      return
+   end subroutine encounter_check_all_triangular_plpl
+
+
+   module subroutine encounter_check_all_triangular_pltp(npl, ntp, xpl, vpl, xtp, vtp, renc, dt, lvdotr, index1, index2, nenc)
+      !! 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)  :: xpl    !! Position vectors of massive bodies
+      real(DP),     dimension(:,:),            intent(in)  :: vpl    !! Velocity vectors of massive bodies
+      real(DP),     dimension(:,:),            intent(in)  :: xtp    !! 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
+      logical,      dimension(:), allocatable, intent(out) :: lvdotr !! Logical flag indicating the sign of v .dot. x
+      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
+      integer(I4B),                            intent(out) :: nenc   !! Total number of encounters
+      ! Internals
+      integer(I4B) :: i, j, nenci, j0, j1
+      real(DP) :: xr, yr, zr, vxr, vyr, vzr, renc1, renc2
+      logical, dimension(ntp) :: lencounteri, lvdotri
+      integer(I4B), dimension(ntp) :: ind_arr
+      type lenctype
+         logical, dimension(:), allocatable :: lvdotr
+         integer(I4B), dimension(:), allocatable :: index2
+         integer(I4B) :: nenc
+      end type
+      type(lenctype), dimension(npl) :: lenc
+
+
+      ind_arr(:) = [(i, i = 1, ntp)]
+      !$omp parallel do default(private) schedule(static)&
+      !$omp shared(npl, ntp, xpl, vpl, xtp, vtp, renc, dt, lenc, ind_arr) &
+      !$omp lastprivate(xr, yr, zr, vxr, vyr, vzr, lencounteri, lvdotri)
+      do i = 1, npl
+         do concurrent(j = 1:npl)
+            xr = xtp(1, j) - xpl(1, i)
+            yr = xtp(2, j) - xpl(2, i)
+            zr = xtp(3, j) - xpl(3, i)
+            vxr = vtp(1, j) - vpl(1, i)
+            vyr = vtp(2, j) - vpl(2, i)
+            vzr = vtp(3, j) - vpl(3, i)
+            call encounter_check_one(xr, yr, zr, vxr, vyr, vzr, renc(i), dt, lencounteri(j), lvdotri(j))
+         end do
+         nenci = count(lencounteri(1:ntp))
+         if (nenci > 0) then
+            allocate(lenc(i)%lvdotr(nenci), lenc(i)%index2(nenci))
+            lenc(i)%nenc = nenci
+            lenc(i)%lvdotr(:) = pack(lvdotri(1:ntp), lencounteri(1:ntp)) 
+            lenc(i)%index2(:) = pack(ind_arr(1:ntp), lencounteri(1:ntp)) 
+         end if
+      end do
+      !$omp end parallel do
+
+      associate(nenc_arr => lenc(:)%nenc)
+         nenc = sum(nenc_arr(1:npl))
+      end associate
+      if (nenc > 0) then
+         allocate(lvdotr(nenc))
+         allocate(index1(nenc))
+         allocate(index2(nenc))
+         j0 = 1
+         do i = 1, npl
+            if (lenc(i)%nenc > 0) then
+               j1 = j0 + lenc(i)%nenc - 1
+               lvdotr(j0:j1) = lenc(i)%lvdotr(:)
+               index1(j0:j1) = i
+               index2(j0:j1) = lenc(i)%index2(:)
+               j0 = j1 + 1
+            end if
+         end do
+      end if
+
+      return
+   end subroutine encounter_check_all_triangular_pltp
+
+
+   module pure subroutine encounter_check_one(xr, yr, zr, vxr, vyr, vzr, dt, renc, lencounter, lvdotr)
+      !$omp declare simd(encounter_check_one)
+      !! 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: encounter_check_one.f90
+      !! Adapted from Hal Levison's Swift routine encounter_check_one.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)  :: dt            !! Step size
+      real(DP), intent(in)  :: renc        !! Square of the critical encounter distance
+      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
+
+      r2 = xr**2 + yr**2 + zr**2
+      r2crit = renc**2
+      lencounter = (r2 < r2crit) 
+      if (lencounter) return 
+
+      vdotr = vxr * xr + vyr * yr + vzr * zr
+      lvdotr = (vdotr < 0.0_DP)
+      if (.not.lvdotr) return
+     
+      v2 = vxr**2 + vyr**2 + vzr**2
+
+      if (-vdotr < v2 * dt) then
+         r2min = r2 - vdotr**2 / v2
+      else
+         r2min = r2 + 2 * vdotr * dt + v2 * dt**2
+      end if
+      lencounter = (r2min <= r2crit)
+
+      return
+   end subroutine encounter_check_one
+
+end submodule s_encounter
diff --git a/src/modules/rmvs_classes.f90 b/src/modules/rmvs_classes.f90
index 2fa6452f9..ee4d7ca1c 100644
--- a/src/modules/rmvs_classes.f90
+++ b/src/modules/rmvs_classes.f90
@@ -103,17 +103,6 @@ module rmvs_classes
    end type rmvs_pl
 
    interface
-      module pure subroutine rmvs_chk_ind(xr, yr, zr, vxr, vyr, vzr, dt, r2crit, lencounter, lvdotr)
-         !$omp declare simd(rmvs_chk_ind)
-         implicit none
-         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)  :: dt            !! Step size
-         real(DP), intent(in)  :: r2crit        !! Square of the critical encounter distance
-         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
-      end subroutine rmvs_chk_ind
-
       module subroutine rmvs_discard_tp(self, system, param)
          use swiftest_classes, only : swiftest_nbody_system, swiftest_parameters
          implicit none
diff --git a/src/modules/swiftest_classes.f90 b/src/modules/swiftest_classes.f90
index cd68e0059..ef6ad5f2c 100644
--- a/src/modules/swiftest_classes.f90
+++ b/src/modules/swiftest_classes.f90
@@ -310,6 +310,7 @@ module swiftest_classes
       real(DP),     dimension(:),   allocatable :: mass    !! Body mass (units MU)
       real(DP),     dimension(:),   allocatable :: Gmass   !! Mass gravitational term G * mass (units GU * MU)
       real(DP),     dimension(:),   allocatable :: rhill   !! Body mass (units MU)
+      real(DP),     dimension(:),   allocatable :: renc    !! Critical radius for close encounters
       real(DP),     dimension(:),   allocatable :: radius  !! Body radius (units DU)
       real(DP),     dimension(:,:), allocatable :: xbeg    !! Position at beginning of step
       real(DP),     dimension(:,:), allocatable :: xend    !! Position at end of step
@@ -344,6 +345,9 @@ module swiftest_classes
       procedure :: set_beg_end  => util_set_beg_end_pl    !! Sets the beginning and ending positions and velocities of planets.
       procedure :: set_mu       => util_set_mu_pl         !! Method used to construct the vectorized form of the central body mass
       procedure :: set_rhill    => util_set_rhill         !! Calculates the Hill's radii for each body
+      procedure :: set_renc_I4B => util_set_renc_I4B      !! Sets the critical radius for encounter given an inpput integer scale factor
+      procedure :: set_renc_DP  => util_set_renc_DP       !! Sets the critical radius for encounter given an input real scale factor
+      generic   :: set_renc     => set_renc_I4B, set_renc_DP 
       procedure :: sort         => util_sort_pl           !! Sorts body arrays by a sortable component
       procedure :: rearrange    => util_sort_rearrange_pl !! Rearranges the order of array elements of body based on an input index array. Used in sorting methods
       procedure :: spill        => util_spill_pl          !! "Spills" bodies from one object to another depending on the results of a mask (uses the PACK intrinsic)
@@ -566,35 +570,58 @@ module pure elemental subroutine drift_one(mu, px, py, pz, vx, vy, vz, dt, iflag
          integer(I4B), intent(out)      :: iflag !! iflag : error status flag for Danby drift (0 = OK, nonzero = ERROR)
       end subroutine drift_one
 
-      module pure subroutine util_flatten_eucl_ij_to_k(n, i, j, k)
-         !$omp declare simd(util_flatten_eucl_ij_to_k)
-         implicit none
-         integer(I4B), intent(in)  :: n !! Number of bodies
-         integer(I4B), intent(in)  :: i !! Index of the ith body
-         integer(I4B), intent(in)  :: j !! Index of the jth body
-         integer(I8B), intent(out) :: k !! Index of the flattened matrix
-      end subroutine util_flatten_eucl_ij_to_k
-
-      module pure subroutine util_flatten_eucl_k_to_ij(n, k, i, j)
-         implicit none
-         integer(I4B), intent(in)  :: n !! Number of bodies
-         integer(I8B), intent(in)  :: k !! Index of the flattened matrix
-         integer(I4B), intent(out) :: i !! Index of the ith body
-         integer(I4B), intent(out) :: j !! Index of the jth body
-      end subroutine util_flatten_eucl_k_to_ij
-
-      module subroutine util_flatten_eucl_plpl(self, param)
-         implicit none
-         class(swiftest_pl),         intent(inout) :: self  !! Swiftest massive body object
-         class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters
-      end subroutine
-
-      module subroutine util_flatten_eucl_pltp(self, pl, param)
-         implicit none
-         class(swiftest_tp),         intent(inout) :: self  !! Swiftest test particle object
-         class(swiftest_pl),         intent(in)    :: pl    !! Swiftest massive body object
-         class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters
-      end subroutine
+      module subroutine encounter_check_all_flat_plpl(nplplm, k_plpl, x, v, renc, dt, lencounter, loc_lvdotr)
+         implicit none
+         integer(I8B),                 intent(in)  :: nplplm     !! Total number of plm-pl encounters to test
+         integer(I4B), dimension(:,:), intent(in)  :: k_plpl     !! List of all pl-pl encounters
+         real(DP),     dimension(:,:), intent(in)  :: x          !! Position vectors of massive bodies
+         real(DP),     dimension(:,:), intent(in)  :: v          !! Velocity vectors of massive bodies
+         real(DP),     dimension(:),   intent(in)  :: renc      !! Hill's radii of massive bodies
+         real(DP),                     intent(in)  :: dt         !! Step size
+         logical,      dimension(:),   intent(out) :: lencounter !! Logical array indicating which pair is in an encounter state
+         logical,      dimension(:),   intent(out) :: loc_lvdotr !! Logical array indicating the sign of v .dot. x for each encounter
+      end subroutine encounter_check_all_flat_plpl
+
+      module subroutine encounter_check_all_triangular_plpl(npl, nplm, x, v, renc, dt, lvdotr, index1, index2, nenc)
+         implicit none
+         integer(I4B),                            intent(in)  :: npl    !! Total number of massive bodies
+         integer(I4B),                            intent(in)  :: nplm   !! Number of fully interacting massive bodies
+         real(DP),     dimension(:,:),            intent(in)  :: x      !! 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
+         logical,      dimension(:), allocatable, intent(out) :: lvdotr !! Logical flag indicating the sign of v .dot. x
+         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
+         integer(I4B),                            intent(out) :: nenc   !! Total number of encounters
+      end subroutine encounter_check_all_triangular_plpl
+
+      module subroutine encounter_check_all_triangular_pltp(npl, ntp, xpl, vpl, xtp, vtp, renc, dt, lvdotr, index1, index2, nenc)
+         implicit none
+         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)  :: xpl    !! Position vectors of massive bodies
+         real(DP),     dimension(:,:),            intent(in)  :: vpl    !! Velocity vectors of massive bodies
+         real(DP),     dimension(:,:),            intent(in)  :: xtp    !! 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
+         logical,      dimension(:), allocatable, intent(out) :: lvdotr !! Logical flag indicating the sign of v .dot. x
+         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
+         integer(I4B),                            intent(out) :: nenc   !! Total number of encounters
+      end subroutine encounter_check_all_triangular_pltp
+
+      module pure subroutine encounter_check_one(xr, yr, zr, vxr, vyr, vzr, dt, renc, lencounter, lvdotr)
+         !$omp declare simd(encounter_check_one)
+         implicit none
+         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)  :: dt            !! Step size
+         real(DP), intent(in)  :: renc          !! Critical encounter distance
+         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
+      end subroutine encounter_check_one
 
       module pure subroutine gr_kick_getaccb_ns_body(self, system, param)
          implicit none
@@ -1394,12 +1421,36 @@ end subroutine util_fill_arr_logical
    end interface
 
    interface
-      module subroutine util_rescale_system(self, param, mscale, dscale, tscale)
+      module pure subroutine util_flatten_eucl_ij_to_k(n, i, j, k)
+         !$omp declare simd(util_flatten_eucl_ij_to_k)
          implicit none
-         class(swiftest_nbody_system), intent(inout) :: self   !! Swiftest nbody system object
-         class(swiftest_parameters),   intent(inout) :: param  !! Current run configuration parameters. Returns with new values of the scale vactors and GU
-         real(DP),                     intent(in)    :: mscale, dscale, tscale !! Scale factors for mass, distance, and time units, respectively. 
-      end subroutine util_rescale_system
+         integer(I4B), intent(in)  :: n !! Number of bodies
+         integer(I4B), intent(in)  :: i !! Index of the ith body
+         integer(I4B), intent(in)  :: j !! Index of the jth body
+         integer(I8B), intent(out) :: k !! Index of the flattened matrix
+      end subroutine util_flatten_eucl_ij_to_k
+
+      module pure subroutine util_flatten_eucl_k_to_ij(n, k, i, j)
+         implicit none
+         integer(I4B), intent(in)  :: n !! Number of bodies
+         integer(I8B), intent(in)  :: k !! Index of the flattened matrix
+         integer(I4B), intent(out) :: i !! Index of the ith body
+         integer(I4B), intent(out) :: j !! Index of the jth body
+      end subroutine util_flatten_eucl_k_to_ij
+
+      module subroutine util_flatten_eucl_plpl(self, param)
+         implicit none
+         class(swiftest_pl),         intent(inout) :: self  !! Swiftest massive body object
+         class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters
+      end subroutine
+
+      module subroutine util_flatten_eucl_pltp(self, pl, param)
+         implicit none
+         class(swiftest_tp),         intent(inout) :: self  !! Swiftest test particle object
+         class(swiftest_pl),         intent(in)    :: pl    !! Swiftest massive body object
+         class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters
+      end subroutine
+
 
       module function util_minimize_bfgs(f, N, x0, eps, maxloop, lerr) result(x1)
          use lambda_function
@@ -1419,6 +1470,14 @@ module subroutine util_peri_tp(self, system, param)
          class(swiftest_nbody_system), intent(inout) :: system !! Swiftest nbody system object
          class(swiftest_parameters),   intent(in)    :: param  !! Current run configuration parameters
       end subroutine util_peri_tp
+
+
+      module subroutine util_rescale_system(self, param, mscale, dscale, tscale)
+         implicit none
+         class(swiftest_nbody_system), intent(inout) :: self   !! Swiftest nbody system object
+         class(swiftest_parameters),   intent(inout) :: param  !! Current run configuration parameters. Returns with new values of the scale vactors and GU
+         real(DP),                     intent(in)    :: mscale, dscale, tscale !! Scale factors for mass, distance, and time units, respectively. 
+      end subroutine util_rescale_system
    end interface
 
 
@@ -1543,6 +1602,18 @@ module subroutine util_set_rhill(self,cb)
          class(swiftest_cb), intent(inout) :: cb   !! Swiftest central body object
       end subroutine util_set_rhill
 
+      module subroutine util_set_renc_I4B(self, scale)
+         implicit none
+         class(swiftest_pl), intent(inout) :: self !! Swiftest massive body object
+         integer(I4B),       intent(in)    :: scale !! Input scale factor (multiplier of Hill's sphere size)
+      end subroutine util_set_renc_I4B
+
+      module subroutine util_set_renc_DP(self, scale)
+         implicit none
+         class(swiftest_pl), intent(inout) :: self !! Swiftest massive body object
+         real(DP),           intent(in)    :: scale !! Input scale factor (multiplier of Hill's sphere size)
+      end subroutine util_set_renc_DP
+
       module subroutine util_set_rhill_approximate(self,cb)
          implicit none
          class(swiftest_pl), intent(inout) :: self !! Swiftest massive body object
diff --git a/src/modules/symba_classes.f90 b/src/modules/symba_classes.f90
index a1c56a456..e6c01bacc 100644
--- a/src/modules/symba_classes.f90
+++ b/src/modules/symba_classes.f90
@@ -67,30 +67,31 @@ module symba_classes
       type(symba_kinship),       dimension(:), allocatable :: kin        !! Array of merger relationship structures that can account for multiple pairwise mergers in a single step
    contains
       procedure :: make_colliders  => symba_collision_make_colliders_pl !! When a single body is involved in more than one collision in a single step, it becomes part of a family
-      procedure :: flatten         => symba_util_flatten_eucl_plpl     !! Sets up the (i, j) -> k indexing used for the single-loop blocking Euclidean distance matrix
-      procedure :: discard         => symba_discard_pl               !! Process massive body discards
-      procedure :: drift           => symba_drift_pl                 !! Method for Danby drift in Democratic Heliocentric coordinates. Sets the mask to the current recursion level
-      procedure :: encounter_check => symba_encounter_check_pl       !! Checks if massive bodies are going through close encounters with each other
-      procedure :: accel_int       => symba_kick_getacch_int_pl      !! Compute direct cross (third) term heliocentric accelerations of massive bodiess, with no mutual interactions between bodies below GMTINY
-      procedure :: accel           => symba_kick_getacch_pl          !! Compute heliocentric accelerations of massive bodies
-      procedure :: setup           => symba_setup_pl                 !! Constructor method - Allocates space for the input number of bodies
-      procedure :: append          => symba_util_append_pl           !! Appends elements from one structure to another
-      procedure :: fill            => symba_util_fill_pl             !! "Fills" bodies from one object into another depending on the results of a mask (uses the UNPACK intrinsic)
-      procedure :: get_peri        => symba_util_peri_pl             !! Determine system pericenter passages for massive bodies
-      procedure :: rearray         => symba_util_rearray_pl          !! Clean up the massive body structures to remove discarded bodies and add new bodies
-      procedure :: reset_kinship   => symba_util_reset_kinship       !! Resets the kinship status of bodies
-      procedure :: resize          => symba_util_resize_pl           !! Checks the current size of a SyMBA massive body against the requested size and resizes it if it is too small.
-      procedure :: sort            => symba_util_sort_pl             !! Sorts body arrays by a sortable componen
-      procedure :: rearrange       => symba_util_sort_rearrange_pl   !! Rearranges the order of array elements of body based on an input index array. Used in sorting methods
-      procedure :: spill           => symba_util_spill_pl            !! "Spills" bodies from one object to another depending on the results of a mask (uses the PACK intrinsic)
+      procedure :: flatten         => symba_util_flatten_eucl_plpl      !! Sets up the (i, j) -> k indexing used for the single-loop blocking Euclidean distance matrix
+      procedure :: discard         => symba_discard_pl                  !! Process massive body discards
+      procedure :: drift           => symba_drift_pl                    !! Method for Danby drift in Democratic Heliocentric coordinates. Sets the mask to the current recursion level
+      procedure :: encounter_check => symba_encounter_check_pl          !! Checks if massive bodies are going through close encounters with each other
+      procedure :: accel_int       => symba_kick_getacch_int_pl         !! Compute direct cross (third) term heliocentric accelerations of massive bodiess, with no mutual interactions between bodies below GMTINY
+      procedure :: accel           => symba_kick_getacch_pl             !! Compute heliocentric accelerations of massive bodies
+      procedure :: setup           => symba_setup_pl                    !! Constructor method - Allocates space for the input number of bodies
+      procedure :: append          => symba_util_append_pl              !! Appends elements from one structure to another
+      procedure :: fill            => symba_util_fill_pl                !! "Fills" bodies from one object into another depending on the results of a mask (uses the UNPACK intrinsic)
+      procedure :: get_peri        => symba_util_peri_pl                !! Determine system pericenter passages for massive bodies
+      procedure :: rearray         => symba_util_rearray_pl             !! Clean up the massive body structures to remove discarded bodies and add new bodies
+      procedure :: reset_kinship   => symba_util_reset_kinship          !! Resets the kinship status of bodies
+      procedure :: resize          => symba_util_resize_pl              !! Checks the current size of a SyMBA massive body against the requested size and resizes it if it is too small.
+      procedure :: set_renc_I4B    => symba_util_set_renc               !! Sets the critical radius for encounter given an input recursion depth
+      procedure :: sort            => symba_util_sort_pl                !! Sorts body arrays by a sortable componen
+      procedure :: rearrange       => symba_util_sort_rearrange_pl      !! Rearranges the order of array elements of body based on an input index array. Used in sorting methods
+      procedure :: spill           => symba_util_spill_pl               !! "Spills" bodies from one object to another depending on the results of a mask (uses the PACK intrinsic)
    end type symba_pl
 
    type, extends(symba_pl) :: symba_merger
       integer(I4B), dimension(:), allocatable :: ncomp
    contains
-      procedure :: append          => symba_util_append_merger       !! Appends elements from one structure to another
-      procedure :: resize          => symba_util_resize_merger       !! Checks the current size of a SyMBA merger list against the requested size and resizes it if it is too small.
-      procedure :: setup           => symba_setup_merger             !! Constructor method - Allocates space for the input number of bodies
+      procedure :: append          => symba_util_append_merger  !! Appends elements from one structure to another
+      procedure :: resize          => symba_util_resize_merger  !! Checks the current size of a SyMBA merger list against the requested size and resizes it if it is too small.
+      procedure :: setup           => symba_setup_merger        !! Constructor method - Allocates space for the input number of bodies
    end type symba_merger
 
    !********************************************************************************************************************************
@@ -323,6 +324,12 @@ module subroutine symba_util_flatten_eucl_plpl(self, param)
          class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters
       end subroutine symba_util_flatten_eucl_plpl
 
+      module subroutine symba_util_set_renc(self, scale)
+         implicit none
+         class(symba_pl), intent(inout) :: self !! SyMBA massive body object
+         integer(I4B),    intent(in)    :: scale !! Current recursion depth
+      end subroutine symba_util_set_renc
+
       module subroutine symba_io_param_reader(self, unit, iotype, v_list, iostat, iomsg) 
          implicit none
          class(symba_parameters), intent(inout) :: self       !! Current run configuration parameters with SyMBA additionss
diff --git a/src/rmvs/rmvs_encounter_check.f90 b/src/rmvs/rmvs_encounter_check.f90
index dd492032b..9a56259ef 100644
--- a/src/rmvs/rmvs_encounter_check.f90
+++ b/src/rmvs/rmvs_encounter_check.f90
@@ -18,10 +18,12 @@ module function rmvs_encounter_check_tp(self, param, system, dt) result(lencount
       ! Result
       logical                                 :: lencounter  !! Returns true if there is at least one close encounter
       ! Internals
-      integer(I4B)                            :: i, j
+      integer(I4B)                            :: i, j, nenc
       real(DP)                                :: xr, yr, zr, vxr, vyr, vzr
-      real(DP), dimension(system%pl%nbody)    :: r2crit
-      logical                                 :: lflag, lvdotr
+      real(DP), dimension(system%pl%nbody)    :: rcrit
+      logical                                 :: lflag
+      logical, dimension(:),      allocatable :: lvdotr
+      integer(I4B), dimension(:), allocatable :: index1, index2
 
       ! The minimization and linear solvers can sometimes lead to floating point exceptions. Rather than halting the code entirely if this occurs, we
       ! can simply fail the attempt and try again. So we need to turn off any floating point exception halting modes temporarily 
@@ -31,27 +33,17 @@ module function rmvs_encounter_check_tp(self, param, system, dt) result(lencount
 
       select type(pl => system%pl)
       class is (rmvs_pl)
-         associate(tp => self, ntp => self%nbody, npl => pl%nbody, rts => system%rts)
-            r2crit(1:npl) = (rts * pl%rhill(1:npl))**2
+         associate(tp => self, ntp => self%nbody, npl => pl%nbody)
             tp%plencP(1:ntp) = 0
-            !$omp parallel do default(private)&
-            !$omp shared(npl, ntp, tp, pl, dt, r2crit)
-            do j = 1, npl
-               do i = 1, ntp
-                  if ((.not.tp%lmask(i)).or.(tp%plencP(i) /= 0)) cycle
-                  xr = tp%xh(1, i) - pl%xbeg(1, j)
-                  yr = tp%xh(2, i) - pl%xbeg(2, j)
-                  zr = tp%xh(3, i) - pl%xbeg(3, j)
-                  vxr = tp%vh(1, i) - pl%vbeg(1, j)
-                  vyr = tp%vh(2, i) - pl%vbeg(2, j)
-                  vzr = tp%vh(3, i) - pl%vbeg(3, j)
-                  call rmvs_chk_ind(xr, yr, zr, vxr, vyr, vzr, dt, r2crit(j), lflag, lvdotr)
-                  if (lflag) tp%plencP(i) = j
+            call encounter_check_all_triangular_pltp(npl, ntp, pl%xbeg, pl%vbeg, tp%xh, tp%vh, pl%renc, dt, lvdotr, index1, index2, nenc)
+
+            lencounter = (nenc > 0)
+            if (lencounter) then
+               tp%plencP(index2(1:nenc)) = index1(1:nenc)
+               do j = 1, npl
+                  pl%nenc(j) = count(tp%plencP(1:ntp) == j)
                end do
-               pl%nenc(j) = count(tp%plencP(1:ntp) == j)
-            end do
-            !$omp end parallel do
-            lencounter = any(pl%nenc(1:npl) > 0)
+            end if
          end associate
       end select
 
@@ -59,43 +51,4 @@ module function rmvs_encounter_check_tp(self, param, system, dt) result(lencount
    end function rmvs_encounter_check_tp
 
 
-   module pure subroutine rmvs_chk_ind(xr, yr, zr, vxr, vyr, vzr, dt, r2crit, lencounter, lvdotr)
-      !$omp declare simd(rmvs_chk_ind)
-      !! 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)  :: dt            !! Step size
-      real(DP), intent(in)  :: r2crit        !! Square of the critical encounter distance
-      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) :: r2min, r2, v2, vdotr
-
-      r2 = xr**2 + yr**2 + zr**2
-      lencounter = (r2 < r2crit) 
-      if (lencounter) return 
-
-      vdotr = vxr * xr + vyr * yr + vzr * zr
-      lvdotr = (vdotr < 0.0_DP)
-      if (.not.lvdotr) return
-     
-      v2 = vxr**2 + vyr**2 + vzr**2
-
-      if (-vdotr < v2 * dt) then
-         r2min = r2 - vdotr**2 / v2
-      else
-         r2min = r2 + 2 * vdotr * dt + v2 * dt**2
-      end if
-      lencounter = (r2min <= r2crit)
-
-      return
-   end subroutine rmvs_chk_ind
-
 end submodule s_rmvs_chk
diff --git a/src/rmvs/rmvs_step.f90 b/src/rmvs/rmvs_step.f90
index 005b6f68e..673188f0e 100644
--- a/src/rmvs/rmvs_step.f90
+++ b/src/rmvs/rmvs_step.f90
@@ -34,7 +34,7 @@ module subroutine rmvs_step_system(self, param, t, dt)
                   allocate(vbeg, source=pl%vh)
                   call pl%set_beg_end(xbeg = xbeg, vbeg = vbeg)
                   ! ****** Check for close encounters ***** !
-                  system%rts = RHSCALE
+                  call pl%set_renc(RHSCALE)
                   lencounter = tp%encounter_check(param, system, dt)
                   if (lencounter) then
                      lfirstpl = pl%lfirst
@@ -172,12 +172,12 @@ subroutine rmvs_step_out(cb, pl, tp, system, param, t, dt)
          elsewhere  
             tp%lperi(1:ntp) = .false.
          end where
+         call pl%set_renc(RHPSCALE)
          do outer_index = 1, NTENC
             outer_time = t + (outer_index - 1) * dto
             call pl%set_beg_end(xbeg = pl%outer(outer_index - 1)%x(:, 1:npl), &
                                 vbeg = pl%outer(outer_index - 1)%v(:, 1:npl), &
                                 xend = pl%outer(outer_index    )%x(:, 1:npl))
-            system%rts = RHPSCALE
             lencounter = tp%encounter_check(param, system, dto) 
             if (lencounter) then
                ! Interpolate planets in inner encounter region
diff --git a/src/setup/setup.f90 b/src/setup/setup.f90
index 2dc57d952..7be3eb2c4 100644
--- a/src/setup/setup.f90
+++ b/src/setup/setup.f90
@@ -310,6 +310,7 @@ module subroutine setup_pl(self, n, param)
       if (allocated(self%mass)) deallocate(self%mass)
       if (allocated(self%Gmass)) deallocate(self%Gmass)
       if (allocated(self%rhill)) deallocate(self%rhill)
+      if (allocated(self%renc)) deallocate(self%renc)
       if (allocated(self%radius)) deallocate(self%radius)
       if (allocated(self%density)) deallocate(self%density)
       if (allocated(self%rot)) deallocate(self%rot)
@@ -323,10 +324,12 @@ module subroutine setup_pl(self, n, param)
       allocate(self%mass(n))
       allocate(self%Gmass(n))
       allocate(self%rhill(n))
+      allocate(self%renc(n))
 
       self%mass(:) = 0.0_DP
       self%Gmass(:) = 0.0_DP
       self%rhill(:) = 0.0_DP
+      self%renc(:) = 0.0_DP
 
       self%nplpl = 0   
 
diff --git a/src/symba/symba_encounter_check.f90 b/src/symba/symba_encounter_check.f90
index 10ea00144..91525b821 100644
--- a/src/symba/symba_encounter_check.f90
+++ b/src/symba/symba_encounter_check.f90
@@ -2,126 +2,6 @@
    use swiftest
 contains
 
-   subroutine symba_encounter_check_all_flat(nplplm, k_plpl, x, v, rhill,  dt, irec, lencounter, loc_lvdotr)
-      !! author: David A. Minton
-      !!
-      !! Check for encounters between massive bodies. Split off from the main subroutine for performance.
-      !! This is the flat (single loop) version.
-      implicit none
-      integer(I8B),                 intent(in)  :: nplplm     !! Total number of plm-pl interactions
-      integer(I4B), dimension(:,:), intent(in)  :: k_plpl     !! List of all pl-pl interactions
-      real(DP),     dimension(:,:), intent(in)  :: x          !! Position vectors of massive bodies
-      real(DP),     dimension(:,:), intent(in)  :: v          !! Velocity vectors of massive bodies
-      real(DP),     dimension(:),   intent(in)  :: rhill      !! Hill's radii of massive bodies
-      real(DP),                     intent(in)  :: dt         !! Step size
-      integer(I4B),                 intent(in)  :: irec       !! Current recursion depth
-      logical,      dimension(:),   intent(out) :: lencounter !! Logical array indicating which pair is in an encounter state
-      logical,      dimension(:),   intent(out) :: loc_lvdotr !! Logical array indicating the sign of v .dot. x for each encounter
-      ! Internals
-      integer(I8B) :: k
-      integer(I4B) :: i, j
-      real(DP) :: xr, yr, zr, vxr, vyr, vzr, rhill1, rhill2
-      
-      !$omp parallel do simd default(private) schedule(static)&
-      !$omp shared(nplplm, k_plpl, x, v, rhill,  dt, irec, lencounter, loc_lvdotr) &
-      !$omp lastprivate(xr, yr, zr, vxr, vyr, vzr, rhill1, rhill2)
-      do k = 1_I8B, nplplm
-         i = k_plpl(1, k)
-         j = k_plpl(2, k)
-         xr = x(1, j) - x(1, i)
-         yr = x(2, j) - x(2, i)
-         zr = x(3, j) - x(3, i)
-         vxr = v(1, j) - v(1, i)
-         vyr = v(2, j) - v(2, i)
-         vzr = v(3, j) - v(3, i)
-         rhill1 = rhill(i)
-         rhill2 = rhill(j)
-         call symba_encounter_check_one(xr, yr, zr, vxr, vyr, vzr, rhill1, rhill2, dt, irec, lencounter(k), loc_lvdotr(k))
-      end do
-      !$omp end parallel do simd
-
-      return
-   end subroutine symba_encounter_check_all_flat
-
-
-   subroutine symba_encounter_check_all_triangular(npl, nplm, x, v, rhill,  dt, irec, lvdotr, index1, index2, nenc)
-      !! 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
-      integer(I4B),                            intent(in)  :: npl    !! Total number of massive bodies
-      integer(I4B),                            intent(in)  :: nplm   !! Number of fully interacting massive bodies
-      real(DP),     dimension(:,:),            intent(in)  :: x      !! Position vectors of massive bodies
-      real(DP),     dimension(:,:),            intent(in)  :: v      !! Velocity vectors of massive bodies
-      real(DP),     dimension(:),              intent(in)  :: rhill  !! Hill's radii of massive bodies
-      real(DP),                                intent(in)  :: dt     !! Step size
-      integer(I4B),                            intent(in)  :: irec   !! Current recursion depth
-      logical,      dimension(:), allocatable, intent(out) :: lvdotr !! Logical flag indicating the sign of v .dot. x
-      integer(I4B), dimension(:), allocatable, intent(out) :: index1 !! List of indices for body 1 in each interaction
-      integer(I4B), dimension(:), allocatable, intent(out) :: index2 !! List of indices for body 2 in each interaction
-      integer(I4B),                            intent(out) :: nenc   !! Total number of interactions
-      ! Internals
-      integer(I4B) :: i, j, nenci, j0, j1
-      real(DP) :: xr, yr, zr, vxr, vyr, vzr, rhill1, rhill2
-      logical, dimension(npl) :: lencounteri, lvdotri
-      integer(I4B), dimension(npl) :: ind_arr
-      type lenctype
-         logical, dimension(:), allocatable :: lvdotr
-         integer(I4B), dimension(:), allocatable :: index2
-         integer(I4B) :: nenc
-      end type
-      type(lenctype), dimension(nplm) :: lenc
-     
-      ind_arr(:) = [(i, i = 1, npl)]
-      !$omp parallel do default(private) schedule(static)&
-      !$omp shared(npl, nplm, x, v, rhill,  dt, irec, lenc, ind_arr) &
-      !$omp lastprivate(xr, yr, zr, vxr, vyr, vzr, rhill1, rhill2, lencounteri, lvdotri)
-      do i = 1, nplm
-         do concurrent(j = i+1:npl)
-            xr = x(1, j) - x(1, i)
-            yr = x(2, j) - x(2, i)
-            zr = x(3, j) - x(3, i)
-            vxr = v(1, j) - v(1, i)
-            vyr = v(2, j) - v(2, i)
-            vzr = v(3, j) - v(3, i)
-            rhill1 = rhill(i)
-            rhill2 = rhill(j)
-            call symba_encounter_check_one(xr, yr, zr, vxr, vyr, vzr, rhill1, rhill2, dt, irec, lencounteri(j), lvdotri(j))
-         end do
-         nenci = count(lencounteri(i+1:npl))
-         if (nenci > 0) then
-            allocate(lenc(i)%lvdotr(nenci), lenc(i)%index2(nenci))
-            lenc(i)%nenc = nenci
-            lenc(i)%lvdotr(:) = pack(lvdotri(i+1:npl), lencounteri(i+1:npl)) 
-            lenc(i)%index2(:) = pack(ind_arr(i+1:npl), lencounteri(i+1:npl)) 
-         end if
-      end do
-      !$omp end parallel do
-
-      associate(nenc_arr => lenc(:)%nenc)
-         nenc = sum(nenc_arr(1:nplm))
-      end associate
-      if (nenc > 0) then
-         allocate(lvdotr(nenc))
-         allocate(index1(nenc))
-         allocate(index2(nenc))
-         j0 = 1
-         do i = 1, nplm
-            if (lenc(i)%nenc > 0) then
-               j1 = j0 + lenc(i)%nenc - 1
-               lvdotr(j0:j1) = lenc(i)%lvdotr(:)
-               index1(j0:j1) = i
-               index2(j0:j1) = lenc(i)%index2(:)
-               j0 = j1 + 1
-            end if
-         end do
-      end if
-
-      return
-   end subroutine symba_encounter_check_all_triangular
-
-
    module function symba_encounter_check_pl(self, param, system, dt, irec) result(lany_encounter)
       !! author: David A. Minton
       !!
@@ -167,7 +47,7 @@ module function symba_encounter_check_pl(self, param, system, dt, irec) result(l
             allocate(lencounter(nplplm))
             allocate(loc_lvdotr(nplplm))
   
-            call symba_encounter_check_all_flat(nplplm, pl%k_plpl, pl%xh, pl%vh, pl%rhill, dt, irec, lencounter, loc_lvdotr)
+            call encounter_check_all_flat_plpl(nplplm, pl%k_plpl, pl%xh, pl%vh, pl%renc, dt, lencounter, loc_lvdotr)
 
             nenc = count(lencounter(:))
 
@@ -187,7 +67,7 @@ module function symba_encounter_check_pl(self, param, system, dt, irec) result(l
             end if					
          else
             nplm = pl%nplm
-            call symba_encounter_check_all_triangular(npl, nplm, pl%xh, pl%vh, pl%rhill, dt, irec, lvdotr, index1, index2, nenc)
+            call encounter_check_all_triangular_plpl(npl, nplm, pl%xh, pl%vh, pl%renc, dt, lvdotr, index1, index2, nenc)
             lany_encounter = nenc > 0
             if (lany_encounter) then
                call plplenc_list%resize(nenc)
@@ -245,7 +125,7 @@ module function symba_encounter_check(self, param, system, dt, irec) result(lany
       integer(I4B)              :: i, j, k, lidx, nenc_enc
       real(DP), dimension(NDIM) :: xr, vr
       logical                   :: isplpl
-      real(DP)                  :: rlim2, rji2
+      real(DP)                  :: rlim2, rji2, rcrit12
       logical, dimension(:), allocatable :: lencmask, lencounter
       integer(I4B), dimension(:), allocatable :: encidx
 
@@ -279,7 +159,8 @@ module function symba_encounter_check(self, param, system, dt, irec) result(lany
                   j = self%index2(k)
                   xr(:) = pl%xh(:,j) - pl%xh(:,i)
                   vr(:) = pl%vb(:,j) - pl%vb(:,i)
-                  call symba_encounter_check_one(xr(1), xr(2), xr(3), vr(1), vr(2), vr(3), pl%rhill(i), pl%rhill(j), dt, irec, lencounter(lidx), self%lvdotr(k))
+                  rcrit12 = pl%renc(i) + pl%renc(j)
+                  call encounter_check_one(xr(1), xr(2), xr(3), vr(1), vr(2), vr(3), rcrit12, dt, lencounter(lidx), self%lvdotr(k))
                   if (lencounter(lidx)) then
                      rlim2 = (pl%radius(i) + pl%radius(j))**2
                      rji2 = dot_product(xr(:), xr(:))! Check to see if these are physically overlapping bodies first, which we should ignore
@@ -293,7 +174,7 @@ module function symba_encounter_check(self, param, system, dt, irec) result(lany
                   j = self%index2(k)
                   xr(:) = tp%xh(:,j) - pl%xh(:,i)
                   vr(:) = tp%vb(:,j) - pl%vb(:,i)
-                  call symba_encounter_check_one(xr(1), xr(2), xr(3), vr(1), vr(2), vr(3), pl%rhill(i), 0.0_DP, dt, irec, lencounter(lidx), self%lvdotr(k))
+                  call encounter_check_one(xr(1), xr(2), xr(3), vr(1), vr(2), vr(3), pl%rhill(i), dt, lencounter(lidx), self%lvdotr(k))
                   if (lencounter(lidx)) then
                      rlim2 = (pl%radius(i))**2
                      rji2 = dot_product(xr(:), xr(:))! Check to see if these are physically overlapping bodies first, which we should ignore
@@ -346,36 +227,24 @@ module function symba_encounter_check_tp(self, param, system, dt, irec) result(l
       real(DP)                                  :: r2crit, vdotr, r2, v2, tmin, r2min, term2
       integer(I4B)                              :: i, j, k,nenc, plind, tpind
       real(DP),     dimension(NDIM)             :: xr, vr
-      logical,      dimension(:,:), allocatable :: lencounter, loc_lvdotr
+      real(DP)                                  :: rshell_irec
+      logical,      dimension(:),   allocatable :: lvdotr
+      integer(I4B), dimension(:),   allocatable :: index1, index2
   
       if (self%nbody == 0) return
 
       associate(tp => self, ntp => self%nbody, pl => system%pl, npl => system%pl%nbody)
-         allocate(lencounter(ntp, npl), loc_lvdotr(ntp, npl))
-         lencounter(:,:) = .false.
+         call encounter_check_all_triangular_pltp(npl, ntp, pl%xh, pl%vh, tp%xh, tp%vh, pl%renc, dt, lvdotr, index1, index2, nenc) 
    
-         do j = 1, npl
-            do i = 1, ntp
-               xr(1) = tp%xh(1, i) - pl%xh(1, j)
-               xr(2) = tp%xh(2, i) - pl%xh(2, j)
-               xr(3) = tp%xh(3, i) - pl%xh(3, j)
-               vr(1) = tp%vh(1, i) - pl%vh(1, j)
-               vr(2) = tp%vh(2, i) - pl%vh(2, j)
-               vr(3) = tp%vh(3, i) - pl%vh(3, j)
-               call symba_encounter_check_one(xr(1), xr(2), xr(3), vr(1), vr(2), vr(3), pl%rhill(j), 0.0_DP, dt, irec, lencounter(i,j), loc_lvdotr(i,j))
-            end do
-         end do
-
-         nenc = count(lencounter(:,:))
          lany_encounter = nenc > 0
          if (lany_encounter) then 
             associate(pltpenc_list => system%pltpenc_list)
                call pltpenc_list%resize(nenc)
                pltpenc_list%status(1:nenc) = ACTIVE
                pltpenc_list%level(1:nenc) = irec
-               pltpenc_list%lvdotr(1:nenc) = pack(loc_lvdotr(1:ntp, 1:npl), lencounter(1:ntp, 1:npl))
-               pltpenc_list%index1(1:nenc) = pack(spread([(i, i = 1, npl)], dim=1, ncopies=ntp), lencounter(1:ntp, 1:npl)) 
-               pltpenc_list%index2(1:nenc) = pack(spread([(i, i = 1, ntp)], dim=2, ncopies=npl), lencounter(1:ntp, 1:npl))
+               call move_alloc(index1, pltpenc_list%index1)
+               call move_alloc(index2, pltpenc_list%index2)
+               call move_alloc(lvdotr, pltpenc_list%lvdotr)
                pltpenc_list%id1(1:nenc) = pl%id(pltpenc_list%index1(1:nenc))
                pltpenc_list%id2(1:nenc) = tp%id(pltpenc_list%index2(1:nenc))
                select type(pl)
@@ -400,34 +269,4 @@ module function symba_encounter_check_tp(self, param, system, dt, irec) result(l
       return
    end function symba_encounter_check_tp
 
-
-   pure subroutine symba_encounter_check_one(xr, yr, zr, vxr, vyr, vzr, rhill1, rhill2, dt, irec, lencounter, lvdotr)
-      !$omp declare simd(symba_encounter_check_one)
-      !! author: David A. Minton
-      !!
-      !! Check for an encounter.
-      !!
-      !! Adapted from David E. Kaufmann's Swifter routine: symba_chk.f90
-      !! Adapted from Hal Levison's Swift routine symba5_chk.f
-      implicit none
-      ! Arguments
-      real(DP),     intent(in)  :: xr, yr, zr, vxr, vyr, vzr
-      real(DP),     intent(in)  :: rhill1, rhill2, dt
-      integer(I4B), intent(in)  :: irec
-      logical,      intent(out) :: lencounter, lvdotr
-      ! Internals
-      real(DP)     :: r2crit, rshell_irec
-      integer(I4B) :: i
-
-      rshell_irec = 1._DP
-      do i = 1, irec
-         rshell_irec = rshell_irec * RSHELL
-      end do
-      r2crit = (rhill1 + rhill2) * RHSCALE * rshell_irec
-      r2crit = r2crit**2
-      call rmvs_chk_ind(xr, yr, zr, vxr, vyr, vzr, dt, r2crit, lencounter, lvdotr)
-
-      return
-   end subroutine symba_encounter_check_one
-
 end submodule s_symba_encounter_check
diff --git a/src/symba/symba_step.f90 b/src/symba/symba_step.f90
index dc07bf561..70317cd70 100644
--- a/src/symba/symba_step.f90
+++ b/src/symba/symba_step.f90
@@ -172,6 +172,7 @@ module recursive subroutine symba_step_recur_system(self, param, t, ireci)
                else
                   nloops = NTENC
                end if
+               call system%pl%set_renc(irecp)
                do j = 1, nloops
                   lencounter = plplenc_list%encounter_check(param, system, dtl, irecp) .or. pltpenc_list%encounter_check(param, system, dtl, irecp)
                    
@@ -250,6 +251,7 @@ module subroutine symba_step_reset_system(self, param)
                      call system%plplenc_list%setup(nenc_old)
                      system%plplenc_list%nenc = 0
                      call system%plplcollision_list%setup(0)
+                     call system%pl%set_renc(0)
                   end if
             
                   if (ntp > 0) then
diff --git a/src/symba/symba_util.f90 b/src/symba/symba_util.f90
index 86250272e..df33f7568 100644
--- a/src/symba/symba_util.f90
+++ b/src/symba/symba_util.f90
@@ -668,6 +668,31 @@ module subroutine symba_util_resize_tp(self, nnew)
       return
    end subroutine symba_util_resize_tp
 
+   
+   module subroutine symba_util_set_renc(self, scale)
+      !! author: David A. Minton
+      !!
+      !! Sets the critical radius for encounter given an input recursion depth
+      !!
+      implicit none
+      ! Arguments
+      class(symba_pl), intent(inout) :: self !! SyMBA massive body object
+      integer(I4B),    intent(in)    :: scale !! Current recursion depth
+      ! Internals
+      integer(I4B) :: i
+      real(DP)     :: rshell_irec
+
+      associate(pl => self, npl => self%nbody)
+         rshell_irec = 1._DP
+         do i = 1, scale
+            rshell_irec = rshell_irec * RSHELL
+         end do
+         pl%renc(1:npl) = pl%rhill(1:npl) * RHSCALE * rshell_irec
+      end associate
+
+      return
+   end subroutine symba_util_set_renc
+
 
    module subroutine symba_util_sort_pl(self, sortby, ascending)
       !! author: David A. Minton
diff --git a/src/util/util_append.f90 b/src/util/util_append.f90
index ffe576c85..c518b3df3 100644
--- a/src/util/util_append.f90
+++ b/src/util/util_append.f90
@@ -271,6 +271,7 @@ module subroutine util_append_pl(self, source, lsource_mask)
             call util_append(self%mass, source%mass, nold, nsrc, lsource_mask)
             call util_append(self%Gmass, source%Gmass, nold, nsrc, lsource_mask)
             call util_append(self%rhill, source%rhill, nold, nsrc, lsource_mask)
+            call util_append(self%renc, source%renc, nold, nsrc, lsource_mask)
             call util_append(self%radius, source%radius, nold, nsrc, lsource_mask)
             call util_append(self%xbeg, source%xbeg, nold, nsrc, lsource_mask)
             call util_append(self%xend, source%xend, nold, nsrc, lsource_mask)
diff --git a/src/util/util_fill.f90 b/src/util/util_fill.f90
index 26b21310e..bd0014be4 100644
--- a/src/util/util_fill.f90
+++ b/src/util/util_fill.f90
@@ -197,6 +197,7 @@ module subroutine util_fill_pl(self, inserts, lfill_list)
             call util_fill(keeps%mass, inserts%mass, lfill_list)
             call util_fill(keeps%Gmass, inserts%Gmass, lfill_list)
             call util_fill(keeps%rhill, inserts%rhill, lfill_list)
+            call util_fill(keeps%renc, inserts%renc, lfill_list)
             call util_fill(keeps%radius, inserts%radius, lfill_list)
             call util_fill(keeps%density, inserts%density, lfill_list)
             call util_fill(keeps%k2, inserts%k2, lfill_list)
diff --git a/src/util/util_resize.f90 b/src/util/util_resize.f90
index dfef8771b..18b874607 100644
--- a/src/util/util_resize.f90
+++ b/src/util/util_resize.f90
@@ -361,6 +361,7 @@ module subroutine util_resize_pl(self, nnew)
       call util_resize(self%mass, nnew)
       call util_resize(self%Gmass, nnew)
       call util_resize(self%rhill, nnew)
+      call util_resize(self%renc, nnew)
       call util_resize(self%radius, nnew)
       call util_resize(self%xbeg, nnew)
       call util_resize(self%xend, nnew)
diff --git a/src/util/util_set.f90 b/src/util/util_set.f90
index 80ed89b1f..b056d692c 100644
--- a/src/util/util_set.f90
+++ b/src/util/util_set.f90
@@ -161,6 +161,42 @@ module subroutine util_set_particle_info(self, name, particle_type, status, orig
    end subroutine util_set_particle_info
 
 
+   module subroutine util_set_renc_I4B(self, scale)
+      !! author: David A. Minton
+      !!
+      !! Sets the critical radius for encounter given an input scale factor
+      !!
+      implicit none
+      ! Arguments
+      class(swiftest_pl), intent(inout) :: self !! Swiftest massive body object
+      integer(I4B),       intent(in)    :: scale !! Input scale factor (multiplier of Hill's sphere size)
+
+      associate(pl => self, npl => self%nbody)
+         pl%renc(1:npl) = pl%rhill(1:npl) * scale
+      end associate
+
+      return
+   end subroutine util_set_renc_I4B
+
+
+   module subroutine util_set_renc_DP(self, scale)
+      !! author: David A. Minton
+      !!
+      !! Sets the critical radius for encounter given an input scale factor
+      !!
+      implicit none
+      ! Arguments
+      class(swiftest_pl), intent(inout) :: self !! Swiftest massive body object
+      real(DP),           intent(in)    :: scale !! Input scale factor (multiplier of Hill's sphere size)
+
+      associate(pl => self, npl => self%nbody)
+         pl%renc(1:npl) = pl%rhill(1:npl) * scale
+      end associate
+
+      return
+   end subroutine util_set_renc_DP
+
+
    module subroutine util_set_rhill(self,cb)
       !! author: David A. Minton
       !!
diff --git a/src/util/util_sort.f90 b/src/util/util_sort.f90
index b69627b0b..0b47de783 100644
--- a/src/util/util_sort.f90
+++ b/src/util/util_sort.f90
@@ -255,6 +255,8 @@ module subroutine util_sort_pl(self, sortby, ascending)
             call util_sort(direction * pl%Gmass(1:npl), ind(1:npl))
          case("rhill")
             call util_sort(direction * pl%rhill(1:npl), ind(1:npl))
+         case("renc")
+            call util_sort(direction * pl%renc(1:npl), ind(1:npl))
          case("radius")
             call util_sort(direction * pl%radius(1:npl), ind(1:npl))
          case("density")
diff --git a/src/util/util_spill.f90 b/src/util/util_spill.f90
index e9d243cb6..fca48d1c1 100644
--- a/src/util/util_spill.f90
+++ b/src/util/util_spill.f90
@@ -414,6 +414,7 @@ module subroutine util_spill_pl(self, discards, lspill_list, ldestructive)
             call util_spill(keeps%mass, discards%mass, lspill_list, ldestructive)
             call util_spill(keeps%Gmass, discards%Gmass, lspill_list, ldestructive)
             call util_spill(keeps%rhill, discards%rhill, lspill_list, ldestructive)
+            call util_spill(keeps%renc, discards%renc, lspill_list, ldestructive)
             call util_spill(keeps%radius, discards%radius, lspill_list, ldestructive)
             call util_spill(keeps%density, discards%density, lspill_list, ldestructive)
             call util_spill(keeps%k2, discards%k2, lspill_list, ldestructive)