diff --git a/src/modules/swiftest_classes.f90 b/src/modules/swiftest_classes.f90
index 8949afdfa..bebf2acfa 100644
--- a/src/modules/swiftest_classes.f90
+++ b/src/modules/swiftest_classes.f90
@@ -324,6 +324,7 @@ module swiftest_classes
       real(DP),     dimension(:),   allocatable :: t      !! Time of encounter
    contains
       procedure :: setup  => setup_encounter       !! A constructor that sets the number of encounters and allocates and initializes all arrays  
+      procedure :: append => util_append_encounter !! Appends elements from one structure to another
       procedure :: copy   => util_copy_encounter   !! Copies elements from the source encounter list into self.
       procedure :: spill  => util_spill_encounter  !! "Spills" bodies from one object to another depending on the results of a mask (uses the PACK intrinsic)
       procedure :: resize => util_resize_encounter !! Checks the current size of the encounter list against the required size and extends it by a factor of 2 more than requested if it is too small.
@@ -883,11 +884,18 @@ end subroutine util_append_arr_logical
    interface
       module subroutine util_append_body(self, source, lsource_mask)
          implicit none
-         class(swiftest_body),            intent(inout) :: self   !! Swiftest body object
-         class(swiftest_body),            intent(in)    :: source  !! Source object to append
+         class(swiftest_body),            intent(inout) :: self          !! Swiftest body object
+         class(swiftest_body),            intent(in)    :: source        !! Source object to append
          logical, dimension(:),           intent(in)    :: lsource_mask  !! Logical mask indicating which elements to append to
       end subroutine util_append_body
 
+      module subroutine util_append_encounter(self, source, lsource_mask)
+         implicit none
+         class(swiftest_encounter), intent(inout) :: self         !! Swiftest encounter list object
+         class(swiftest_encounter), intent(in)    :: source       !! Source object to append
+         logical, dimension(:),     intent(in)    :: lsource_mask !! Logical mask indicating which elements to append to
+      end subroutine util_append_encounter
+
       module subroutine util_append_pl(self, source, lsource_mask)
          implicit none
          class(swiftest_pl),              intent(inout) :: self         !! Swiftest massive body object
@@ -1198,6 +1206,43 @@ module subroutine util_sort_index_dp(arr,ind)
       end subroutine util_sort_index_dp
    end interface util_sort
 
+   interface util_sort_rearrange
+      module subroutine util_sort_rearrange_arr_char_string(arr, ind, n)
+         implicit none
+         character(len=STRMAX), dimension(:), allocatable, intent(inout) :: arr !! Destination array 
+         integer(I4B),          dimension(:),              intent(in)    :: ind !! Index to rearrange against
+         integer(I4B),                                     intent(in)    :: n   !! Number of elements in arr and ind to rearrange
+      end subroutine util_sort_rearrange_arr_char_string
+
+      module subroutine util_sort_rearrange_arr_DP(arr, ind, n)
+         implicit none
+         real(DP),     dimension(:), allocatable, intent(inout) :: arr !! Destination array 
+         integer(I4B), dimension(:),              intent(in)  :: ind !! Index to rearrange against
+         integer(I4B),                            intent(in)  :: n   !! Number of elements in arr and ind to rearrange
+      end subroutine util_sort_rearrange_arr_DP
+
+      module subroutine util_sort_rearrange_arr_DPvec(arr, ind, n)
+         implicit none
+         real(DP),     dimension(:,:), allocatable, intent(inout) :: arr !! Destination array 
+         integer(I4B), dimension(:),                intent(in)    :: ind !! Index to rearrange against
+         integer(I4B),                              intent(in)    :: n   !! Number of elements in arr and ind to rearrange
+      end subroutine util_sort_rearrange_arr_DPvec
+
+      module subroutine util_sort_rearrange_arr_I4B(arr, ind, n)
+         implicit none
+         integer(I4B), dimension(:), allocatable, intent(inout) :: arr !! Destination array 
+         integer(I4B), dimension(:),              intent(in)    :: ind !! Index to rearrange against
+         integer(I4B),                             intent(in)    :: n   !! Number of elements in arr and ind to rearrange
+      end subroutine util_sort_rearrange_arr_I4B
+
+      module subroutine util_sort_rearrange_arr_logical(arr, ind, n)
+         implicit none
+         logical,      dimension(:), allocatable, intent(inout) :: arr !! Destination array 
+         integer(I4B), dimension(:),              intent(in)    :: ind !! Index to rearrange against
+         integer(I4B),                            intent(in)    :: n   !! Number of elements in arr and ind to rearrange
+      end subroutine util_sort_rearrange_arr_logical
+   end interface util_sort_rearrange
+
    interface
       module subroutine util_sort_rearrange_body(self, ind)
          implicit none
diff --git a/src/modules/symba_classes.f90 b/src/modules/symba_classes.f90
index a07ce36fc..8516e3861 100644
--- a/src/modules/symba_classes.f90
+++ b/src/modules/symba_classes.f90
@@ -137,10 +137,11 @@ module symba_classes
       integer(I4B), dimension(:),   allocatable :: level  !! encounter recursion level
    contains
       procedure :: collision_check => symba_collision_check_encounter !! Checks if a test particle is going to collide with a massive body
-      procedure :: encounter_check => symba_encounter_check !! Checks if massive bodies are going through close encounters with each other
+      procedure :: encounter_check => symba_encounter_check           !! Checks if massive bodies are going through close encounters with each other
       procedure :: kick            => symba_kick_encounter            !! Kick barycentric velocities of active test particles within SyMBA recursion
       procedure :: setup           => symba_setup_encounter           !! A constructor that sets the number of encounters and allocates and initializes all arrays  
       procedure :: spill           => symba_util_spill_encounter      !! "Spills" bodies from one object to another depending on the results of a mask (uses the PACK intrinsic)
+      procedure :: append          => symba_util_append_encounter     !! Appends elements from one structure to another
    end type symba_encounter
 
    !********************************************************************************************************************************
@@ -149,6 +150,7 @@ module symba_classes
    !> SyMBA class for tracking pl-tp close encounters in a step
    type, extends(symba_encounter) :: symba_pltpenc
    contains
+      procedure :: resolve_collision => symba_collision_resolve_pltpenc !! Process the pl-tp collision list, then modifiy the massive bodies based on the outcome of the c
    end type symba_pltpenc
 
    !********************************************************************************************************************************
@@ -160,6 +162,7 @@ module symba_classes
       procedure :: extract_collisions     => symba_collision_encounter_extract_collisions !! Processes the pl-pl encounter list remove only those encounters that led to a collision
       procedure :: resolve_fragmentations => symba_collision_resolve_fragmentations       !! Process list of collisions, determine the collisional regime, and then create fragments
       procedure :: resolve_mergers        => symba_collision_resolve_mergers              !! Process list of collisions and merge colliding bodies together
+      procedure :: resolve_collision      => symba_collision_resolve_plplenc              !! Process the pl-pl collision list, then modifiy the massive bodies based on the outcome of the c
    end type symba_plplenc
 
    !********************************************************************************************************************************
@@ -221,6 +224,22 @@ module subroutine symba_collision_resolve_mergers(self, system, param)
          class(symba_parameters),   intent(in)    :: param  !! Current run configuration parameters with SyMBA additions
       end subroutine symba_collision_resolve_mergers
 
+      module subroutine symba_collision_resolve_plplenc(self, system, param, t)
+         implicit none
+         class(symba_plplenc),       intent(inout) :: self   !! SyMBA pl-pl encounter list
+         class(symba_nbody_system),  intent(inout) :: system !! SyMBA nbody system object
+         class(swiftest_parameters), intent(inout) :: param  !! Current run configuration parameters with SyMBA additions
+         real(DP),                   intent(in)    :: t      !! Current simulation time
+      end subroutine symba_collision_resolve_plplenc
+   
+      module subroutine symba_collision_resolve_pltpenc(self, system, param, t)
+         implicit none
+         class(symba_pltpenc),       intent(inout) :: self   !! SyMBA pl-tp encounter list
+         class(symba_nbody_system),  intent(inout) :: system !! SyMBA nbody system object
+         class(swiftest_parameters), intent(inout) :: param  !! Current run configuration parameters with SyMBA additions
+         real(DP),                   intent(in)    :: t      !! Current simulation time
+      end subroutine symba_collision_resolve_pltpenc
+
       module subroutine symba_discard_pl(self, system, param)
          use swiftest_classes, only : swiftest_nbody_system, swiftest_parameters
          implicit none
@@ -502,6 +521,13 @@ end subroutine symba_util_append_arr_kin
    end interface
 
    interface
+      module subroutine symba_util_append_encounter(self, source, lsource_mask)
+         implicit none
+         class(symba_encounter),    intent(inout) :: self         !! SyMBA encounter list object
+         class(swiftest_encounter), intent(in)    :: source       !! Source object to append
+         logical, dimension(:),     intent(in)    :: lsource_mask !! Logical mask indicating which elements to append to
+      end subroutine symba_util_append_encounter
+
       module subroutine symba_util_append_merger(self, source, lsource_mask)
          use swiftest_classes, only : swiftest_body
          implicit none
@@ -622,7 +648,25 @@ module subroutine symba_util_sort_tp(self, sortby, ascending)
          character(*),    intent(in)    :: sortby    !! Sorting attribute
          logical,         intent(in)    :: ascending !! Logical flag indicating whether or not the sorting should be in ascending or descending order
       end subroutine symba_util_sort_tp
+   end interface
+
+   interface util_sort_rearrange
+      module subroutine symba_util_sort_rearrange_arr_info(arr, ind, n)
+         implicit none
+         type(symba_particle_info),   dimension(:), allocatable, intent(inout) :: arr !! Destination array 
+         integer(I4B), dimension(:),              intent(in)    :: ind !! Index to rearrange against
+         integer(I4B),                            intent(in)    :: n   !! Number of elements in arr and ind to rearrange
+      end subroutine symba_util_sort_rearrange_arr_info
 
+      module subroutine symba_util_sort_rearrange_arr_kin(arr, ind, n)
+         implicit none
+         type(symba_kinship),         dimension(:), allocatable, intent(inout) :: arr !! Destination array 
+         integer(I4B), dimension(:),              intent(in)    :: ind !! Index to rearrange against
+         integer(I4B),                            intent(in)    :: n   !! Number of elements in arr and ind to rearrange
+      end subroutine symba_util_sort_rearrange_arr_kin
+   end interface util_sort_rearrange
+
+   interface
       module subroutine symba_util_sort_rearrange_pl(self, ind)
          implicit none
          class(symba_pl),               intent(inout) :: self !! SyMBA massive body object
diff --git a/src/symba/symba_collision.f90 b/src/symba/symba_collision.f90
index 2fb99ccee..0e6c69440 100644
--- a/src/symba/symba_collision.f90
+++ b/src/symba/symba_collision.f90
@@ -333,7 +333,7 @@ module function symba_collision_casesupercatastrophic(system, param, family, x,
       allocate(Ip_frag(NDIM, nfrag))
 
       mtot = sum(mass(:))
-      xcom(:) = (mass(1) * x(:,1) + mass(2) * x(:,2)) / mtot
+      xcom(:) = (mass(1) * x(:,1) + mass(2) * x(:,2)) / mtot 
       vcom(:) = (mass(1) * v(:,1) + mass(2) * v(:,2)) / mtot
 
       ! Get mass weighted mean of Ip and average density
@@ -393,7 +393,7 @@ module function symba_collision_check_encounter(self, system, param, t, dt, irec
       !! Adapted from Hal Levison's Swift routine symba5_merge.f
       implicit none
       ! Arguments
-      class(symba_encounter),       intent(inout) :: self           !! SyMBA pl-tp encounter list object
+      class(symba_encounter),     intent(inout) :: self           !! SyMBA pl-tp encounter list object
       class(symba_nbody_system),  intent(inout) :: system         !! SyMBA nbody system object
       class(swiftest_parameters), intent(in)    :: param          !! Current run configuration parameters 
       real(DP),                   intent(in)    :: t              !! current time
@@ -451,7 +451,6 @@ module function symba_collision_check_encounter(self, system, param, t, dt, irec
                   end do
                end if
 
-
                do k = 1, nenc
                   if (lcollision(k)) self%status(k) = COLLISION
                   self%t(k) = t
@@ -485,6 +484,14 @@ module function symba_collision_check_encounter(self, system, param, t, dt, irec
 
       lany_collision = any(lcollision(:))
 
+      ! Extract the pl-pl encounter list and return the plplcollision_list
+      if (lany_collision) then
+         select type(plplenc_list => self)
+         class is (symba_plplenc)
+            call plplenc_list%extract_collisions(system, param)
+         end select
+      end if
+
       return 
    end function symba_collision_check_encounter
 
@@ -604,7 +611,7 @@ function symba_collision_consolidate_familes(pl, param, idx_parent, family, x, v
 
       ! Find the barycenter of each body along with its children, if it has any
       do j = 1, 2
-         x(:, j)  = pl%xb(:, idx_parent(j))
+         x(:, j)  = pl%xh(:, idx_parent(j))
          v(:, j)  = pl%vb(:, idx_parent(j))
          ! Assume principal axis rotation about axis corresponding to highest moment of inertia (3rd Ip)
          if (param%lrotation) then
@@ -617,7 +624,7 @@ function symba_collision_consolidate_familes(pl, param, idx_parent, family, x, v
                idx_child = parent_child_index_array(j)%idx(i + 1)
                if (.not. pl%lcollision(idx_child)) cycle
                mchild = pl%mass(idx_child)
-               xchild(:) = pl%xb(:, idx_child)
+               xchild(:) = pl%xh(:, idx_child)
                vchild(:) = pl%vb(:, idx_child)
                volchild = (4.0_DP / 3.0_DP) * PI * pl%radius(idx_child)**3
                volume(j) = volume(j) + volchild
@@ -947,6 +954,10 @@ module subroutine symba_collision_resolve_fragmentations(self, system, param)
                if (.not. lgoodcollision) cycle
                if (any(pl%status(idx_parent(:)) /= COLLISION)) cycle ! One of these two bodies has already been resolved
 
+               ! Convert from DH to barycentric
+               x(:,1) = x(:,1) + cb%xb(:)
+               x(:,2) = x(:,2) + cb%xb(:)
+
                ! Convert all quantities to SI units and determine which of the pair is the projectile vs. target before sending them 
                ! to symba_regime
                if (mass(1) > mass(2)) then
@@ -1020,7 +1031,7 @@ module subroutine symba_collision_resolve_mergers(self, system, param)
       logical                                     :: lgoodcollision
       integer(I4B)                                :: i, status
 
-      associate(plpl_collisions => self, ncollisions => self%nenc, idx1 => self%index1, idx2 => self%index2)
+      associate(plpl_collisions => self, ncollisions => self%nenc, idx1 => self%index1, idx2 => self%index2, cb => system%cb)
          select type(pl => system%pl)
          class is (symba_pl)
             do i = 1, ncollisions
@@ -1029,6 +1040,11 @@ module subroutine symba_collision_resolve_mergers(self, system, param)
                lgoodcollision = symba_collision_consolidate_familes(pl, param, idx_parent, family, x, v, mass, radius, L_spin, Ip)
                if (.not. lgoodcollision) cycle
                if (any(pl%status(idx_parent(:)) /= COLLISION)) cycle ! One of these two bodies has already been resolved
+
+               ! Convert from DH to barycentric
+               x(:,1) = x(:,1) + cb%xb(:)
+               x(:,2) = x(:,2) + cb%xb(:)
+
                status = symba_collision_casemerge(system, param, family, x, v, mass, radius, L_spin, Ip) 
             end do
          end select
@@ -1037,4 +1053,77 @@ module subroutine symba_collision_resolve_mergers(self, system, param)
       return
    end subroutine symba_collision_resolve_mergers
 
+
+   module subroutine symba_collision_resolve_plplenc(self, system, param, t)
+      !! author: David A. Minton
+      !! 
+      !! Process the pl-pl collision list, then modifiy the massive bodies based on the outcome of the collision
+      !! 
+      implicit none
+      ! Arguments
+      class(symba_plplenc),       intent(inout) :: self   !! SyMBA pl-pl encounter list
+      class(symba_nbody_system),  intent(inout) :: system !! SyMBA nbody system object
+      class(swiftest_parameters), intent(inout) :: param  !! Current run configuration parameters with SyMBA additions
+      real(DP),                   intent(in)    :: t      !! Current simulation time
+      ! Internals
+      real(DP) :: Eorbit_before, Eorbit_after
+   
+      associate(plplenc_list => self, plplcollision_list => system%plplcollision_list)
+         select type(pl => system%pl)
+         class is (symba_pl)
+            select type(param)
+            class is (symba_parameters)
+               if (plplcollision_list%nenc == 0) return ! No collisions to resolve
+
+               write(*, *) "Collision between massive bodies detected at time t = ", t
+               if (param%lfragmentation) then
+                  call plplcollision_list%resolve_fragmentations(system, param)
+               else
+                  call plplcollision_list%resolve_mergers(system, param)
+               end if
+
+               ! Destroy the collision list now that the collisions are resolved
+               call plplcollision_list%setup(0)
+
+               ! Get the energy before the collision is resolved
+               if (param%lenergy) then
+                  call system%get_energy_and_momentum(param)
+                  Eorbit_before = system%te
+               end if
+
+               call pl%rearray(system, param)
+
+               if (param%lenergy) then
+                  call system%get_energy_and_momentum(param)
+                  Eorbit_after = system%te
+                  system%Ecollisions = system%Ecollisions + (Eorbit_after - Eorbit_before)
+               end if
+
+            end select 
+         end select
+      end associate
+
+      return
+   end subroutine symba_collision_resolve_plplenc
+
+
+   module subroutine symba_collision_resolve_pltpenc(self, system, param, t)
+      !! author: David A. Minton
+      !! 
+      !! Process the pl-tp collision list, then modifiy the massive bodies based on the outcome of the collision
+      !! 
+      implicit none
+      ! Arguments
+      class(symba_pltpenc),       intent(inout) :: self   !! SyMBA pl-pl encounter list
+      class(symba_nbody_system),  intent(inout) :: system !! SyMBA nbody system object
+      class(swiftest_parameters), intent(inout) :: param  !! Current run configuration parameters with SyMBA additions
+      real(DP),                   intent(in)    :: t      !! Current simulation tim
+  
+      call system%tp%discard(system, param)
+
+      return
+   end subroutine symba_collision_resolve_pltpenc
+
+
+
 end submodule s_symba_collision
\ No newline at end of file
diff --git a/src/symba/symba_discard.f90 b/src/symba/symba_discard.f90
index c35a1565c..70c0898a5 100644
--- a/src/symba/symba_discard.f90
+++ b/src/symba/symba_discard.f90
@@ -276,8 +276,7 @@ end subroutine symba_discard_peri_pl
    module subroutine symba_discard_pl(self, system, param)
       !! author: David A. Minton
       !!
-      !! Call the various flavors of discards for massive bodies in SyMBA runs, including discards due to colling with the central body, 
-      !! escaping the system, or colliding with each other.
+      !! Call the various flavors of discards for massive bodies in SyMBA runs, including discards due to colliding with the central body or escaping the system
       implicit none
       ! Arguments
       class(symba_pl),              intent(inout) :: self   !! SyMBA test particle object
@@ -291,38 +290,25 @@ module subroutine symba_discard_pl(self, system, param)
          select type(param)
          class is (symba_parameters)
             associate(pl => self, plplenc_list => system%plplenc_list, plplcollision_list => system%plplcollision_list)
-               call pl%h2b(system%cb) 
+               call plplenc_list%write(pl, pl, param)
 
-               ! First deal with the non pl-pl collisions
                call symba_discard_nonplpl(self, system, param)
 
-               ! Extract the pl-pl encounter list and return the plplcollision_list
-               call plplenc_list%extract_collisions(system, param)
-               call plplenc_list%write(pl, pl, param)
+               if (.not.any(pl%ldiscard(:))) return
 
-               if ((plplcollision_list%nenc > 0) .and. any(pl%lcollision(:))) then 
-                  write(*, *) "Collision between massive bodies detected at time t = ",param%t
-                  if (param%lfragmentation) then
-                     call plplcollision_list%resolve_fragmentations(system, param)
-                  else
-                     call plplcollision_list%resolve_mergers(system, param)
-                  end if
-                  ! Destroy the collision list now that the collisions are resolved
-                  call plplcollision_list%setup(0)
+               if (param%lenergy) then
+                  call system%get_energy_and_momentum(param)
+                  Eorbit_before = system%te
                end if
 
-               if (any(pl%ldiscard(:))) then
-                  if (param%lenergy) then
-                     call system%get_energy_and_momentum(param)
-                     Eorbit_before = system%te
-                  end if
-                  call symba_discard_nonplpl_conservation(self, system, param)
-                  call pl%rearray(system, param)
-                  if (param%lenergy) then
-                     call system%get_energy_and_momentum(param)
-                     Eorbit_after = system%te
-                     system%Ecollisions = system%Ecollisions + (Eorbit_after - Eorbit_before)
-                  end if
+               call symba_discard_nonplpl_conservation(self, system, param)
+
+               call pl%rearray(system, param)
+
+               if (param%lenergy) then
+                  call system%get_energy_and_momentum(param)
+                  Eorbit_after = system%te
+                  system%Ecollisions = system%Ecollisions + (Eorbit_after - Eorbit_before)
                end if
 
             end associate
diff --git a/src/symba/symba_step.f90 b/src/symba/symba_step.f90
index 69d7e2cbd..ca5aa43cf 100644
--- a/src/symba/symba_step.f90
+++ b/src/symba/symba_step.f90
@@ -108,32 +108,17 @@ module subroutine symba_step_set_recur_levels_system(self, ireci)
       ! Internals
       integer(I4B) :: k, irecp
 
-      associate(system => self, plplenc_list => self%plplenc_list, pltpenc_list => self%pltpenc_list)
+      associate(system => self, plplenc_list => self%plplenc_list, pltpenc_list => self%pltpenc_list, npl => self%pl%nbody, ntp => self%tp%nbody)
          select type(pl => self%pl)
          class is (symba_pl)
             select type(tp => self%tp)
             class is (symba_tp)
                irecp = ireci + 1
 
-               if (plplenc_list%nenc > 0) then
-                  do k = 1, plplenc_list%nenc
-                     associate(i => plplenc_list%index1(k), j => plplenc_list%index2(k))
-                        if (pl%levelg(i) == irecp) pl%levelg(i) = ireci
-                        if (pl%levelg(j) == irecp) pl%levelg(j) = ireci
-                     end associate
-                  end do
-                  where(plplenc_list%level(1:plplenc_list%nenc) == irecp) plplenc_list%level(1:plplenc_list%nenc) = ireci
-               end if
-
-               if (pltpenc_list%nenc > 0) then
-                  do k = 1, pltpenc_list%nenc
-                     associate(i => pltpenc_list%index1(k), j => pltpenc_list%index2(k))
-                        if (pl%levelg(i) == irecp) pl%levelg(i) = ireci
-                        if (tp%levelg(j) == irecp) tp%levelg(j) = ireci
-                     end associate
-                  end do
-                  where(pltpenc_list%level(1:pltpenc_list%nenc) == irecp) pltpenc_list%level(1:pltpenc_list%nenc) = ireci
-               end if
+               if (npl >0) where(pl%levelg(1:npl) == irecp) pl%levelg(1:npl) = ireci
+               if (ntp > 0) where(tp%levelg(1:ntp) == irecp) tp%levelg(1:ntp) = ireci
+               if (plplenc_list%nenc > 0) where(plplenc_list%level(1:plplenc_list%nenc) == irecp) plplenc_list%level(1:plplenc_list%nenc) = ireci
+               if (pltpenc_list%nenc > 0) where(pltpenc_list%level(1:pltpenc_list%nenc) == irecp) pltpenc_list%level(1:pltpenc_list%nenc) = ireci
 
                system%irec = ireci
 
@@ -212,8 +197,8 @@ module recursive subroutine symba_step_recur_system(self, param, t, ireci)
                      lplpl_collision = plplenc_list%collision_check(system, param, t+dtl, dtl, ireci) 
                      lpltp_collision = pltpenc_list%collision_check(system, param, t+dtl, dtl, ireci) 
 
-                     if (lplpl_collision) call pl%discard(system, param)
-                     if (lpltp_collision) call tp%discard(system, param)
+                     if (lplpl_collision) call plplenc_list%resolve_collision(system, param, t+dtl)
+                     if (lpltp_collision) call pltpenc_list%resolve_collision(system, param, t+dtl)
                   end if
 
                   call self%set_recur_levels(ireci)
diff --git a/src/symba/symba_util.f90 b/src/symba/symba_util.f90
index 05ee19f5e..8340f6e14 100644
--- a/src/symba/symba_util.f90
+++ b/src/symba/symba_util.f90
@@ -58,6 +58,29 @@ module subroutine symba_util_append_arr_kin(arr, source, nold, nsrc, lsource_mas
    end subroutine symba_util_append_arr_kin
 
 
+   module subroutine symba_util_append_encounter(self, source, lsource_mask)
+      !! author: David A. Minton
+      !!
+      !! Append components from one encounter list (pl-pl or pl-tp) body object to another. 
+      !! This method will automatically resize the destination body if it is too small
+      implicit none
+      ! Arguments
+      class(symba_encounter),    intent(inout) :: self         !! SyMBA encounter list object
+      class(swiftest_encounter), intent(in)    :: source       !! Source object to append
+      logical, dimension(:),     intent(in)    :: lsource_mask !! Logical mask indicating which elements to append to
+
+      associate(nold => self%nenc, nsrc => source%nenc)
+         select type(source)
+         class is (symba_encounter)
+            call util_append(self%level, source%level, nold, nsrc, lsource_mask)
+         end select
+         call util_append_encounter(self, source, lsource_mask) 
+      end associate
+
+      return
+   end subroutine symba_util_append_encounter
+
+
    module subroutine symba_util_append_pl(self, source, lsource_mask)
       !! author: David A. Minton
       !!
@@ -374,6 +397,8 @@ module subroutine symba_util_rearray_pl(self, system, param)
       class(symba_pl), allocatable :: tmp !! The discarded body list.
       integer(I4B) :: i
       logical, dimension(:), allocatable :: lmask
+      class(symba_plplenc), allocatable :: plplenc_old
+      logical :: lencounter
 
       associate(pl => self, pl_adds => system%pl_adds)
          allocate(tmp, mold=pl)
@@ -391,19 +416,24 @@ module subroutine symba_util_rearray_pl(self, system, param)
 
          ! Add in any new bodies
          if (pl_adds%nbody > 0) then
+            ! First store the original plplenc list so we don't remove any of the original encounters
+            allocate(plplenc_old, source=system%plplenc_list)
+
+            ! Append the adds to the main pl object
             call pl%append(pl_adds, lsource_mask=[(.true., i=1, pl_adds%nbody)])
+
             allocate(lmask(pl%nbody)) 
             lmask(:) = pl%status(1:pl%nbody) == NEW_PARTICLE
             call symba_io_dump_particle_info(system, param, plidx=pack([(i, i=1, pl%nbody)], lmask))
             where(pl%status(:) /= INACTIVE) pl%status(:) = ACTIVE
-         end if 
 
-         ! If there are still bodies in the system, sort by mass in descending order and re-index
-         if (pl%nbody > 0) then
             call pl%sort("mass", ascending=.false.)
             pl%lmtiny(:) = pl%Gmass(:) > param%GMTINY
             pl%nplm = count(pl%lmtiny(:))
+
+            ! Reindex
             call pl%eucl_index()
+
          end if
 
       end associate
@@ -637,6 +667,55 @@ module subroutine symba_util_sort_tp(self, sortby, ascending)
    end subroutine symba_util_sort_tp
 
 
+   module subroutine symba_util_sort_rearrange_arr_info(arr, ind, n)
+      !! author: David A. Minton
+      !!
+      !! Rearrange a single array of particle information type in-place from an index list.
+      implicit none
+      ! Arguments
+      type(symba_particle_info),  dimension(:), allocatable, intent(inout) :: arr !! Destination array 
+      integer(I4B),              dimension(:),              intent(in)    :: ind !! Index to rearrange against
+      integer(I4B),                            intent(in)    :: n   !! Number of elements in arr and ind to rearrange
+      ! Internals
+      type(symba_particle_info), dimension(:), allocatable                :: tmp !! Temporary copy of array used during rearrange operation
+
+      if (.not. allocated(arr) .or. n <= 0) return
+      allocate(tmp, mold=arr)
+      tmp(1:n) = arr(ind(1:n))
+      call move_alloc(tmp, arr)
+
+      return
+   end subroutine symba_util_sort_rearrange_arr_info
+
+
+   module subroutine symba_util_sort_rearrange_arr_kin(arr, ind, n)
+      !! author: David A. Minton
+      !!
+      !! Rearrange a single array of particle kinship type in-place from an index list.
+      implicit none
+      ! Arguments
+      type(symba_kinship),  dimension(:), allocatable, intent(inout) :: arr !! Destination array 
+      integer(I4B),         dimension(:),              intent(in)    :: ind !! Index to rearrange against
+      integer(I4B),                                    intent(in)    :: n   !! Number of elements in arr and ind to rearrange
+      ! Internals
+      type(symba_kinship),  dimension(:), allocatable                :: tmp !! Temporary copy of array used during rearrange operation
+      integer(I4B) :: i,j
+
+      if (.not. allocated(arr) .or. n <= 0) return
+      allocate(tmp, mold=arr)
+      tmp(1:n) = arr(ind(1:n))
+
+      do i = 1, n
+         do j = 1, tmp(i)%nchild
+            tmp(i)%child(j) = ind(tmp(i)%child(j))
+         end do
+      end do
+
+      call move_alloc(tmp, arr)
+      return
+   end subroutine symba_util_sort_rearrange_arr_kin
+
+
    module subroutine symba_util_sort_rearrange_pl(self, ind)
       !! author: David A. Minton
       !!
@@ -647,32 +726,23 @@ module subroutine symba_util_sort_rearrange_pl(self, ind)
       class(symba_pl),               intent(inout) :: self !! SyMBA massive body object
       integer(I4B),    dimension(:), intent(in)    :: ind  !! Index array used to restructure the body (should contain all 1:n index values in the desired order)
       ! Internals
-      class(symba_pl), allocatable :: pl_sorted  !! Temporary holder for sorted body
       integer(I4B) :: i, j
 
       associate(pl => self, npl => self%nbody)
+         call util_sort_rearrange(pl%lcollision, ind, npl)
+         call util_sort_rearrange(pl%lencounter, ind, npl)
+         call util_sort_rearrange(pl%lmtiny,     ind, npl)
+         call util_sort_rearrange(pl%nplenc,     ind, npl)
+         call util_sort_rearrange(pl%ntpenc,     ind, npl)
+         call util_sort_rearrange(pl%levelg,     ind, npl)
+         call util_sort_rearrange(pl%levelm,     ind, npl)
+         call util_sort_rearrange(pl%isperi,     ind, npl)
+         call util_sort_rearrange(pl%peri,       ind, npl)
+         call util_sort_rearrange(pl%atp,        ind, npl)
+         call util_sort_rearrange(pl%info,       ind, npl)
+         call util_sort_rearrange(pl%kin,        ind, npl)
+
          call util_sort_rearrange_pl(pl,ind)
-         allocate(pl_sorted, source=self)
-         if (allocated(pl%lcollision)) pl%lcollision(1:npl) = pl_sorted%lcollision(ind(1:npl))
-         if (allocated(pl%lencounter)) pl%lencounter(1:npl) = pl_sorted%lencounter(ind(1:npl))
-         if (allocated(pl%lmtiny))     pl%lmtiny(1:npl) = pl_sorted%lmtiny(ind(1:npl))
-         if (allocated(pl%nplenc))     pl%nplenc(1:npl) = pl_sorted%nplenc(ind(1:npl))
-         if (allocated(pl%ntpenc))     pl%ntpenc(1:npl) = pl_sorted%ntpenc(ind(1:npl))
-         if (allocated(pl%levelg))     pl%levelg(1:npl) = pl_sorted%levelg(ind(1:npl))
-         if (allocated(pl%levelm))     pl%levelm(1:npl) = pl_sorted%levelm(ind(1:npl))
-         if (allocated(pl%isperi))     pl%isperi(1:npl) = pl_sorted%isperi(ind(1:npl))
-         if (allocated(pl%peri))       pl%peri(1:npl) = pl_sorted%peri(ind(1:npl))
-         if (allocated(pl%atp))        pl%atp(1:npl) = pl_sorted%atp(ind(1:npl))
-         if (allocated(pl%info))       pl%info(1:npl) = pl_sorted%info(ind(1:npl))
-         if (allocated(pl%kin)) then
-            pl%kin(1:npl) = pl_sorted%kin(ind(1:npl))
-            do i = 1, npl
-               do j = 1, pl%kin(i)%nchild
-                  pl%kin(i)%child(j) = ind(pl%kin(i)%child(j))
-               end do
-            end do
-         end if
-         deallocate(pl_sorted)
       end associate
 
       return
@@ -688,17 +758,14 @@ module subroutine symba_util_sort_rearrange_tp(self, ind)
       ! Arguments
       class(symba_tp),               intent(inout) :: self !! SyMBA test particle object
       integer(I4B),    dimension(:), intent(in)    :: ind  !! Index array used to restructure the body (should contain all 1:n index values in the desired order)
-      ! Internals
-      class(symba_tp), allocatable :: tp_sorted  !! Temporary holder for sorted body
 
       associate(tp => self, ntp => self%nbody)
+         call util_sort_rearrange(tp%nplenc, ind, ntp)
+         call util_sort_rearrange(tp%levelg, ind, ntp)
+         call util_sort_rearrange(tp%levelm, ind, ntp)
+         call util_sort_rearrange(tp%info,   ind, ntp)
+
          call util_sort_rearrange_tp(tp,ind)
-         allocate(tp_sorted, source=self)
-         if (allocated(tp%nplenc)) tp%nplenc(1:ntp) = tp_sorted%nplenc(ind(1:ntp))
-         if (allocated(tp%levelg)) tp%levelg(1:ntp) = tp_sorted%levelg(ind(1:ntp))
-         if (allocated(tp%levelm)) tp%levelm(1:ntp) = tp_sorted%levelm(ind(1:ntp))
-         if (allocated(tp%info))   tp%info(1:ntp)   = tp_sorted%info(ind(1:ntp))
-         deallocate(tp_sorted)
       end associate
       
       return
diff --git a/src/util/util_append.f90 b/src/util/util_append.f90
index cb15fa18c..6f35de54e 100644
--- a/src/util/util_append.f90
+++ b/src/util/util_append.f90
@@ -184,6 +184,34 @@ module subroutine util_append_body(self, source, lsource_mask)
    end subroutine util_append_body
 
 
+   module subroutine util_append_encounter(self, source, lsource_mask)
+      !! author: David A. Minton
+      !!
+      !! Append components from one Swiftest body object to another. 
+      !! This method will automatically resize the destination body if it is too small
+      implicit none
+      ! Arguments
+      class(swiftest_encounter), intent(inout) :: self         !! Swiftest encounter list object
+      class(swiftest_encounter), intent(in)    :: source       !! Source object to append
+      logical, dimension(:),     intent(in)    :: lsource_mask !! Logical mask indicating which elements to append to
+
+      associate(nold => self%nenc, nsrc => source%nenc)
+         call util_append(self%lvdotr, source%lvdotr, nold, nsrc, lsource_mask)
+         call util_append(self%status, source%status, nold, nsrc, lsource_mask)
+         call util_append(self%index1, source%index1, nold, nsrc, lsource_mask)
+         call util_append(self%index2, source%index2, nold, nsrc, lsource_mask)
+         call util_append(self%x1, source%x1, nold, nsrc, lsource_mask)
+         call util_append(self%x2, source%x2, nold, nsrc, lsource_mask)
+         call util_append(self%v1, source%v1, nold, nsrc, lsource_mask)
+         call util_append(self%v2, source%v2, nold, nsrc, lsource_mask)
+         call util_append(self%t, source%t, nold, nsrc, lsource_mask)
+         self%nenc = nold + count(lsource_mask(:))
+      end associate
+
+      return
+   end subroutine util_append_encounter
+
+
    module subroutine util_append_pl(self, source, lsource_mask)
       !! author: David A. Minton
       !!
diff --git a/src/util/util_sort.f90 b/src/util/util_sort.f90
index 752e78ab7..b2a5464aa 100644
--- a/src/util/util_sort.f90
+++ b/src/util/util_sort.f90
@@ -55,6 +55,175 @@ module subroutine util_sort_body(self, sortby, ascending)
    end subroutine util_sort_body
 
 
+
+   module subroutine util_sort_dp(arr)
+      !! author: David A. Minton
+      !!
+      !! Sort input double precision array in place into ascending numerical order using insertion sort.
+      !! This algorithm works well for partially sorted arrays (which is usually the case here)
+      !!
+      implicit none
+      ! Arguments
+      real(DP), dimension(:), intent(inout) :: arr
+      ! Internals
+      real(DP)                :: tmp
+      integer(I4B)            :: n, i, j
+
+      n = size(arr)
+      do i = 2, n
+         tmp = arr(i)
+         do j = i - 1, 1, -1
+            if (arr(j) <= tmp) exit
+            arr(j + 1) = arr(j)
+         end do
+         arr(j + 1) = tmp
+      end do
+
+      return
+   end subroutine util_sort_dp
+
+
+   module subroutine util_sort_index_dp(arr, ind)
+      !! author: David A. Minton
+      !!
+      !! Sort input double precision array by index in ascending numerical order using insertion sort.
+      !! This algorithm works well for partially sorted arrays (which is usually the case here)
+      !!
+      implicit none
+      ! Arguments
+      real(DP), dimension(:), intent(in)  :: arr
+      integer(I4B), dimension(:), intent(out) :: ind
+      ! Internals
+      real(DP) :: tmp
+      integer(I4B) :: n, i, j
+
+      n = size(arr)
+      ind = [(i, i=1, n)]
+      do i = 2, n
+         tmp = arr(ind(i))
+         do j = i - 1, 1, -1
+            if (arr(ind(j)) <= tmp) exit
+            ind(j + 1) = ind(j)
+         end do
+         ind(j + 1) = i
+      end do
+
+      return
+   end subroutine util_sort_index_dp
+
+
+   module subroutine util_sort_i4b(arr)
+      !! author: David A. Minton
+      !!
+      !! Sort input integer array in place into ascending numerical order using insertion sort.
+      !! This algorithm works well for partially sorted arrays (which is usually the case here)
+      !!
+      implicit none
+      ! Arguments
+      integer(I4B), dimension(:), intent(inout) :: arr
+      ! Internals
+      integer(I4B) :: tmp
+      integer(I4B) :: n, i, j
+
+      n = size(arr)
+      do i = 2, n
+         tmp = arr(i)
+         do j = i - 1, 1, -1
+            if (arr(j) <= tmp) exit
+            arr(j + 1) = arr(j)
+         end do
+         arr(j + 1) = tmp
+      end do
+
+      return
+   end subroutine util_sort_i4b
+
+
+   module subroutine util_sort_index_i4b(arr, ind)
+      !! author: David A. Minton
+      !!
+      !! Sort input integer array by index in ascending numerical order using insertion sort.
+      !! This algorithm works well for partially sorted arrays (which is usually the case here)
+      !!
+      implicit none
+      ! Arguments
+      integer(I4B), dimension(:), intent(in)  :: arr
+      integer(I4B), dimension(:), intent(out) :: ind
+      ! Internals
+      integer(I4B) :: tmp
+      integer(I4B) :: n, i, j
+
+      n = size(arr)
+      ind = [(i, i=1, n)]
+      do i = 2, n
+         tmp = arr(ind(i))
+         do j = i - 1, 1, -1
+            if (arr(ind(j)) <= tmp) exit
+            ind(j + 1) = ind(j)
+         end do
+         ind(j + 1) = i
+      end do
+
+      return
+   end subroutine util_sort_index_i4b
+
+
+   module subroutine util_sort_sp(arr)
+      !! author: David A. Minton
+      !!
+      !! Sort input single precision array in place into ascending numerical order using insertion sort.
+      !! This algorithm works well for partially sorted arrays (which is usually the case here)
+      !
+      implicit none
+      ! Arguments
+      real(SP), dimension(:), intent(inout) :: arr
+      ! Internals
+      real(SP) :: tmp
+      integer(I4B) :: n, i, j
+
+      n = size(arr)
+      do i = 2, n
+         tmp = arr(i)
+         do j = i - 1, 1, -1
+            if (arr(j) <= tmp) exit
+            arr(j + 1) = arr(j)
+         end do
+         arr(j + 1) = tmp
+      end do
+
+      return
+   end subroutine util_sort_sp
+
+
+   module subroutine util_sort_index_sp(arr, ind)
+      !! author: David A. Minton
+      !!
+      !! Sort input single precision array by index in ascending numerical order using insertion sort.
+      !! This algorithm works well for partially sorted arrays (which is usually the case here)
+      !!
+      implicit none
+      ! Arguments
+      real(SP), dimension(:), intent(in)  :: arr
+      integer(I4B), dimension(:), intent(out) :: ind
+      ! Internals
+      real(SP) :: tmp
+      integer(I4B) :: n, i, j
+
+      n = size(arr)
+      ind = [(i, i=1, n)]
+      do i = 2, n
+         tmp = arr(ind(i))
+         do j = i - 1, 1, -1
+            if (arr(ind(j)) <= tmp) exit
+            ind(j + 1) = ind(j)
+         end do
+         ind(j + 1) = i
+      end do
+
+      return
+   end subroutine util_sort_index_sp
+
+
    module subroutine util_sort_pl(self, sortby, ascending)
       !! author: David A. Minton
       !!
@@ -156,264 +325,189 @@ module subroutine util_sort_rearrange_body(self, ind)
       ! Arguments
       class(swiftest_body),               intent(inout) :: self !! Swiftest body object
       integer(I4B),         dimension(:), intent(in)    :: ind  !! Index array used to restructure the body (should contain all 1:n index values in the desired order)
-      ! Internals
-      class(swiftest_body), allocatable :: body_sorted  !! Temporary holder for sorted body
 
       associate(n => self%nbody)
-         allocate(body_sorted, source=self)
-         if (allocated(self%id))       self%id(1:n) = body_sorted%id(ind(1:n))
-         if (allocated(self%name))     self%name(1:n) = body_sorted%name(ind(1:n))
-         if (allocated(self%status))   self%status(1:n) = body_sorted%status(ind(1:n))
-         if (allocated(self%ldiscard)) self%ldiscard(1:n) = body_sorted%ldiscard(ind(1:n))
-         if (allocated(self%xh))       self%xh(:,1:n) = body_sorted%xh(:,ind(1:n))
-         if (allocated(self%vh))       self%vh(:,1:n) = body_sorted%vh(:,ind(1:n))
-         if (allocated(self%xb))       self%xb(:,1:n) = body_sorted%xb(:,ind(1:n))
-         if (allocated(self%vb))       self%vb(:,1:n) = body_sorted%vb(:,ind(1:n))
-         if (allocated(self%ah))       self%ah(:,1:n) = body_sorted%ah(:,ind(1:n))
-         if (allocated(self%ir3h))     self%ir3h(1:n) = body_sorted%ir3h(ind(1:n))
-         if (allocated(self%mu))       self%mu(1:n) = body_sorted%mu(ind(1:n))
-         if (allocated(self%lmask))    self%lmask(1:n) = body_sorted%lmask(ind(1:n))
-         if (allocated(self%a))        self%a(1:n) = body_sorted%a(ind(1:n))
-         if (allocated(self%e))        self%e(1:n) = body_sorted%e(ind(1:n))
-         if (allocated(self%inc))      self%inc(1:n) = body_sorted%inc(ind(1:n))
-         if (allocated(self%capom))    self%capom(1:n) = body_sorted%capom(ind(1:n))
-         if (allocated(self%omega))    self%omega(1:n) = body_sorted%omega(ind(1:n))
-         if (allocated(self%capm))     self%capm(1:n) = body_sorted%capm(ind(1:n))
-         if (allocated(self%aobl))     self%aobl(:,1:n) = body_sorted%aobl(:,ind(1:n))
-         if (allocated(self%atide))    self%atide(:,1:n) = body_sorted%atide(:,ind(1:n))
-         if (allocated(self%agr))      self%agr(:,1:n) = body_sorted%agr(:,ind(1:n))
-         deallocate(body_sorted)
+         call util_sort_rearrange(self%id,       ind, n)
+         call util_sort_rearrange(self%name,     ind, n)
+         call util_sort_rearrange(self%status,   ind, n)
+         call util_sort_rearrange(self%ldiscard, ind, n)
+         call util_sort_rearrange(self%xh,       ind, n)
+         call util_sort_rearrange(self%vh,       ind, n)
+         call util_sort_rearrange(self%xb,       ind, n)
+         call util_sort_rearrange(self%vb,       ind, n)
+         call util_sort_rearrange(self%ah,       ind, n)
+         call util_sort_rearrange(self%ir3h,     ind, n)
+         call util_sort_rearrange(self%mu,       ind, n)
+         call util_sort_rearrange(self%lmask,    ind, n)
+         call util_sort_rearrange(self%a,        ind, n)
+         call util_sort_rearrange(self%e,        ind, n)
+         call util_sort_rearrange(self%inc,      ind, n)
+         call util_sort_rearrange(self%capom,    ind, n)
+         call util_sort_rearrange(self%omega,    ind, n)
+         call util_sort_rearrange(self%capm,     ind, n)
+         call util_sort_rearrange(self%aobl,     ind, n)
+         call util_sort_rearrange(self%atide,    ind, n)
+         call util_sort_rearrange(self%agr,      ind, n)
       end associate
 
       return
    end subroutine util_sort_rearrange_body
 
 
-   module subroutine util_sort_rearrange_pl(self, ind)
+   module subroutine util_sort_rearrange_arr_char_string(arr, ind, n)
       !! author: David A. Minton
       !!
-      !! Rearrange Swiftest massive body structure in-place from an index list.
-      !! This is a helper utility used to make polymorphic sorting work on Swiftest structures.
-      implicit none
-      class(swiftest_pl),               intent(inout) :: self !! Swiftest massive body object
-      integer(I4B),       dimension(:), intent(in)    :: ind  !! Index array used to restructure the body (should contain all 1:n index values in the desired order)
-      ! Internals
-      class(swiftest_pl), allocatable :: pl_sorted  !! Temporary holder for sorted body
-
-      associate(pl => self, npl => self%nbody)
-         call util_sort_rearrange_body(pl,ind)
-         allocate(pl_sorted, source=self)
-         if (allocated(pl%mass))    pl%mass(1:npl) = pl_sorted%mass(ind(1:npl))
-         if (allocated(pl%Gmass))   pl%Gmass(1:npl) = pl_sorted%Gmass(ind(1:npl))
-         if (allocated(pl%rhill))   pl%rhill(1:npl) = pl_sorted%rhill(ind(1:npl))
-         if (allocated(pl%xbeg))    pl%xbeg(:,1:npl) = pl_sorted%xbeg(:,ind(1:npl))
-         if (allocated(pl%xend))    pl%xend(:,1:npl) = pl_sorted%xend(:,ind(1:npl))
-         if (allocated(pl%vbeg))    pl%vbeg(:,1:npl) = pl_sorted%vbeg(:,ind(1:npl))
-         if (allocated(pl%radius))  pl%radius(1:npl) = pl_sorted%radius(ind(1:npl))
-         if (allocated(pl%density)) pl%density(1:npl) = pl_sorted%density(ind(1:npl))
-         if (allocated(pl%Ip))      pl%Ip(:,1:npl) = pl_sorted%Ip(:,ind(1:npl))
-         if (allocated(pl%rot))     pl%rot(:,1:npl) = pl_sorted%rot(:,ind(1:npl))
-         if (allocated(pl%k2))      pl%k2(1:npl) = pl_sorted%k2(ind(1:npl))
-         if (allocated(pl%Q))       pl%Q(1:npl) = pl_sorted%Q(ind(1:npl))
-         if (allocated(pl%tlag))    pl%tlag(1:npl) = pl_sorted%tlag(ind(1:npl))
-
-         deallocate(pl_sorted)
-      end associate
-
-      return
-   end subroutine util_sort_rearrange_pl
-
-
-   module subroutine util_sort_rearrange_tp(self, ind)
-      !! author: David A. Minton
-      !!
-      !! Rearrange Swiftest massive body structure in-place from an index list.
-      !! This is a helper utility used to make polymorphic sorting work on Swiftest structures.
+      !! Rearrange a single array of character string in-place from an index list.
       implicit none
       ! Arguments
-      class(swiftest_tp),                 intent(inout) :: self !! Swiftest test particle object
-      integer(I4B),         dimension(:), intent(in)    :: ind  !! Index array used to restructure the body (should contain all 1:n index values in the desired order)
+      character(len=STRMAX), dimension(:), allocatable, intent(inout) :: arr !! Destination array 
+      integer(I4B),          dimension(:),              intent(in)    :: ind !! Index to rearrange against
+      integer(I4B),                                     intent(in)    :: n   !! Number of elements in arr and ind to rearrange
       ! Internals
-      class(swiftest_tp), allocatable :: tp_sorted  !! Temporary holder for sorted body
+      character(len=STRMAX), dimension(:), allocatable                :: tmp !! Temporary copy of arry used during rearrange operation
 
-      associate(tp => self, ntp => self%nbody)
-         call util_sort_rearrange_body(tp,ind)
-         allocate(tp_sorted, source=self)
-         if (allocated(tp%isperi))  tp%isperi(1:ntp) = tp_sorted%isperi(ind(1:ntp))
-         if (allocated(tp%peri))    tp%peri(1:ntp) = tp_sorted%peri(ind(1:ntp))
-         if (allocated(tp%atp))     tp%atp(1:ntp) = tp_sorted%atp(ind(1:ntp))
-         deallocate(tp_sorted)
-      end associate
+      if (.not. allocated(arr) .or. n <= 0) return
+      allocate(tmp, mold=arr)
+      tmp(1:n) = arr(ind(1:n))
+      call move_alloc(tmp, arr)
 
       return
-   end subroutine util_sort_rearrange_tp
+   end subroutine util_sort_rearrange_arr_char_string
 
 
-   module subroutine util_sort_dp(arr)
+   module subroutine util_sort_rearrange_arr_DP(arr, ind, n)
       !! author: David A. Minton
       !!
-      !! Sort input double precision array in place into ascending numerical order using insertion sort.
-      !! This algorithm works well for partially sorted arrays (which is usually the case here)
-      !!
+      !! Rearrange a single array of DP type in-place from an index list.
       implicit none
       ! Arguments
-      real(DP), dimension(:), intent(inout) :: arr
+      real(DP),     dimension(:), allocatable, intent(inout) :: arr !! Destination array 
+      integer(I4B), dimension(:),              intent(in)  :: ind !! Index to rearrange against
+      integer(I4B),                            intent(in)  :: n   !! Number of elements in arr and ind to rearrange
       ! Internals
-      real(DP)                :: tmp
-      integer(I4B)            :: n, i, j
+      real(DP), dimension(:), allocatable :: tmp !! Temporary copy of array used during rearrange operation
 
-      n = size(arr)
-      do i = 2, n
-         tmp = arr(i)
-         do j = i - 1, 1, -1
-            if (arr(j) <= tmp) exit
-            arr(j + 1) = arr(j)
-         end do
-         arr(j + 1) = tmp
-      end do
+      if (.not. allocated(arr) .or. n <= 0) return
+      allocate(tmp, mold=arr)
+      tmp(1:n) = arr(ind(1:n))
+      call move_alloc(tmp, arr)
 
       return
-   end subroutine util_sort_dp
+   end subroutine util_sort_rearrange_arr_DP
 
 
-   module subroutine util_sort_index_dp(arr, ind)
+   module subroutine util_sort_rearrange_arr_DPvec(arr, ind, n)
       !! author: David A. Minton
       !!
-      !! Sort input double precision array by index in ascending numerical order using insertion sort.
-      !! This algorithm works well for partially sorted arrays (which is usually the case here)
-      !!
+      !! Rearrange a single array of (NDIM,n) DP-type vectors in-place from an index list.
       implicit none
       ! Arguments
-      real(DP), dimension(:), intent(in)  :: arr
-      integer(I4B), dimension(:), intent(out) :: ind
+      real(DP),     dimension(:,:), allocatable, intent(inout) :: arr !! Destination array 
+      integer(I4B), dimension(:),                intent(in)    :: ind !! Index to rearrange against
+      integer(I4B),                              intent(in)    :: n   !! Number of elements in arr and ind to rearrange
       ! Internals
-      real(DP) :: tmp
-      integer(I4B) :: n, i, j
+      real(DP), dimension(:,:), allocatable :: tmp !! Temporary copy of array used during rearrange operation
 
-      n = size(arr)
-      ind = [(i, i=1, n)]
-      do i = 2, n
-         tmp = arr(ind(i))
-         do j = i - 1, 1, -1
-            if (arr(ind(j)) <= tmp) exit
-            ind(j + 1) = ind(j)
-         end do
-         ind(j + 1) = i
-      end do
+      if (.not. allocated(arr) .or. n <= 0) return
+      allocate(tmp, mold=arr)
+      tmp(:,1:n) = arr(:, ind(1:n))
+      call move_alloc(tmp, arr)
 
       return
-   end subroutine util_sort_index_dp
+   end subroutine util_sort_rearrange_arr_DPvec
 
 
-   module subroutine util_sort_i4b(arr)
+   module subroutine util_sort_rearrange_arr_I4B(arr, ind, n)
       !! author: David A. Minton
       !!
-      !! Sort input integer array in place into ascending numerical order using insertion sort.
-      !! This algorithm works well for partially sorted arrays (which is usually the case here)
-      !!
+      !! Rearrange a single array of integers in-place from an index list.
       implicit none
       ! Arguments
-      integer(I4B), dimension(:), intent(inout) :: arr
+      integer(I4B), dimension(:), allocatable, intent(inout) :: arr !! Destination array 
+      integer(I4B), dimension(:),              intent(in)    :: ind !! Index to rearrange against
+      integer(I4B),                             intent(in)    :: n   !! Number of elements in arr and ind to rearrange
       ! Internals
-      integer(I4B) :: tmp
-      integer(I4B) :: n, i, j
+      integer(I4B), dimension(:), allocatable                :: tmp !! Temporary copy of array used during rearrange operation
 
-      n = size(arr)
-      do i = 2, n
-         tmp = arr(i)
-         do j = i - 1, 1, -1
-            if (arr(j) <= tmp) exit
-            arr(j + 1) = arr(j)
-         end do
-         arr(j + 1) = tmp
-      end do
+      if (.not. allocated(arr) .or. n <= 0) return
+      allocate(tmp, mold=arr)
+      tmp(1:n) = arr(ind(1:n))
+      call move_alloc(tmp, arr)
 
       return
-   end subroutine util_sort_i4b
+   end subroutine util_sort_rearrange_arr_I4B
 
 
-   module subroutine util_sort_index_i4b(arr, ind)
+   module subroutine util_sort_rearrange_arr_logical(arr, ind, n)
       !! author: David A. Minton
       !!
-      !! Sort input integer array by index in ascending numerical order using insertion sort.
-      !! This algorithm works well for partially sorted arrays (which is usually the case here)
-      !!
+      !! Rearrange a single array of logicals in-place from an index list.
       implicit none
       ! Arguments
-      integer(I4B), dimension(:), intent(in)  :: arr
-      integer(I4B), dimension(:), intent(out) :: ind
+      logical,      dimension(:), allocatable, intent(inout) :: arr !! Destination array 
+      integer(I4B), dimension(:),              intent(in)    :: ind !! Index to rearrange against
+      integer(I4B),                            intent(in)    :: n   !! Number of elements in arr and ind to rearrange
       ! Internals
-      integer(I4B) :: tmp
-      integer(I4B) :: n, i, j
+      logical, dimension(:), allocatable                :: tmp !! Temporary copy of array used during rearrange operation
 
-      n = size(arr)
-      ind = [(i, i=1, n)]
-      do i = 2, n
-         tmp = arr(ind(i))
-         do j = i - 1, 1, -1
-            if (arr(ind(j)) <= tmp) exit
-            ind(j + 1) = ind(j)
-         end do
-         ind(j + 1) = i
-      end do
+      if (.not. allocated(arr) .or. n <= 0) return
+      allocate(tmp, mold=arr)
+      tmp(1:n) = arr(ind(1:n))
+      call move_alloc(tmp, arr)
 
       return
-   end subroutine util_sort_index_i4b
+   end subroutine util_sort_rearrange_arr_logical
 
 
-   module subroutine util_sort_sp(arr)
+   module subroutine util_sort_rearrange_pl(self, ind)
       !! author: David A. Minton
       !!
-      !! Sort input single precision array in place into ascending numerical order using insertion sort.
-      !! This algorithm works well for partially sorted arrays (which is usually the case here)
-      !
+      !! Rearrange Swiftest massive body structure in-place from an index list.
+      !! This is a helper utility used to make polymorphic sorting work on Swiftest structures.
       implicit none
-      ! Arguments
-      real(SP), dimension(:), intent(inout) :: arr
-      ! Internals
-      real(SP) :: tmp
-      integer(I4B) :: n, i, j
+      class(swiftest_pl),               intent(inout) :: self !! Swiftest massive body object
+      integer(I4B),       dimension(:), intent(in)    :: ind  !! Index array used to restructure the body (should contain all 1:n index values in the desired order)
 
-      n = size(arr)
-      do i = 2, n
-         tmp = arr(i)
-         do j = i - 1, 1, -1
-            if (arr(j) <= tmp) exit
-            arr(j + 1) = arr(j)
-         end do
-         arr(j + 1) = tmp
-      end do
+      associate(pl => self, npl => self%nbody)
+         call util_sort_rearrange(pl%mass,    ind, npl)
+         call util_sort_rearrange(pl%Gmass,   ind, npl)
+         call util_sort_rearrange(pl%rhill,   ind, npl)
+         call util_sort_rearrange(pl%xbeg,    ind, npl)
+         call util_sort_rearrange(pl%vbeg,    ind, npl)
+         call util_sort_rearrange(pl%radius,  ind, npl)
+         call util_sort_rearrange(pl%density, ind, npl)
+         call util_sort_rearrange(pl%Ip,      ind, npl)
+         call util_sort_rearrange(pl%rot,     ind, npl)
+         call util_sort_rearrange(pl%k2,      ind, npl)
+         call util_sort_rearrange(pl%Q,       ind, npl)
+         call util_sort_rearrange(pl%tlag,    ind, npl)
+
+         call util_sort_rearrange_body(pl, ind)
+      end associate
 
       return
-   end subroutine util_sort_sp
+   end subroutine util_sort_rearrange_pl
 
 
-   module subroutine util_sort_index_sp(arr, ind)
+   module subroutine util_sort_rearrange_tp(self, ind)
       !! author: David A. Minton
       !!
-      !! Sort input single precision array by index in ascending numerical order using insertion sort.
-      !! This algorithm works well for partially sorted arrays (which is usually the case here)
-      !!
+      !! Rearrange Swiftest massive body structure in-place from an index list.
+      !! This is a helper utility used to make polymorphic sorting work on Swiftest structures.
       implicit none
       ! Arguments
-      real(SP), dimension(:), intent(in)  :: arr
-      integer(I4B), dimension(:), intent(out) :: ind
-      ! Internals
-      real(SP) :: tmp
-      integer(I4B) :: n, i, j
+      class(swiftest_tp),                 intent(inout) :: self !! Swiftest test particle object
+      integer(I4B),         dimension(:), intent(in)    :: ind  !! Index array used to restructure the body (should contain all 1:n index values in the desired order)
 
-      n = size(arr)
-      ind = [(i, i=1, n)]
-      do i = 2, n
-         tmp = arr(ind(i))
-         do j = i - 1, 1, -1
-            if (arr(ind(j)) <= tmp) exit
-            ind(j + 1) = ind(j)
-         end do
-         ind(j + 1) = i
-      end do
+      associate(tp => self, ntp => self%nbody)
+         call util_sort_rearrange(tp%isperi, ind, ntp)
+         call util_sort_rearrange(tp%peri,   ind, ntp)
+         call util_sort_rearrange(tp%atp,    ind, ntp)
+
+         call util_sort_rearrange_body(tp, ind)
+      end associate
 
       return
-   end subroutine util_sort_index_sp
+   end subroutine util_sort_rearrange_tp
 
 end submodule s_util_sort
diff --git a/src/whm/whm_util.f90 b/src/whm/whm_util.f90
index cc84ba3d5..537866d0e 100644
--- a/src/whm/whm_util.f90
+++ b/src/whm/whm_util.f90
@@ -165,21 +165,17 @@ module subroutine whm_util_sort_rearrange_pl(self, ind)
       ! Arguments
       class(whm_pl),               intent(inout) :: self !! WHM massive body object
       integer(I4B),  dimension(:), intent(in)    :: ind  !! Index array used to restructure the body (should contain all 1:n index values in the desired order)
-      ! Internals
-      class(whm_pl), allocatable :: pl_sorted  !! Temporary holder for sorted body
-      integer(I4B) :: i
 
       if (self%nbody == 0) return
 
       associate(pl => self, npl => self%nbody)
+         call util_sort_rearrange(pl%eta,  ind, npl)
+         call util_sort_rearrange(pl%xj,   ind, npl)
+         call util_sort_rearrange(pl%vj,   ind, npl)
+         call util_sort_rearrange(pl%muj,  ind, npl)
+         call util_sort_rearrange(pl%ir3j, ind, npl)
+
          call util_sort_rearrange_pl(pl,ind)
-         allocate(pl_sorted, source=self)
-         if (allocated(pl%eta))  pl%eta(1:npl) = pl_sorted%eta(ind(1:npl))
-         if (allocated(pl%xj))   pl%xj(:,1:npl) = pl_sorted%xj(:,ind(1:npl))
-         if (allocated(pl%vj))   pl%vj(:,1:npl) = pl_sorted%vj(:,ind(1:npl))
-         if (allocated(pl%muj))  pl%muj(1:npl) = pl_sorted%muj(ind(1:npl))
-         if (allocated(pl%ir3j)) pl%ir3j(1:npl) = pl_sorted%ir3j(ind(1:npl))
-         deallocate(pl_sorted)
       end associate
 
       return