From 2247213be45cdbc6119a563691a1e9ffe20dd14e Mon Sep 17 00:00:00 2001 From: David A Minton Date: Tue, 10 Jan 2023 15:01:54 -0500 Subject: [PATCH] Fixed bug where Gmass was used instead of mass to compute central body spin angular momentum --- src/swiftest/swiftest_io.f90 | 7 ++++--- src/symba/symba_discard.f90 | 9 ++++++--- 2 files changed, 10 insertions(+), 6 deletions(-) diff --git a/src/swiftest/swiftest_io.f90 b/src/swiftest/swiftest_io.f90 index 7e25b316f..7bf7c1117 100644 --- 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 index 88d674fb5..674772768 100644 --- 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(:)