From 14d37ee8fca33b563ceae5559e104f98b171be7a Mon Sep 17 00:00:00 2001 From: David A Minton Date: Sun, 8 Aug 2021 17:15:50 -0400 Subject: [PATCH] Added missing central body terms to energy and momentum calculations --- src/util/util_get_energy_momentum.f90 | 40 +++++++++++++++++++++------ 1 file changed, 31 insertions(+), 9 deletions(-) diff --git a/src/util/util_get_energy_momentum.f90 b/src/util/util_get_energy_momentum.f90 index 5cf8d20c7..90f0d2242 100644 --- a/src/util/util_get_energy_momentum.f90 +++ b/src/util/util_get_energy_momentum.f90 @@ -16,6 +16,7 @@ module subroutine util_get_energy_momentum_system(self, param) integer(I4B) :: i, j integer(I8B) :: k real(DP) :: rmag, v2, rot2, oblpot, hx, hy, hz, hsx, hsy, hsz + real(DP) :: kecb, kespincb, Lcborbitx, Lcborbity, Lcborbitz, Lcbspinx, Lcbspiny, Lcbspinz real(DP), dimension(self%pl%nbody) :: irh, kepl, kespinpl, pecb real(DP), dimension(self%pl%nbody) :: Lplorbitx, Lplorbity, Lplorbitz real(DP), dimension(self%pl%nbody) :: Lplspinx, Lplspiny, Lplspinz @@ -38,6 +39,13 @@ module subroutine util_get_energy_momentum_system(self, param) Lplspinz(:) = 0.0_DP lstatus(1:npl) = pl%status(1:npl) /= INACTIVE call pl%h2b(cb) + kecb = cb%mass * dot_product(cb%vb(:), cb%vb(:)) + hx = cb%xb(2) * cb%vb(3) - cb%xb(3) * cb%vb(2) + hy = cb%xb(3) * cb%vb(1) - cb%xb(1) * cb%vb(3) + hz = cb%xb(1) * cb%vb(2) - cb%xb(2) * cb%vb(1) + Lcborbitx = cb%mass * hx + Lcborbity = cb%mass * hy + Lcborbitz = cb%mass * hz !!$omp simd private(v2, rot2, hx, hy, hz) do i = 1, npl v2 = dot_product(pl%vb(:,i), pl%vb(:,i)) @@ -55,6 +63,17 @@ module subroutine util_get_energy_momentum_system(self, param) end do if (param%lrotation) then + kespincb = cb%mass * cb%Ip(3) * cb%radius**2 * dot_product(cb%rot(:), cb%rot(:)) + ! For simplicity, we always assume that the rotation pole is the 3rd principal axis + hsx = cb%Ip(3) * cb%radius**2 * cb%rot(1) + hsy = cb%Ip(3) * cb%radius**2 * cb%rot(2) + hsz = cb%Ip(3) * cb%radius**2 * cb%rot(3) + + ! Angular momentum from spin + Lcbspinx = cb%mass * hsx + Lcbspiny = cb%mass * hsy + Lcbspinz = cb%mass * hsz + do i = 1, npl rot2 = dot_product(pl%rot(:,i), pl%rot(:,i)) ! For simplicity, we always assume that the rotation pole is the 3rd principal axis @@ -69,6 +88,7 @@ module subroutine util_get_energy_momentum_system(self, param) kespinpl(i) = pl%mass(i) * pl%Ip(3, i) * pl%radius(i)**2 * rot2 end do else + kespincb = 0.0_DP kespinpl(:) = 0.0_DP end if @@ -88,8 +108,8 @@ module subroutine util_get_energy_momentum_system(self, param) end associate end do - system%ke_orbit = 0.5_DP * sum(kepl(1:npl), lstatus(:)) - if (param%lrotation) system%ke_spin = 0.5_DP * sum(kespinpl(1:npl), lstatus(:)) + system%ke_orbit = 0.5_DP * (kecb + sum(kepl(1:npl), lstatus(:))) + if (param%lrotation) system%ke_spin = 0.5_DP * (kespincb + sum(kespinpl(1:npl), lstatus(:))) system%pe = sum(pepl(:), lstatpl(:)) + sum(pecb(1:npl), lstatus(1:npl)) @@ -103,13 +123,15 @@ module subroutine util_get_energy_momentum_system(self, param) system%pe = system%pe + oblpot end if - system%Lorbit(1) = sum(Lplorbitx(1:npl), lstatus(1:npl)) - system%Lorbit(2) = sum(Lplorbity(1:npl), lstatus(1:npl)) - system%Lorbit(3) = sum(Lplorbitz(1:npl), lstatus(1:npl)) - - system%Lspin(1) = sum(Lplspinx(1:npl), lstatus(1:npl)) - system%Lspin(2) = sum(Lplspiny(1:npl), lstatus(1:npl)) - system%Lspin(3) = sum(Lplspinz(1:npl), lstatus(1:npl)) + system%Lorbit(1) = Lcborbitx + sum(Lplorbitx(1:npl), lstatus(1:npl)) + system%Lorbit(2) = Lcborbity + sum(Lplorbity(1:npl), lstatus(1:npl)) + system%Lorbit(3) = Lcborbitz + sum(Lplorbitz(1:npl), lstatus(1:npl)) + + if (param%lrotation) then + system%Lspin(1) = Lcbspinx + sum(Lplspinx(1:npl), lstatus(1:npl)) + system%Lspin(2) = Lcbspiny + sum(Lplspiny(1:npl), lstatus(1:npl)) + system%Lspin(3) = Lcbspinz + sum(Lplspinz(1:npl), lstatus(1:npl)) + end if end associate return