Refactored msys to system%Gmtot. Added back the massive body pericent…

…er check (originally symba_peri.f90 and symba_discard_peri_pl.f90 from Swifter. Now symba_util_peri_pl bound to the get_peri method call, and symba_discard_peri_pl

src/discard/discard.f90

@@ -94,7 +94,7 @@ subroutine discard_cb_tp(tp, system, param)
       integer(I4B)        :: i
       real(DP)            :: energy, vb2, rb2, rh2, rmin2, rmax2, rmaxu2
 
-      associate(ntp => tp%nbody, cb => system%cb, t => param%t, msys => system%msys)
+      associate(ntp => tp%nbody, cb => system%cb, t => param%t, Gmtot => system%Gmtot)
          rmin2 = max(param%rmin * param%rmin, cb%radius * cb%radius)
          rmax2 = param%rmax**2
          rmaxu2 = param%rmaxu**2
@@ -114,7 +114,7 @@ subroutine discard_cb_tp(tp, system, param)
                else if (param%rmaxu >= 0.0_DP) then
                   rb2 = dot_product(tp%xb(:, i),  tp%xb(:, i))
                   vb2 = dot_product(tp%vb(:, i), tp%vb(:, i))
-                  energy = 0.5_DP * vb2 - msys / sqrt(rb2)
+                  energy = 0.5_DP * vb2 - Gmtot / sqrt(rb2)
                   if ((energy > 0.0_DP) .and. (rb2 > rmaxu2)) then
                      tp%status(i) = DISCARDED_RMAXU
                      write(*, *) "Particle ", tp%id(i), " is unbound and too far from barycenter at t = ", t

src/helio/helio_coord.f90

@@ -68,14 +68,14 @@ module subroutine helio_coord_vh2vb_pl(self, cb)
       class(swiftest_cb),     intent(inout) :: cb   !! Swiftest central body object
       ! Internals
       integer(I4B)  :: i
-      real(DP)      :: msys
+      real(DP)      :: Gmtot
 
       if (self%nbody == 0) return
 
       associate(pl => self, npl => self%nbody)
-         msys = cb%Gmass + sum(pl%Gmass(1:npl))
+         Gmtot = cb%Gmass + sum(pl%Gmass(1:npl))
          do i = 1, NDIM
-            cb%vb(i) = -sum(pl%Gmass(1:npl) * pl%vh(i, 1:npl)) / msys
+            cb%vb(i) = -sum(pl%Gmass(1:npl) * pl%vh(i, 1:npl)) / Gmtot
             pl%vb(i, 1:npl) = pl%vh(i, 1:npl) + cb%vb(i)
          end do
       end associate

src/modules/swiftest_classes.f90

@@ -259,7 +259,7 @@ module swiftest_classes
       class(swiftest_pl),            allocatable :: pl                   !! Massive body data structure
       class(swiftest_tp),            allocatable :: tp                   !! Test particle data structure
       class(swiftest_tp),            allocatable :: tp_discards          !! Discarded test particle data structure
-      real(DP)                                   :: msys = 0.0_DP        !! Total system mass - used for barycentric coordinate conversion
+      real(DP)                                   :: Gmtot = 0.0_DP        !! Total system mass - used for barycentric coordinate conversion
       real(DP)                                   :: ke = 0.0_DP          !! System kinetic energy
       real(DP)                                   :: pe = 0.0_DP          !! System potential energy
       real(DP)                                   :: te = 0.0_DP          !! System total energy

src/modules/symba_classes.f90

@@ -94,6 +94,7 @@ module symba_classes
       procedure :: setup           => symba_setup_pl                 !! Constructor method - Allocates space for number of particle
       procedure :: append          => symba_util_append_pl           !! Appends elements from one structure to another
       procedure :: fill            => symba_util_fill_pl             !! "Fills" bodies from one object into another depending on the results of a mask (uses the UNPACK intrinsic)
+      procedure :: get_peri        => symba_util_peri_pl             !! Determine system pericenter passages for massive bodies
       procedure :: resize          => symba_util_resize_pl           !! Checks the current size of a Swiftest body against the requested size and resizes it if it is too small.
       procedure :: sort            => symba_util_sort_pl             !! Sorts body arrays by a sortable componen
       procedure :: rearrange       => symba_util_sort_rearrange_pl   !! Rearranges the order of array elements of body based on an input index array. Used in sorting methods
@@ -479,6 +480,14 @@ module subroutine symba_util_fill_tp(self, inserts, lfill_list)
          class(swiftest_body),  intent(in)    :: inserts    !! Inserted object 
          logical, dimension(:), intent(in)    :: lfill_list !! Logical array of bodies to merge into the keeps
       end subroutine symba_util_fill_tp
+
+      module subroutine symba_util_peri_pl(self, system, param)
+         use swiftest_classes, only : swiftest_nbody_system, swiftest_parameters
+         implicit none
+         class(symba_pl),              intent(inout) :: self   !! SyMBA massive body object
+         class(swiftest_nbody_system), intent(inout) :: system !! Swiftest nbody system object
+         class(swiftest_parameters),   intent(in)    :: param  !! Current run configuration parameters
+      end subroutine symba_util_peri_pl
    end interface
 
    interface util_resize

src/symba/symba_discard.f90

@@ -2,62 +2,6 @@
    use swiftest
 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(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
-            call pl%h2b(cb) 
-         end if
-         if ((param%rmin >= 0.0_DP) .or. (param%rmax >= 0.0_DP) .or.  (param%rmaxu >= 0.0_DP)) then
-            call symba_discard_cb_pl(pl, system, param)
-         end if
-         if (param%qmin >= 0.0_DP .and. npl > 0) call symba_discard_peri_pl(pl, system, param)
-      end associate
-
-      return
-   end subroutine symba_discard_nonplpl
-
-
-
    subroutine symba_discard_cb_pl(pl, system, param)
       !! author: David A. Minton
       !!
@@ -98,7 +42,7 @@ subroutine symba_discard_cb_pl(pl, system, param)
                else if (param%rmaxu >= 0.0_DP) then
                   rb2 = dot_product(pl%xb(:,i), pl%xb(:,i))
                   vb2 = dot_product(pl%vb(:,i), pl%vb(:,i))
-                  energy = 0.5_DP * vb2 - system%msys / sqrt(rb2)
+                  energy = 0.5_DP * vb2 - system%Gmtot / sqrt(rb2)
                   if ((energy > 0.0_DP) .and. (rb2 > rmaxu2)) then
                      pl%ldiscard(i) = .true.
                      pl%lcollision(i) = .false. 
@@ -114,19 +58,104 @@ subroutine symba_discard_cb_pl(pl, system, param)
    end subroutine symba_discard_cb_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(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
+            call pl%h2b(cb) 
+         end if
+         if ((param%rmin >= 0.0_DP) .or. (param%rmax >= 0.0_DP) .or.  (param%rmaxu >= 0.0_DP)) then
+            call symba_discard_cb_pl(pl, system, param)
+         end if
+         if (param%qmin >= 0.0_DP .and. npl > 0) call symba_discard_peri_pl(pl, system, param)
+      end associate
+
+      return
+   end subroutine symba_discard_nonplpl
+
+
    subroutine symba_discard_peri_pl(pl, system, param)
       !! author: David A. Minton
       !!
       !! Check to see if a test particle should be discarded because its perihelion distance becomes too small
       !!
-      !! Adapted from David E. Kaufmann's Swifter routine: discard_peri.f90
-      !! Adapted from Hal Levison's Swift routine discard_peri.f
+      !! Adapted from David E. Kaufmann's Swifter routine: symba_discard_peri_pl.f90
+      !! Adapted from Hal Levison's Swift routine discard_mass_peri.f
       implicit none
       ! Arguments
       class(symba_pl),              intent(inout) :: pl     !! SyMBA massive body object
       class(swiftest_nbody_system), intent(inout) :: system !! Swiftest nbody system object
       class(swiftest_parameters),   intent(in)    :: param  !! Current run configuration parameters
+      ! Internals
+      logical, save      :: lfirst = .true.
+      logical            :: lfirst_orig
+      integer(I4B)       :: i
+
+
+      lfirst_orig = pl%lfirst
+      pl%lfirst = lfirst
+      if (lfirst) then
+         call pl%get_peri(system, param)          
+         lfirst = .false.
+      else
+         call pl%get_peri(system, param)          
+         do i = 1, pl%nbody
+            if (pl%status(i) == ACTIVE) then
+               if ((pl%isperi(i) == 0) .and. (pl%nplenc(i)== 0)) then
+                  if ((pl%atp(i) >= param%qmin_alo) .and. (pl%atp(i) <= param%qmin_ahi) .and. (pl%peri(i) <= param%qmin)) then
+                     pl%ldiscard(i) = .true.
+                     pl%lcollision(i) = .false.
+                     pl%status(i) = DISCARDED_PERI
+                     write(*, *) "Particle ", pl%id(i), " perihelion distance too small at t = ", param%t
+                  end if
+               end if
+            end if
+         end do
+      end if
+      pl%lfirst = lfirst_orig
+
       return
+
    end subroutine symba_discard_peri_pl
 
+
+   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
+
 end submodule s_symba_discard

src/symba/symba_util.f90

@@ -271,6 +271,92 @@ module subroutine symba_util_fill_tp(self, inserts, lfill_list)
    end subroutine symba_util_fill_tp
 
 
+   module subroutine symba_util_peri_pl(self, system, param)
+      !! author: David A. Minton
+      !!
+      !! Determine system pericenter passages for planets in SyMBA
+      !!
+      !! Adapted from David E. Kaufmann's Swifter routine: symba_peri.f90
+      !! Adapted from Hal Levison's Swift routine util_mass_peri.f
+      implicit none
+      ! Arguments
+      class(symba_pl),              intent(inout) :: self   !! SyMBA massive body object
+      class(swiftest_nbody_system), intent(inout) :: system !! Swiftest nbody system object
+      class(swiftest_parameters),   intent(in)    :: param  !! Current run configuration parameters
+      ! Internals
+      integer(I4B)       :: i
+      real(DP)           :: vdotr, e
+
+      associate(pl => self, npl => self%nbody)
+         if (pl%lfirst) then
+            if (param%qmin_coord == "HELIO") then
+               do i = 1, npl
+                  if (pl%status(i) == ACTIVE) then
+                     vdotr = dot_product(pl%xh(:,i), pl%vh(:,i))
+                     if (vdotr > 0.0_DP) then
+                        pl%isperi(i) = 1
+                     else
+                        pl%isperi(i) = -1
+                     end if
+                  end if
+               end do
+            else
+               do i = 1, npl
+                  if (pl%status(i) == ACTIVE) then
+                     vdotr = dot_product(pl%xb(:,i), pl%vb(:,i))
+                     if (vdotr > 0.0_DP) then
+                        pl%isperi(i) = 1
+                     else
+                        pl%isperi(i) = -1
+                     end if
+                  end if
+               end do
+            end if
+         else
+            if (param%qmin_coord == "HELIO") then
+               do i = 1, npl
+                  if (pl%status(i) == ACTIVE) then
+                     vdotr = dot_product(pl%xh(:,i), pl%vh(:,i))
+                     if (pl%isperi(i) == -1) then
+                        if (vdotr >= 0.0_DP) then
+                           pl%isperi(i) = 0
+                           CALL orbel_xv2aeq(pl%mu(i), pl%xh(:,i), pl%vh(:,i), pl%atp(i), e, pl%peri(i))
+                        end if
+                     else
+                        if (vdotr > 0.0_DP) then
+                           pl%isperi(i) = 1
+                        else
+                           pl%isperi(i) = -1
+                        end if
+                     end if
+                  end if
+               end do
+            else
+               do i = 1, npl
+                  if (pl%status(i) == ACTIVE) then
+                     vdotr = dot_product(pl%xb(:,i), pl%vb(:,i))
+                     if (pl%isperi(i) == -1) then
+                        if (vdotr >= 0.0_DP) then
+                           pl%isperi(i) = 0
+                           CALL orbel_xv2aeq(system%Gmtot, pl%xb(:,i), pl%vb(:,i), pl%atp(i), e, pl%peri(i))
+                        end if
+                     else
+                        if (vdotr > 0.0_DP) then
+                           pl%isperi(i) = 1
+                        else
+                           pl%isperi(i) = -1
+                        end if
+                     end if
+                  end if
+               end do
+            end if
+         end if
+      end associate
+
+     return
+   end subroutine symba_util_peri_pl
+
+
    module subroutine symba_util_resize_arr_info(arr, nnew)
       !! author: David A. Minton
       !!
@@ -385,6 +471,7 @@ module subroutine symba_util_resize_tp(self, nnew)
       return
    end subroutine symba_util_resize_tp
 
+
    module subroutine symba_util_sort_pl(self, sortby, ascending)
       !! author: David A. Minton
       !!

src/util/util_coord.f90

@@ -15,21 +15,21 @@ module subroutine util_coord_h2b_pl(self, cb)
       class(swiftest_cb), intent(inout) :: cb   !! Swiftest central body object
       ! Internals
       integer(I4B)  :: i
-      real(DP)      :: msys
+      real(DP)      :: Gmtot
       real(DP), dimension(NDIM) :: xtmp, vtmp
 
       if (self%nbody == 0) return
       associate(pl => self, npl => self%nbody)
-         msys = cb%Gmass
+         Gmtot = cb%Gmass
          xtmp(:) = 0.0_DP
          vtmp(:) = 0.0_DP
          do i = 1, npl
-            msys = msys + pl%Gmass(i)
+            Gmtot = Gmtot + pl%Gmass(i)
             xtmp(:) = xtmp(:) + pl%Gmass(i) * pl%xh(:,i)
             vtmp(:) = vtmp(:) + pl%Gmass(i) * pl%vh(:,i)
          end do
-         cb%xb(:) = -xtmp(:) / msys
-         cb%vb(:) = -vtmp(:) / msys
+         cb%xb(:) = -xtmp(:) / Gmtot
+         cb%vb(:) = -vtmp(:) / Gmtot
          do i = 1, npl
             pl%xb(:,i) = pl%xh(:,i) + cb%xb(:)
             pl%vb(:,i) = pl%vh(:,i) + cb%vb(:)

src/util/util_peri.f90

@@ -45,7 +45,7 @@ module subroutine util_peri_tp(self, system, param)
                if (tp%isperi(i) == -1) then
                   if (vdotr(i) >= 0.0_DP) then
                      tp%isperi(i) = 0
-                     call orbel_xv2aeq(system%msys, tp%xb(:, i), tp%vb(:, i), tp%atp(i), e, tp%peri(i))
+                     call orbel_xv2aeq(system%Gmtot, tp%xb(:, i), tp%vb(:, i), tp%atp(i), e, tp%peri(i))
                   end if
                else
                   if (vdotr(i) > 0.0_DP) then

src/util/util_set.f90

@@ -62,7 +62,7 @@ module subroutine util_set_msys(self)
       ! Arguments
       class(swiftest_nbody_system),  intent(inout) :: self    !! Swiftest nobdy system object
 
-      self%msys = self%cb%Gmass + sum(self%pl%Gmass(1:self%pl%nbody), self%pl%status(1:self%pl%nbody) /= INACTIVE)
+      self%Gmtot = self%cb%Gmass + sum(self%pl%Gmass(1:self%pl%nbody), self%pl%status(1:self%pl%nbody) /= INACTIVE)
 
       return
    end subroutine util_set_msys