Added missing G term to potential energy calculations

MintonGroup · Aug 6, 2021 · 6706e16 · 6706e16
1 parent a9b56b3
commit 6706e16
Show file tree

Hide file tree

Showing 2 changed files with 22 additions and 22 deletions.
diff --git a/src/fragmentation/fragmentation.f90 b/src/fragmentation/fragmentation.f90
@@ -525,21 +525,21 @@ subroutine set_fragment_tan_vel(lerr)
 
          call calculate_system_energy(linclude_fragments=.true.)
          ke_frag_budget = -dEtot - Qloss
-         !write(*,*) '***************************************************'
-         !write(*,*) 'Original dis   : ',norm2(x(:,2) - x(:,1))
-         !write(*,*) 'r_max          : ',r_max
-         !write(*,*) 'f_spin         : ',f_spin
-         !write(*,*) '***************************************************'
-         !write(*,*) 'Energy balance so far: '
-         !write(*,*) 'ke_frag_budget : ',ke_frag_budget
-         !write(*,*) 'ke_orbit_before: ',ke_orbit_before 
-         !write(*,*) 'ke_orbit_after : ',ke_orbit_after  
-         !write(*,*) 'ke_spin_before : ',ke_spin_before 
-         !write(*,*) 'ke_spin_after  : ',ke_spin_after  
-         !write(*,*) 'pe_before      : ',pe_before 
-         !write(*,*) 'pe_after       : ',pe_after  
-         !write(*,*) 'Qloss          : ',Qloss
-         !write(*,*) '***************************************************'
+         write(*,*) '***************************************************'
+         write(*,*) 'Original dis   : ',norm2(x(:,2) - x(:,1))
+         write(*,*) 'r_max          : ',r_max
+         write(*,*) 'f_spin         : ',f_spin
+         write(*,*) '***************************************************'
+         write(*,*) 'Energy balance so far: '
+         write(*,*) 'ke_frag_budget : ',ke_frag_budget
+         write(*,*) 'ke_orbit_before: ',ke_orbit_before 
+         write(*,*) 'ke_orbit_after : ',ke_orbit_after  
+         write(*,*) 'ke_spin_before : ',ke_spin_before 
+         write(*,*) 'ke_spin_after  : ',ke_spin_after  
+         write(*,*) 'pe_before      : ',pe_before 
+         write(*,*) 'pe_after       : ',pe_after  
+         write(*,*) 'Qloss          : ',Qloss
+         write(*,*) '***************************************************'
          if (ke_frag_budget < 0.0_DP) then
             write(*,*) 'Negative ke_frag_budget: ',ke_frag_budget
             r_max_start = r_max_start / 2 
@@ -604,11 +604,11 @@ subroutine set_fragment_tan_vel(lerr)
 
          ! If we are over the energy budget, flag this as a failure so we can try again
          lerr = (ke_radial < 0.0_DP)
-         !write(*,*) 'Tangential'
-         !write(*,*) 'ke_frag_budget: ',ke_frag_budget
-         !write(*,*) 'ke_frag_orbit : ',ke_frag_orbit
-         !write(*,*) 'ke_frag_spin  : ',ke_frag_spin
-         !write(*,*) 'ke_radial     : ',ke_radial
+         write(*,*) 'Tangential'
+         write(*,*) 'ke_frag_budget: ',ke_frag_budget
+         write(*,*) 'ke_frag_orbit : ',ke_frag_orbit
+         write(*,*) 'ke_frag_spin  : ',ke_frag_spin
+         write(*,*) 'ke_radial     : ',ke_radial
 
          return
 

diff --git a/src/util/util_get_energy_momentum.f90 b/src/util/util_get_energy_momentum.f90
@@ -76,7 +76,7 @@ module subroutine util_get_energy_momentum_system(self, param)
          !$omp simd 
          do i = 1, npl
             associate(px => pl%xh(1,i), py => pl%xh(2,i), pz => pl%xh(3,i))
-               pecb(i) = -cb%mass * pl%mass(i) / sqrt(px**2 + py**2 + pz**2)
+               pecb(i) = -param%GU * cb%mass * pl%mass(i) / sqrt(px**2 + py**2 + pz**2)
             end associate
          end do
 
@@ -100,7 +100,7 @@ module subroutine util_get_energy_momentum_system(self, param)
                irh(i) = 1.0_DP / norm2(pl%xh(:,i))
             end do
             call obl_pot(npl, cb%mass, pl%mass, cb%j2rp2, cb%j4rp4, pl%xh, irh, oblpot)
-            system%pe = system%pe + oblpot
+            system%pe = system%pe + param%GU * oblpot
          end if
 
          system%Lorbit(1) = sum(Lplorbitx(1:npl), lstatus(1:npl))