diff --git a/src/discard/discard.f90 b/src/discard/discard.f90 index 5f87ee752..f4244145a 100644 --- a/src/discard/discard.f90 +++ b/src/discard/discard.f90 @@ -38,6 +38,7 @@ module subroutine discard_pl(self, system, param) class(swiftest_nbody_system), intent(inout) :: system !! Swiftest nbody system object class(swiftest_parameters), intent(in) :: param !! Current run configuration parameter + if (self%nbody == 0) return self%ldiscard(:) = .false. return @@ -58,17 +59,17 @@ module subroutine discard_tp(self, system, param) class(swiftest_parameters), intent(in) :: param !! Current run configuration parameter associate(tp => self, ntp => self%nbody, cb => system%cb, pl => system%pl, npl => system%pl%nbody) - if (ntp == 0) return + if ((ntp == 0) .or. (npl ==0)) return + if ((param%rmin >= 0.0_DP) .or. (param%rmax >= 0.0_DP) .or. & (param%rmaxu >= 0.0_DP) .or. ((param%qmin >= 0.0_DP) .and. (param%qmin_coord == "BARY"))) then - if (npl > 0) call pl%h2b(cb) - if (ntp > 0) call tp%h2b(cb) - end if - if ((param%rmin >= 0.0_DP) .or. (param%rmax >= 0.0_DP) .or. (param%rmaxu >= 0.0_DP)) then - if (ntp > 0) call discard_cb_tp(tp, system, param) + call pl%h2b(cb) + call tp%h2b(cb) end if - if (param%qmin >= 0.0_DP .and. ntp > 0) call discard_peri_tp(tp, system, param) - if (param%lclose .and. ntp > 0) call discard_pl_tp(tp, system, param) + + if ((param%rmin >= 0.0_DP) .or. (param%rmax >= 0.0_DP) .or. (param%rmaxu >= 0.0_DP)) call discard_cb_tp(tp, system, param) + if (param%qmin >= 0.0_DP) call discard_peri_tp(tp, system, param) + if (param%lclose) call discard_pl_tp(tp, system, param) if (any(tp%ldiscard)) call tp%spill(system%tp_discards, tp%ldiscard, ldestructive=.true.) end associate diff --git a/src/symba/symba_collision.f90 b/src/symba/symba_collision.f90 index bb67508fd..66aa1ef30 100644 --- a/src/symba/symba_collision.f90 +++ b/src/symba/symba_collision.f90 @@ -142,6 +142,7 @@ pure elemental function symba_collision_check_one(xr, yr, zr, vxr, vyr, vzr, Gmt return end function symba_collision_check_one + module subroutine symba_collision_make_family_pl(self, idx) !! author: Jennifer L.L. Pouplin, Carlisle A. wishard, and David A. Minton !! diff --git a/src/symba/symba_discard.f90 b/src/symba/symba_discard.f90 index f9066fdc6..7d15d996b 100644 --- a/src/symba/symba_discard.f90 +++ b/src/symba/symba_discard.f90 @@ -3,14 +3,45 @@ contains module subroutine symba_discard_pl(self, system, param) + !! author: David A. Minton + !! + !! Call the various flavors of discards for massive bodies in SyMBA runs, including discards due to colling with the central body, + !! escaping the system, or colliding with each other. implicit none ! Arguments class(symba_pl), intent(inout) :: self !! SyMBA test particle object class(swiftest_nbody_system), intent(inout) :: system !! Swiftest nbody system object class(swiftest_parameters), intent(in) :: param !! Current run configuration parameters + call symba_discard_nonplpl(self, system, param) + + select type(param) + class is (symba_parameters) + if (param%lfragmentation) then + + end if + + end select + + return + end subroutine symba_discard_pl + + subroutine symba_discard_nonplpl(pl, system, param) + !! author: David A. Minton + !! + !! Check to see if planets should be discarded based on their positions or because they are unbound + !s + !! + !! Adapted from David E. Kaufmann's Swifter routine: symba_discard_pl.f90 + !! Adapted from Hal Levison's Swift routine discard_massive5.f + implicit none + ! Arguments + class(symba_pl), intent(inout) :: pl !! SyMBA test particle object + class(swiftest_nbody_system), intent(inout) :: system !! Swiftest nbody system object + class(swiftest_parameters), intent(in) :: param !! Current run configuration parameters + ! First check for collisions with the central body - associate(pl => self, npl => self%nbody, cb => system%cb) + associate(npl => pl%nbody, cb => system%cb) if (npl == 0) return if ((param%rmin >= 0.0_DP) .or. (param%rmax >= 0.0_DP) .or. & (param%rmaxu >= 0.0_DP) .or. ((param%qmin >= 0.0_DP) .and. (param%qmin_coord == "BARY"))) then @@ -22,16 +53,9 @@ module subroutine symba_discard_pl(self, system, param) if (param%qmin >= 0.0_DP .and. npl > 0) call symba_discard_peri_pl(pl, system, param) end associate - select type(param) - class is (symba_parameters) - if (param%lfragmentation) then - - end if - - end select - return - end subroutine symba_discard_pl + end subroutine symba_discard_nonplpl + subroutine symba_discard_cb_pl(pl, system, param) @@ -42,7 +66,7 @@ subroutine symba_discard_cb_pl(pl, system, param) !! its collisional status will be revoked. Discards due to colliding with or escaping the central body take precedence !! over pl-pl collisions !! - !! Adapted from David E. Kaufmann's Swifter routine: symba_discard_cb.f90 + !! Adapted from David E. Kaufmann's Swifter routine: symba_discard_sun.f90 !! Adapted from Hal Levison's Swift routine discard_massive5.f implicit none ! Arguments diff --git a/src/util/util_coord.f90 b/src/util/util_coord.f90 index 938d04951..ab5fe85df 100644 --- a/src/util/util_coord.f90 +++ b/src/util/util_coord.f90 @@ -18,6 +18,7 @@ module subroutine util_coord_h2b_pl(self, cb) real(DP) :: msys real(DP), dimension(NDIM) :: xtmp, vtmp + if (self%nbody == 0) return associate(pl => self, npl => self%nbody) msys = cb%Gmass xtmp(:) = 0.0_DP @@ -51,6 +52,7 @@ module subroutine util_coord_h2b_tp(self, cb) class(swiftest_tp), intent(inout) :: self !! Swiftest test particle object class(swiftest_cb), intent(in) :: cb !! Swiftest central body object + if (self%nbody == 0) return associate(ntp => self%nbody, xbcb => cb%xb, vbcb => cb%vb, status => self%status, & xb => self%xb, xh => self%xh, vb => self%vb, vh => self%vh) @@ -83,6 +85,8 @@ module subroutine util_coord_b2h_pl(self, cb) ! Internals integer(I4B) :: i + if (self%nbody == 0) return + associate(npl => self%nbody, xbcb => cb%xb, vbcb => cb%vb, xb => self%xb, xh => self%xh, & vb => self%vb, vh => self%vh) do i = 1, NDIM @@ -107,6 +111,8 @@ module subroutine util_coord_b2h_tp(self, cb) class(swiftest_tp), intent(inout) :: self !! Swiftest massive body object class(swiftest_cb), intent(in) :: cb !! Swiftest central body object + if (self%nbody == 0) return + associate(ntp => self%nbody, xbcb => cb%xb, vbcb => cb%vb, xb => self%xb, xh => self%xh, & vb => self%vb, vh => self%vh, status => self%status) where(status(1:ntp) /= INACTIVE)