Added total energy tracker and update the collisional energy bookkeep…

…ing term after discards

src/symba/symba_discard.f90

@@ -271,7 +271,6 @@ subroutine symba_discard_peri_pl(pl, system, param)
       pl%lfirst = lfirst_orig
 
       return
-
    end subroutine symba_discard_peri_pl
 
 
@@ -285,6 +284,8 @@ module subroutine symba_discard_pl(self, system, param)
       class(symba_pl),              intent(inout) :: self   !! SyMBA test particle object
       class(swiftest_nbody_system), intent(inout) :: system !! Swiftest nbody system object
       class(swiftest_parameters),   intent(inout) :: param  !! Current run configuration parameters 
+      ! Internals
+      real(DP) :: Eorbit_before, Eorbit_after
 
       select type(system)
       class is (symba_nbody_system)
@@ -309,8 +310,17 @@ module subroutine symba_discard_pl(self, system, param)
                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 = Eorbit_after - Eorbit_before
+                  end if
                end if
 
             end associate

src/util/util_get_energy_momentum.f90

@@ -129,6 +129,8 @@ module subroutine util_get_energy_momentum_system(self, param)
             system%Lspin(2) = Lcbspin(2) + sum(Lplspiny(1:npl), lstatus(1:npl)) 
             system%Lspin(3) = Lcbspin(3) + sum(Lplspinz(1:npl), lstatus(1:npl)) 
          end if
+
+         system%te = system%ke_orbit + system%ke_spin + system%pe
       end associate
 
       return