Added symba_discard_conserve_mtm code but have not troubleshooted all…

… of the probems yet

src/modules/symba_classes.f90

@@ -95,6 +95,7 @@ module symba_classes
       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 :: 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
@@ -533,6 +534,13 @@ module subroutine symba_util_peri_pl(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 symba_util_peri_pl
+
+      module subroutine symba_util_rearray_pl(self, system, param)
+         implicit none
+         class(symba_pl),              intent(inout) :: self   !! SyMBA massive body object
+         class(swiftest_nbody_system), intent(inout) :: system !! Swiftest nbody system object
+         class(swiftest_parameters),   intent(in)    :: param  !! Current run configuration parameters
+      end subroutine symba_util_rearray_pl
    end interface
 
    interface util_resize

src/symba/symba_collision.f90

@@ -208,8 +208,8 @@ function symba_collision_consolidate_familes(pl, param, idx_parent, family, x, v
       fam_size = 2 + sum(nchild(:))
       allocate(family(fam_size))
       family = [parent_child_index_array(1)%idx(:),parent_child_index_array(2)%idx(:)]
-      fam_size = count(pl%status(family(:)) == COLLISION)
-      family = pack(family(:), pl%status(family(:)) == COLLISION)
+      fam_size = count(pl%lcollision(family(:)))
+      family = pack(family(:), pl%lcollision(family(:)))
       L_spin(:,:) = 0.0_DP
       Ip(:,:) = 0.0_DP
 
@@ -226,8 +226,8 @@ function symba_collision_consolidate_familes(pl, param, idx_parent, family, x, v
          if (nchild(j) > 0) then
             do i = 1, nchild(j) ! Loop over all children and take the mass weighted mean of the properties
                idx_child = parent_child_index_array(j)%idx(i + 1)
-               if ((idx_child) /= COLLISION) cycle
-               mchild = pl%Gmass(idx_child)
+               if (.not. pl%lcollision(idx_child)) cycle
+               mchild = pl%mass(idx_child)
                xchild(:) = pl%xb(:, idx_child)
                vchild(:) = pl%vb(:, idx_child)
                volchild = (4.0_DP / 3.0_DP) * PI * pl%radius(idx_child)**3
@@ -266,6 +266,7 @@ function symba_collision_consolidate_familes(pl, param, idx_parent, family, x, v
       return
    end function symba_collision_consolidate_familes
 
+
    module subroutine symba_collision_encounter_scrub(self, system, param)
       !! author: David A. Minton
       !! 

src/symba/symba_discard.f90

@@ -58,6 +58,108 @@ subroutine symba_discard_cb_pl(pl, system, param)
    end subroutine symba_discard_cb_pl
 
 
+   subroutine symba_discard_conserve_mtm(pl, system, param, ipl, lescape_body)
+      !! author: David A. Minton
+      !! 
+      !! Conserves system momentum when a body is lost from the system or collides with central body
+      implicit none
+      ! Arguments
+      class(symba_pl),           intent(inout) :: pl
+      class(symba_nbody_system), intent(inout) :: system
+      class(symba_parameters),   intent(inout) :: param
+      integer(I4B),              intent(in)    :: ipl
+      logical,                   intent(in)         :: lescape_body
+      ! Internals
+      real(DP), dimension(NDIM) :: Lpl, Ltot, Lcb, xcom, vcom
+      real(DP)                  :: pe, ke_orbit, ke_spin
+      integer(I4B)              :: i, oldstat
+
+      select type(cb => system%cb)
+      class is (symba_cb)
+
+         ! Add the potential and kinetic energy of the lost body to the records
+         pe = -cb%mass * pl%mass(ipl) / norm2(pl%xb(:, ipl) - pl%xb(:, 1))
+         ke_orbit = 0.5_DP * pl%mass(ipl) * dot_product(pl%vb(:, ipl), pl%vb(:, ipl)) 
+         if (param%lrotation) then
+            ke_spin  = 0.5_DP * pl%mass(ipl) * pl%radius(ipl)**2 * pl%Ip(3, ipl) * dot_product(pl%rot(:, ipl), pl%rot(:, ipl))
+         else
+            ke_spin = 0.0_DP
+         end if
+
+         ! Add the pre-collision ke of the central body to the records
+         ! Add planet mass to central body accumulator
+         if (lescape_body) then
+            system%Mescape = system%Mescape + pl%mass(ipl)
+            do i = 1, npl
+               if (i == ipl) cycle
+               pe = pe - pl%mass(i) * pl%mass(ipl) / norm2(xb(:, ipl) - xb(:, i))
+            end do
+
+            Ltot(:) = 0.0_DP
+            do i = 1, npl
+               Lpl(:) = mass(i) * pl%xb(:,i) .cross. pl%vb(:, i)
+               Ltot(:) = Ltot(:) + Lpl(:)
+            end do
+            Ltot(:) = Ltot(:) + cb%mass * cb%xb(:) .cross. cb%vb(:)
+            call pl%b2h(cb)
+            oldstat = status(ipl)
+            pl%status(ipl) = INACTIVE
+            call pl%h2b(cb)
+            pl%status(ipl) = oldstat
+            do i = 1, npl
+               if (i == ipl) cycle
+               Lpl(:) = mass(i) * pl%xb(:,i) .cross. pl%vb(:, i)
+               Ltot(:) = Ltot(:) - Lpl(:) 
+            end do 
+            Ltot(:) = Ltot(:) - cb%mass * cb%xb(:) .cross. cb%vb(:)
+            system%Lescape(:) = system%Lescape(:) + system%Ltot(:) 
+            if (param%lrotation) system%Lescape(:) = system%Lescape + pl%mass(ipl) * pl%radius(ipl)**2 * pl%Ip(3, ipl) * pl%rot(:, ipl)
+
+         else
+            xcom(:) = (pl%mass(ipl) * pl%xb(:, ipl) + cb%mass * cb%xb(:)) / (cb%mass + pl%mass(ipl))
+            vcom(:) = (pl%mass(ipl) * pl%vb(:, ipl) + cb%mass * cb%vb(:)) / (cb%mass + pl%mass(ipl))
+            Lpl(:) = (pl%xb(:,ipl) - xcom(:)) .cross. pL%vb(:,ipl) - vcom(:)
+            if (param%lrotation) Lpl(:) = pl%mass(ipl) * (Lpl(:) + pl%radius(ipl)**2 * pl%Ip(3,ipl) * pl%rot(:, ipl))
+
+            Lcb(:) = cb%mass * (cb%xb(:) - xcom(:)) .cross. (cb%vb(:) - vcom(:))
+
+            ke_orbit = ke_orbit + 0.5_DP * cb%mass * dot_product(cb%vb(:), cb%vb(:)) 
+            if (param%lrotation) ke_spin = ke_spin + 0.5_DP * cb%mass * cb%radius**2 * cb%Ip(3) * dot_product(cb%rot(:), cb%rot(:))
+            ! Update mass of central body to be consistent with its total mass
+            cb%dM = cb%dM + pl%mass(ipl)
+            cb%dR = cb%dR + 1.0_DP / 3.0_DP * (pl%radius(ipl) / cb%radius)**3 - 2.0_DP / 9.0_DP * (pl%radius(ipl) / cb%radius)**6
+            cb%mass = cb%M0 + cb%dM
+            cb%Gmass = param%GU * cb%mass 
+            cb%radius = cb%R0 + cb%dR
+            param%rmin = cb%radius
+            ! Add planet angular momentum to central body accumulator
+            cb%dL(:) = Lpl(:) + Lcb(:) + cb%dL(:)
+            ! Update rotation of central body to by consistent with its angular momentum 
+            if (param%lrotation) then
+               cb%rot(:) = (cb%L0(:) + cb%dL(:)) / (cb%Ip(3) * cb%mass * cb%radius**2)        
+               ke_spin  = ke_spin - 0.5_DP * cb%mass * cb%radius**2 * cb%Ip(3) * dot_product(cb%rot(:), cb%rot(:)) 
+            end if
+            cb%xb(:) = xcom(:)
+            cb%vb(:) = vcom(:)
+            ke_orbit = ke_orbit - 0.5_DP * cb%mass * dot_product(cb%vb(:), cb%vb(:)) 
+         end if
+         call pl%b2h(cb)
+
+         ! We must do this for proper book-keeping, since we can no longer track this body's contribution to energy directly
+         if (lescape_body) then
+            system%Ecollisions  = system%Ecollisions + ke_orbit + ke_spin + pe
+            system%Euntracked  = system%Euntracked - (ke_orbit + ke_spin + pe)
+         else
+            system%Ecollisions  = system%Ecollisions + pe 
+            system%Euntracked = system%Euntracked - pe
+         end if
+
+      end select
+      return
+
+   end subroutine symba_discard_conserve_mtm
+
+
    subroutine symba_discard_nonplpl(pl, system, param)
       !! author: David A. Minton
       !!
@@ -89,6 +191,45 @@ subroutine symba_discard_nonplpl(pl, system, param)
    end subroutine symba_discard_nonplpl
 
 
+   subroutine symba_discard_nonplpl_conservation(pl, system, param)
+      !! author: David A. Minton
+      !!
+      !! If there are any bodies that are removed due to either colliding with the central body or escaping the systme,
+      !! we need to track the conserved quantities with the system bookkeeping terms.
+      implicit none
+      ! Arguments
+      class(symba_pl),              intent(inout) :: pl     !! SyMBA test particle object
+      class(swiftest_nbody_system), intent(inout) :: system !! Swiftest nbody system object
+      class(swiftest_parameters),   intent(in)    :: param  !! Current run configuration parameters 
+      ! Internals
+      integer(I4B)                            :: i, ndiscard, dstat
+      logical                                 :: lescape
+      logical, dimension(pl%nbody)            :: discard_l_pl
+      integer(I4B), dimension(:), allocatable :: discard_index_list
+
+      associate(npl => pl%nbody)
+         discard_l_pl(1:npl) = pl%ldiscard(1:npl) .and. .not. pl%lcollision(1:npl) ! These are bodies that are discarded but not flagged as pl-pl collision
+         ndiscard = count(discard_l_pl(:)) 
+         allocate(discard_index_list(ndiscard))
+         discard_index_list(:) = pack([(i, i = 1, npl)], discard_l_pl(1:npl))
+         do i = 1, ndiscard
+            dstat = pl%status(discard_index_list(i)) 
+            if ((dstat == DISCARDED_RMIN) .or. (dstat == DISCARDED_PERI)) then
+               lescape = .false.
+            else if ((dstat == DISCARDED_RMAX) .or. (dstat == DISCARDED_RMAXU)) then
+               lescape = .true.
+            else 
+               cycle
+            end if
+            ! Conserve all the quantities
+            call symba_discard_conserve_mtm(pl, system, param, discard_index_list(i), lescape)
+         end do
+      end associate
+
+      return
+   end subroutine symba_discard_nonplpl_conservation
+
+
    subroutine symba_discard_peri_pl(pl, system, param)
       !! author: David A. Minton
       !!
@@ -150,17 +291,28 @@ module subroutine symba_discard_pl(self, system, param)
          select type(param)
          class is (symba_parameters)
             associate(pl => self, plplenc_list => system%plplenc_list)
+               call pl%h2b(system%cb) 
+
+               ! First deal with the non pl-pl collisions
                call symba_discard_nonplpl(self, system, param)
+
+               ! Scrub the pl-pl encounter list of any encounters that did not lead to a collision
                call plplenc_list%scrub_non_collision(system, param)
-               if (plplenc_list%nenc == 0) return ! No collisions to resolve
-               write(*, *) "Collision detected at time t = ",param%t
 
-               call pl%h2b(system%cb) 
-               if (param%lfragmentation) then
-                  call plplenc_list%resolve_fragmentations(system, param)
-               else
-                  call plplenc_list%resolve_mergers(system, param)
+               if ((plplenc_list%nenc > 0) .and. any(pl%lcollision(:))) then 
+                  write(*, *) "Collision between massive bodies detected at time t = ",param%t
+                  if (param%lfragmentation) then
+                     call plplenc_list%resolve_fragmentations(system, param)
+                  else
+                     call plplenc_list%resolve_mergers(system, param)
+                  end if
                end if
+
+               if (any(pl%ldiscard(:))) then
+                  call symba_discard_nonplpl_conservation(self, system, param)
+                  call pl%rearray(self, system, param)
+               end if
+
             end associate
          end select 
       end select

src/symba/symba_fragmentation.f90

@@ -73,7 +73,7 @@ module function symba_fragmentation_casemerge(system, param, family, x, v, mass,
             pe = 0.0_DP
             do j = 1, nfamily
                do i = j + 1, nfamily
-                  pe = pe - pl%Gmass(i) * pl%Gmass(j) / norm2(pl%xb(:, i) - pl%xb(:, j))
+                  pe = pe - pl%mass(i) * pl%mass(j) / norm2(pl%xb(:, i) - pl%xb(:, j))
                end do
             end do
             system%Ecollisions  = system%Ecollisions + pe 

src/symba/symba_util.f90

@@ -369,6 +369,23 @@ module subroutine symba_util_peri_pl(self, system, param)
    end subroutine symba_util_peri_pl
 
 
+   module subroutine symba_util_rearray_pl(self, system, param)
+      !! Author: the Purdue Swiftest Team -  David A. Minton, Carlisle A. Wishard, Jennifer L.L. Pouplin, and Jacob R. Elliott
+      !!
+      !! Clean up the massive body structures to remove discarded bodies and add new bodies
+      implicit none
+      ! Arguments
+      class(symba_pl),              intent(inout) :: self   !! SyMBA massive body object
+      class(swiftest_nbody_system), intent(inout) :: system !! Swiftest nbody system object
+      class(swiftest_parameters),   intent(in)    :: param  !! Current run configuration parameters
+
+      ! First 
+
+      return
+   end subroutine symba_util_rearray_pl
+
+
+
    module subroutine symba_util_resize_arr_info(arr, nnew)
       !! author: David A. Minton
       !!