Fixed bug where Gmass was used instead of mass to compute central bod…

…y spin angular momentum
MintonGroup · Jan 10, 2023 · 2247213 · 2247213
1 parent a34f851
commit 2247213
Show file tree

Hide file tree

Showing 2 changed files with 10 additions and 6 deletions.
diff --git a/src/swiftest/swiftest_io.f90 b/src/swiftest/swiftest_io.f90
@@ -567,7 +567,7 @@ module subroutine swiftest_io_netcdf_get_t0_values_system(self, param)
       real(DP), dimension(:), allocatable       :: vals
       real(DP), dimension(1)                    :: rtemp
       real(DP), dimension(NDIM)                 :: rot0, Ip0
-      real(DP) :: KE_orb_orig, KE_spin_orig, PE_orig, BE_orig
+      real(DP) :: KE_orb_orig, KE_spin_orig, PE_orig, BE_orig, mass0
 
       associate (nc => self%system_history%nc, cb => self%cb)
          call nc%open(param, readonly=.true.)
@@ -610,7 +610,8 @@ module subroutine swiftest_io_netcdf_get_t0_values_system(self, param)
             if (param%lrotation) then
                call netcdf_io_check( nf90_get_var(nc%id, nc%rot_varid, rot0, start=[1,1,tslot], count=[NDIM,1,1]), "netcdf_io_get_t0_values_system rot_varid" )
                call netcdf_io_check( nf90_get_var(nc%id, nc%Ip_varid, Ip0, start=[1,1,tslot], count=[NDIM,1,1]), "netcdf_io_get_t0_values_system Ip_varid" )
-               cb%L0(:) = Ip0(3) * cb%GM0 * cb%R0**2 * rot0(:)
+               mass0 = cb%GM0 / param%GU
+               cb%L0(:) = Ip0(3) * mass0 * cb%R0**2 * rot0(:)
             end if
 
             ! Retrieve the current bookkeeping variables
@@ -1119,7 +1120,7 @@ module function swiftest_io_netcdf_read_frame_system(self, nc, param) result(ier
             end do
 
             ! Set initial central body angular momentum for bookkeeping
-            cb%L0(:) = cb%Ip(3) * cb%GM0 * cb%R0**2 * cb%rot(:)         
+            cb%L0(:) = cb%Ip(3) * cb%mass * cb%R0**2 * cb%rot(:)         
          end if
 
          ! if (param%ltides) then

diff --git a/src/symba/symba_discard.f90 b/src/symba/symba_discard.f90
@@ -116,7 +116,7 @@ subroutine symba_discard_conserve_mtm(pl, nbody_system, param, ipl, lescape_body
       integer(I4B),              intent(in)    :: ipl
       logical,                   intent(in)         :: lescape_body
       ! Internals
-      real(DP), dimension(NDIM) :: Lpl, L_total, Lcb, xcom, vcom
+      real(DP), dimension(NDIM) :: Lpl, L_total, Lcb, xcom, vcom, drot0, drot1
       real(DP)                  :: pe, be, ke_orbit, ke_spin
       integer(I4B)              :: i, oldstat
 
@@ -167,7 +167,8 @@ subroutine symba_discard_conserve_mtm(pl, nbody_system, param, ipl, lescape_body
             xcom(:) = (pl%mass(ipl) * pl%rb(:, ipl) + cb%mass * cb%rb(:)) / (cb%mass + pl%mass(ipl))
             vcom(:) = (pl%mass(ipl) * pl%vb(:, ipl) + cb%mass * cb%vb(:)) / (cb%mass + pl%mass(ipl))
             Lpl(:) = (pl%rb(:,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))
+            if (param%lrotation) Lpl(:) = Lpl(:) + pl%radius(ipl)**2 * pl%Ip(3,ipl) * pl%rot(:, ipl)
+            Lpl(:) = pl%mass(ipl) * Lpl(:)
 
             Lcb(:) = cb%mass * ((cb%rb(:) - xcom(:)) .cross. (cb%vb(:) - vcom(:)))
 
@@ -184,7 +185,9 @@ subroutine symba_discard_conserve_mtm(pl, nbody_system, param, ipl, lescape_body
             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)        
+               drot0(:) = cb%L0(:)/ (cb%Ip(3) * cb%mass * cb%radius**2)  
+               drot1(:) = cb%dL(:) / (cb%Ip(3) * cb%mass * cb%radius**2)        
+               cb%rot(:) = drot0(:) + drot1(:)
                ke_spin  = ke_spin - 0.5_DP * cb%mass * cb%radius**2 * cb%Ip(3) * dot_product(cb%rot(:), cb%rot(:)) 
             end if
             cb%rb(:) = xcom(:)