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).

MintonGroup · May 25, 2021 · 5004848 · 5004848
1 parent 2e7df1d
commit 5004848
Show file tree

Hide file tree

Showing 5 changed files with 3 additions and 9 deletions.
diff --git 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
@@ -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
@@ -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
@@ -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
@@ -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