From 5004848f6b489d7670e2a697adbda082a167ee2a Mon Sep 17 00:00:00 2001 From: David Minton Date: Mon, 24 May 2021 23:20:34 -0400 Subject: [PATCH] Fixed problem in which angular momentum would be miscalculated if a collision and sun/escape discard happened on the same step. The coordinate conversion subroutines were too strict with status flags (they just need to not be INACTIVE, rather than only being ACTIVE). --- src/coord/coord_h2b.f90 | 2 +- src/coord/coord_vb2vh.f90 | 2 +- src/coord/coord_vh2vb.f90 | 3 +-- src/symba/symba_discard_conserve_mtm.f90 | 2 -- src/symba/symba_discard_pl.f90 | 3 --- 5 files changed, 3 insertions(+), 9 deletions(-) diff --git a/src/coord/coord_h2b.f90 b/src/coord/coord_h2b.f90 index a5396842e..ba1569498 100644 --- a/src/coord/coord_h2b.f90 +++ b/src/coord/coord_h2b.f90 @@ -25,7 +25,7 @@ subroutine coord_h2b(npl, swiftest_plA, msys) mass => swiftest_plA%mass, status => swiftest_plA%status, & dMcb => swiftest_plA%dMcb, Mcb_initial => swiftest_plA%Mcb_initial) - lstatus(2:npl) = status(2:npl) == ACTIVE + lstatus(2:npl) = status(2:npl) /= INACTIVE xbcb(:) = 0.0_DP vbcb(:) = 0.0_DP diff --git a/src/coord/coord_vb2vh.f90 b/src/coord/coord_vb2vh.f90 index 84f4a586b..49ba9b1b9 100644 --- a/src/coord/coord_vb2vh.f90 +++ b/src/coord/coord_vb2vh.f90 @@ -23,7 +23,7 @@ subroutine coord_vb2vh(npl, swiftest_plA) vb => swiftest_plA%vb, vh => swiftest_plA%vh, & mass => swiftest_plA%mass, status => swiftest_plA%status) - lstatus(2:npl) = status(2:npl) == ACTIVE + lstatus(2:npl) = status(2:npl) /= INACTIVE vbcb(:) = 0.0_DP do i = 2, npl if (.not.lstatus(i)) cycle diff --git a/src/coord/coord_vh2vb.f90 b/src/coord/coord_vh2vb.f90 index b8db6cb63..57fa3ab24 100644 --- a/src/coord/coord_vh2vb.f90 +++ b/src/coord/coord_vh2vb.f90 @@ -26,7 +26,7 @@ subroutine coord_vh2vb(npl, swiftest_plA, msys) mass => swiftest_plA%mass, status => swiftest_plA%status, & dMcb => swiftest_plA%dMcb, Mcb_initial => swiftest_plA%Mcb_initial) - lstatus(2:npl) = status(2:npl) == ACTIVE + lstatus(2:npl) = status(2:npl) /= INACTIVE vbcb(:) = 0.0_DP do i = 2,npl @@ -37,7 +37,6 @@ subroutine coord_vh2vb(npl, swiftest_plA, msys) msys = dMcb + sum(mass(2:npl), lstatus(2:npl)) + Mcb_initial vbcb(:) = -vbcb(:) / msys - do i = 2,npl vb(:,i) = vh(:,i) + vbcb(:) end do diff --git a/src/symba/symba_discard_conserve_mtm.f90 b/src/symba/symba_discard_conserve_mtm.f90 index d4717fd0a..d19325607 100644 --- a/src/symba/symba_discard_conserve_mtm.f90 +++ b/src/symba/symba_discard_conserve_mtm.f90 @@ -77,8 +77,6 @@ subroutine symba_discard_conserve_mtm(param, swiftest_plA, ipl, lescape_body) Euntracked = Euntracked - pe end if - ! Update the heliocentric coordinates of everything else - call coord_b2h(npl, swiftest_plA) end associate return diff --git a/src/symba/symba_discard_pl.f90 b/src/symba/symba_discard_pl.f90 index 2607dfc07..d1df88b44 100644 --- a/src/symba/symba_discard_pl.f90 +++ b/src/symba/symba_discard_pl.f90 @@ -42,9 +42,6 @@ subroutine symba_discard_pl(t, npl, ntp, symba_plA, symba_tpA, rmin, rmax, rmaxu if (allocated(discard_stat_list)) deallocate(discard_stat_list) allocate(discard_stat_list(ndiscard)) discard_stat_list = pack(status(1:npl), discard_l_pl) - where(discard_l_pl(1:npl)) - status(:) = ACTIVE - end where end if end associate