Restructured discards in SyMBA to separate types of discards from eac…

…h other. Also added particle number checks to coordinate conversion methods in order to simplify calls to them

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
 

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
       !!

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

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)