From f5632252daacd54c49c97930d24aa9467f41e44c Mon Sep 17 00:00:00 2001
From: David A Minton <daminton@purdue.edu>
Date: Mon, 30 Aug 2021 20:38:16 -0400
Subject: [PATCH] Rearranged how the non-plpl discards work to use the same
 pattern as the plpl discards

---
 src/symba/symba_discard.f90 | 31 ++++++++++++++++++++++---------
 1 file changed, 22 insertions(+), 9 deletions(-)

diff --git a/src/symba/symba_discard.f90 b/src/symba/symba_discard.f90
index 97a7f84ac..24f40b1bd 100644
--- a/src/symba/symba_discard.f90
+++ b/src/symba/symba_discard.f90
@@ -184,19 +184,32 @@ subroutine symba_discard_nonplpl(pl, system, param)
       class(swiftest_parameters),   intent(in)    :: param  !! Current run configuration parameters 
       ! Internals
       logical, dimension(pl%nbody) ::  ldiscard
+      integer(I4B) :: i, nstart, nend, nsub
+      class(symba_pl), allocatable            :: plsub
     
       ! First check for collisions with the central body
       associate(npl => pl%nbody, cb => system%cb)
          if (npl == 0) return 
-         ldiscard(1:npl) = pl%ldiscard(1:npl) ! Don't include any bodies that were previously flagged for discard in here
-         if ((param%rmin >= 0.0_DP) .or. (param%rmax >= 0.0_DP) .or.  (param%rmaxu >= 0.0_DP)) then
-            call symba_discard_cb_pl(pl, system, param)
-         end if
-         if (param%qmin >= 0.0_DP .and. npl > 0) call symba_discard_peri_pl(pl, system, param)
-         if (any(.not.ldiscard(1:npl) .and. pl%ldiscard(1:npl))) then
-            ldiscard(1:npl) = .not.ldiscard(1:npl) .and. pl%ldiscard(1:npl)
-            call system%pl_discards%append(pl, ldiscard) 
-         end if
+         select type(pl_discards => system%pl_discards)
+         class is (symba_merger)
+            if ((param%rmin >= 0.0_DP) .or. (param%rmax >= 0.0_DP) .or.  (param%rmaxu >= 0.0_DP)) then
+               call symba_discard_cb_pl(pl, system, param)
+            end if
+            if (param%qmin >= 0.0_DP) call symba_discard_peri_pl(pl, system, param)
+            if (any(pl%ldiscard(1:npl))) then
+               ldiscard(1:npl) = pl%ldiscard(1:npl)
+                  
+               allocate(plsub, mold=pl)
+               call pl%spill(plsub, ldiscard, ldestructive=.false.)
+               nsub = plsub%nbody
+               nstart = pl_discards%nbody + 1
+               nend = pl_discards%nbody + nsub
+               call pl_discards%append(plsub, lsource_mask=[(.true., i = 1, nsub)])
+   
+               ! Record how many bodies were subtracted in this event
+               pl_discards%ncomp(nstart:nend) = nsub
+            end if
+         end select
       end associate
 
       return