Fixed a bunch of issues to get energy values computed in mergers

src/fraggle/fraggle_set.f90

@@ -287,7 +287,6 @@ module subroutine fraggle_set_original_scale_factors(self, colliders)
          do i = 1, 2
             colliders%rot(:,i) = colliders%L_spin(:,i) * (colliders%mass(i) * colliders%radius(i)**2 * colliders%Ip(3, i))
          end do
-         frag%Qloss = frag%Qloss * frag%Escale
 
          frag%mtot = frag%mtot * frag%mscale
          frag%mass = frag%mass * frag%mscale
@@ -300,7 +299,9 @@ module subroutine fraggle_set_original_scale_factors(self, colliders)
             frag%rb(:, i) = frag%x_coll(:, i) + frag%rbcom(:)
             frag%vb(:, i) = frag%v_coll(:, i) + frag%vbcom(:)
          end do
-
+
+         frag%Qloss = frag%Qloss * frag%Escale
+
          frag%Lorbit_before(:) = frag%Lorbit_before * frag%Lscale
          frag%Lspin_before(:) = frag%Lspin_before * frag%Lscale
          frag%Ltot_before(:) = frag%Ltot_before * frag%Lscale

src/modules/fraggle_classes.f90

@@ -86,12 +86,12 @@ module fraggle_classes
       real(DP)                  :: Etot_before, Etot_after         !! Before/after total system energy
 
       ! Scale factors used to scale dimensioned quantities to a more "natural" system where important quantities (like kinetic energy, momentum) are of order ~1
-      real(DP) :: dscale !! Distance dimension scale factor
-      real(DP) :: mscale !! Mass scale factor
-      real(DP) :: tscale !! Time scale factor
-      real(DP) :: vscale !! Velocity scale factor (a convenience unit that is derived from dscale and tscale)
-      real(DP) :: Escale !! Energy scale factor (a convenience unit that is derived from dscale, tscale, and mscale)
-      real(DP) :: Lscale !! Angular momentum scale factor (a convenience unit that is derived from dscale, tscale, and mscale)
+      real(DP) :: dscale = 1.0_DP !! Distance dimension scale factor
+      real(DP) :: mscale = 1.0_DP !! Mass scale factor
+      real(DP) :: tscale = 1.0_DP !! Time scale factor
+      real(DP) :: vscale = 1.0_DP !! Velocity scale factor (a convenience unit that is derived from dscale and tscale)
+      real(DP) :: Escale = 1.0_DP !! Energy scale factor (a convenience unit that is derived from dscale, tscale, and mscale)
+      real(DP) :: Lscale = 1.0_DP  !! Angular momentum scale factor (a convenience unit that is derived from dscale, tscale, and mscale)
    contains
       procedure :: generate_fragments      => fraggle_generate_fragments         !! Generates a system of fragments in barycentric coordinates that conserves energy and momentum.
       procedure :: accel                   => fraggle_placeholder_accel          !! Placeholder subroutine to fulfill requirement for an accel method

src/symba/symba_collision.f90

@@ -185,6 +185,10 @@ module function symba_collision_casemerge(system, param, t)  result(status)
          class is (symba_pl)
 
             call fragments%set_mass_dist(colliders, param)
+
+            ! Calculate the initial energy of the system without the collisional family
+            call fragments%get_energy_and_momentum(colliders, system, param, lbefore=.true.)
+
             ibiggest = colliders%idx(maxloc(pl%Gmass(colliders%idx(:)), dim=1))
             fragments%id(1) = pl%id(ibiggest)
             fragments%rb(:,1) = fragments%rbcom(:)
@@ -197,18 +201,16 @@ module function symba_collision_casemerge(system, param, t)  result(status)
                ! Assume prinicpal axis rotation on 3rd Ip axis
                fragments%rot(:,1) = L_spin_new(:) / (fragments%Ip(3,1) * fragments%mass(1) * fragments%radius(1)**2)
             else ! If spin is not enabled, we will consider the lost pre-collision angular momentum as "escaped" and add it to our bookkeeping variable
-               param%Lescape(:) = param%Lescape(:) + colliders%L_orbit(:,1) + colliders%L_orbit(:,2) 
+               system%Lescape(:) = system%Lescape(:) + colliders%L_orbit(:,1) + colliders%L_orbit(:,2) 
             end if
 
             ! Keep track of the component of potential energy due to the pre-impact colliders%idx for book-keeping
-            pe = 0.0_DP
-            do j = 1, colliders%ncoll
-               do i = j + 1, colliders%ncoll
-                  pe = pe - pl%Gmass(i) * pl%mass(j) / norm2(pl%rb(:, i) - pl%rb(:, j))
-               end do
-            end do
-            system%Ecollisions = system%Ecollisions + pe 
-            system%Euntracked  = system%Euntracked - pe 
+            ! Get the energy of the system after the collision
+            call fragments%get_energy_and_momentum(colliders, system, param, lbefore=.false.)
+            pe = fragments%pe_after - fragments%pe_before 
+            system%Ecollisions = system%Ecollisions - pe 
+            system%Euntracked  = system%Euntracked + pe 
+
 
             ! Update any encounter lists that have the removed bodies in them so that they instead point to the new 
             do k = 1, system%plplenc_list%nenc