Added a rearray_tp method to the swiftest_tp class in order to better…

… resolve pl-tp collisions

src/collision/collision_resolve.f90

@@ -220,7 +220,7 @@ module subroutine collision_resolve_extract_plpl(self, nbody_system, param)
             unique_parent_idx = pack(collision_idx(:), lplpl_unique_parent(:))
 
             ! Scrub all pl-pl collisions involving unique pairs of parents, which will remove all duplicates and leave behind
-            ! all pairs that have themselves as parents but are not part of the unique parent list. This can hapepn in rare cases
+            ! all pairs that have themselves as parents but are not part of the unique parent list. This can happen in rare cases
             ! due to restructuring of parent/child relationships when there are large numbers of multi-body collisions in a single
             ! step
             lplpl_unique_parent(:) = .true.
@@ -713,13 +713,19 @@ module subroutine collision_resolve_pltp(self, nbody_system, param, t, dt, irec)
       class is (swiftest_nbody_system)
       select type(param)
       class is (swiftest_parameters)
-         call nbody_system%pl%vb2vh(nbody_system%cb)
-         call nbody_system%tp%vb2vh(nbody_system%cb%vb)
-         call nbody_system%pl%b2h(nbody_system%cb)
-         call nbody_system%tp%b2h(nbody_system%cb)
-
-         ! Discard the collider
-         call nbody_system%tp%discard(nbody_system, param)
+         associate(pltp_collision => nbody_system%pltp_collision, &
+            collision_history => nbody_system%collision_history, pl => nbody_system%pl, cb => nbody_system%cb, &
+            tp => nbody_system%tp, collider => nbody_system%collider, impactors => nbody_system%collider%impactors)
+            call pl%vb2vh(nbody_system%cb)
+            call tp%vb2vh(nbody_system%cb%vb)
+            call pl%b2h(nbody_system%cb)
+            call tp%b2h(nbody_system%cb)
+
+            call tp%rearray(nbody_system, param)
+
+            ! Discard the collider
+            !call nbody_system%tp%discard(nbody_system, param)
+         end associate
       end select
       end select
 

src/swiftest/swiftest_module.f90

@@ -321,6 +321,7 @@ module swiftest
       procedure :: rh2rb     => swiftest_util_coord_rh2rb_tp    !! Convert test particles from heliocentric to barycentric coordinates (position only)
       procedure :: dealloc   => swiftest_util_dealloc_tp        !! Deallocates all allocatable arrays
       procedure :: fill      => swiftest_util_fill_tp           !! "Fills" bodies from one object into another depending on the results of a mask (uses the UNPACK intrinsic)
+      procedure :: rearray   => swiftest_util_rearray_tp        !! Clean up the test particle structures to remove discarded bodies 
       procedure :: resize    => swiftest_util_resize_tp         !! Checks the current size of a Swiftest body against the requested size and resizes it if it is too small.
       procedure :: set_mu    => swiftest_util_set_mu_tp         !! Method used to construct the vectorized form of the central body mass
       procedure :: sort      => swiftest_util_sort_tp           !! Sorts body arrays by a sortable component
@@ -1485,11 +1486,18 @@ end subroutine swiftest_util_peri_tp
 
       module subroutine swiftest_util_rearray_pl(self, nbody_system, param)
          implicit none
-         class(swiftest_pl),           intent(inout) :: self   !! SyMBA massive body object
-         class(swiftest_nbody_system), intent(inout) :: nbody_system !! SyMBA nbody system object
-         class(swiftest_parameters),   intent(inout) :: param  !! Current run configuration parameters with SyMBA additions
+         class(swiftest_pl),           intent(inout) :: self   !! Swiftest massive body object
+         class(swiftest_nbody_system), intent(inout) :: nbody_system !! Swiftest nbody system object
+         class(swiftest_parameters),   intent(inout) :: param  !! Current run configuration parameters 
       end subroutine swiftest_util_rearray_pl
 
+      module subroutine swiftest_util_rearray_tp(self, nbody_system, param)
+         implicit none
+         class(swiftest_tp),           intent(inout) :: self   !! Swiftest test particle object
+         class(swiftest_nbody_system), intent(inout) :: nbody_system !! Swiftest nbody system object
+         class(swiftest_parameters),   intent(inout) :: param  !! Current run configuration parameters 
+      end subroutine swiftest_util_rearray_tp
+
       module subroutine swiftest_util_rescale_system(self, param, mscale, dscale, tscale)
          implicit none
          class(swiftest_nbody_system), intent(inout) :: self   !! Swiftest nbody system object

src/swiftest/swiftest_util.f90

@@ -1704,10 +1704,13 @@ module subroutine swiftest_util_rearray_pl(self, nbody_system, param)
 
          if (npl == 0) then
             if (param%lmtiny_pl) pl%nplm = 0
+            ! There are no more massive bodies. Reset the encounter lists and move on
+            if (allocated(nbody_system%plpl_encounter)) call nbody_system%plpl_encounter%setup(0_I8B)
+            if (allocated(nbody_system%pltp_encounter)) call nbody_system%pltp_encounter%setup(0_I8B)
             return
          end if
 
-         ! Reset all of the status flags for this body
+         ! Reset all of the status flags for the remaining bodies
          pl%status(1:npl) = ACTIVE
          do i = 1, npl
             call pl%info(i)%set_value(status="ACTIVE")
@@ -1842,6 +1845,82 @@ module subroutine swiftest_util_rearray_pl(self, nbody_system, param)
    end subroutine swiftest_util_rearray_pl
 
 
+
+   module subroutine swiftest_util_rearray_tp(self, nbody_system, param)
+      !! Author: David A. Minton
+      !!
+      !! Clean up the test particle structures to remove discarded bodies
+      use symba
+      implicit none
+      ! Arguments
+      class(swiftest_tp),           intent(inout) :: self   !! Swiftest test particle object
+      class(swiftest_nbody_system), intent(inout) :: nbody_system !! Swiftest nbody system object
+      class(swiftest_parameters),   intent(inout) :: param  !! Current run configuration parameters 
+      ! Internals
+      class(swiftest_tp), allocatable :: tmp !! The discarded body list.
+      integer(I4B) :: i, ntp
+      integer(I8B) :: k, nenc_old, nencmin
+      logical, dimension(:), allocatable :: lmask
+      class(encounter_list), allocatable :: pltpenc_old
+      logical :: lencounter
+
+      associate(pl => self, tp => nbody_system%tp, cb => nbody_system%cb, pl_adds => nbody_system%pl_adds)
+
+         ntp = tp%nbody
+         if (ntp == 0) return
+
+         ! Remove the discards and destroy the list, as the nbody_system already tracks tp_discards elsewhere
+         allocate(lmask(ntp))
+         lmask(1:ntp) = tp%ldiscard(1:ntp)
+         if (count(lmask(:)) > 0) then
+            allocate(tmp, mold=self)
+            call tp%spill(tmp, lspill_list=lmask, ldestructive=.true.)
+            ntp = tp%nbody
+            call tmp%setup(0,param)
+            deallocate(tmp)
+            deallocate(lmask)
+         end if
+         ntp = tp%nbody
+         if (ntp == 0) then
+            ! There are no more test particles. Reset the encounter list and move on
+            if (allocated(nbody_system%pltp_encounter)) call nbody_system%pltp_encounter%setup(0_I8B)
+            return
+         end if
+
+         ! Reset all of thes tatus flags for the remaining bodies
+         tp%status(1:ntp) = ACTIVE
+         do i = 1, ntp
+            call tp%info(i)%set_value(status="ACTIVE")
+         end do
+         tp%ldiscard(1:ntp) = .false.
+         tp%lcollision(1:ntp) = .false.
+         tp%lmask(1:ntp) = .true.
+
+         if (allocated(nbody_system%pltp_encounter)) then
+            ! Store the original pltpenc list so we don't remove any of the original encounters
+            nenc_old = nbody_system%pltp_encounter%nenc
+            if (nenc_old > 0_I8B) then 
+               allocate(pltpenc_old, source=nbody_system%pltp_encounter)
+               call pltpenc_old%copy(nbody_system%pltp_encounter)
+            end if
+
+            ! Re-build the encounter list
+            ! Be sure to get the level info if this is a SyMBA nbody_system
+            select type(nbody_system)
+            class is (symba_nbody_system)
+            select type(tp)
+            class is (symba_tp)
+               lencounter = tp%encounter_check(param, nbody_system, param%dt, nbody_system%irec)
+            end select
+            end select
+         end if
+
+      end associate
+
+      return
+   end subroutine swiftest_util_rearray_tp
+
+
    module subroutine swiftest_util_rescale_system(self, param, mscale, dscale, tscale)
       !! author: David A. Minton
       !!