Wrote the helio drift subroutine for swiftest_body type and then wrap…

…ped it with interfaces so that they can be used as type-bound procedures for helio_pl and helio_tp

src/helio/helio_drift.f90

@@ -1,18 +1,17 @@
 submodule (helio_classes) s_helio_drift
    use swiftest
 contains
-   module subroutine helio_drift_pl(self, system, param, dt, mask)
+   module subroutine helio_drift_body(self, system, param, dt, mask)
 
       !! author: David A. Minton
       !!
-      !! Loop through massive bodies and call Danby drift routine
-      !! New vectorized version included
+      !! Loop through bodies and call Danby drift routine on democratic heliocentric coordinates
       !!
       !! Adapted from David E. Kaufmann's Swifter routine helio_drift.f90
       !! Adapted from Hal Levison's Swift routine drift.f
       implicit none
       ! Arguments
-      class(helio_pl),              intent(inout) :: self   !! Helio massive body object
+      class(swiftest_body),         intent(inout) :: self   !! Swiftest body object
       class(swiftest_nbody_system), intent(inout) :: system !! Swiftest nbody system object
       class(swiftest_parameters),   intent(in)    :: param  !! Current run configuration parameters 
       real(DP),                     intent(in)    :: dt     !! Stepsize
@@ -51,7 +50,7 @@ module subroutine helio_drift_pl(self, system, param, dt, mask)
          if (any(iflag(1:npl) /= 0)) then
             do i = 1, npl
                if (iflag(i) /= 0) then
-                  write(*, *) " Planet ", self%id(i), " is lost!!!!!!!!!!"
+                  write(*, *) " Body", self%id(i), " is lost!!!!!!!!!!"
                   write(*, *) pl%xh(:,i)
                   write(*, *) pl%vb(:,i)
                   write(*, *) " stopping "
@@ -63,7 +62,39 @@ module subroutine helio_drift_pl(self, system, param, dt, mask)
 
       return
 
+   end subroutine helio_drift_body
+
+   module subroutine helio_drift_pl(self, system, param, dt, mask)
+      !! author: David A. Minton
+      !!
+      !! Wrapper function used to call the body drift routine from a helio_pl structure
+      implicit none
+      ! Arguments
+      class(helio_pl),              intent(inout) :: self   !! Helio massive body object
+      class(swiftest_nbody_system), intent(inout) :: system !! Swiftest nbody system object
+      class(swiftest_parameters),   intent(in)    :: param  !! Current run configuration parameters 
+      real(DP),                     intent(in)    :: dt     !! Stepsize
+      logical, dimension(:),        intent(in)    :: mask   !! Logical mask of size self%nbody that determines which bodies to drift.
+
+      call helio_drift_body(self, system, param, dt, mask)
+      return
    end subroutine helio_drift_pl
+
+   module subroutine helio_drift_tp(self, system, param, dt, mask)
+      !! author: David A. Minton
+      !!
+      !! Wrapper function used to call the body drift routine from a helio_pl structure
+      implicit none
+      ! Arguments
+      class(helio_tp),              intent(inout) :: self   !! Helio massive body object
+      class(swiftest_nbody_system), intent(inout) :: system !! Swiftest nbody system object
+      class(swiftest_parameters),   intent(in)    :: param  !! Current run configuration parameters 
+      real(DP),                     intent(in)    :: dt     !! Stepsize
+      logical, dimension(:),        intent(in)    :: mask   !! Logical mask of size self%nbody that determines which bodies to drift.
+
+      call helio_drift_body(self, system, param, dt, mask)
+      return
+   end subroutine helio_drift_tp
 
    module subroutine helio_drift_linear_pl(self, cb, dt, lbeg)
       !! author: David A. Minton

src/helio/helio_getacch.f90

@@ -51,9 +51,10 @@ module subroutine helio_getacch_tp(self, system, param, t, lbeg)
          real(DP),                     intent(in)    :: t      !! Current time
          logical, optional,            intent(in)    :: lbeg   !! Optional argument that determines whether or not this is the beginning or end of the step
 
-         associate(tp => self, ntp => self%nbody, cb => system%cb, pl => system%pl, npl => system%pl%nbody)
+         associate(tp => self, cb => system%cb, pl => system%pl, npl => system%pl%nbody)
+            tp%ah(:,:) = 0.0_DP
             if (present(lbeg)) system%lbeg = lbeg
-            if (lbeg) then
+            if (system%lbeg) then
                call tp%accel_int(pl%Gmass(:), pl%xbeg(:,:), npl)
             else
                call tp%accel_int(pl%Gmass(:), pl%xend(:,:), npl)

src/kick/kick.f90

@@ -48,17 +48,19 @@ module pure subroutine kick_getacch_int_tp(self, GMpl, xhp, npl)
       integer(I4B),             intent(in)    :: npl  !! Number of active massive bodies
       ! Internals
       integer(I4B)              :: i, j
-      real(DP)                  :: rji2, irij3, fac
+      real(DP)                  :: rji2, irij3, fac, r2
       real(DP), dimension(NDIM) :: dx, acc
 
       associate(tp => self, ntp => self%nbody)
          do concurrent(i = 1:ntp, tp%status(i) == ACTIVE)
             acc(:) = 0.0_DP
             do j = 1, npl
                dx(:) = tp%xh(:,i) - xhp(:, j)
-               rji2 = dot_product(dx(:), dx(:))
-               irij3 = 1.0_DP / (rji2 * sqrt(rji2))
-               fac = GMpl(j) * irij3
+               !rji2 = dot_product(dx(:), dx(:))
+               !irij3 = 1.0_DP / (rji2 * sqrt(rji2))
+               !fac = GMpl(j) * irij3
+               r2 = dot_product(dx(:), dx(:))
+               fac = GMpl(j) / (r2 * sqrt(r2))
                acc(:) = acc(:) - fac * dx(:)
             end do
             tp%ah(:, i) = tp%ah(:, i) + acc(:)

src/modules/helio_classes.f90

@@ -53,6 +53,7 @@ module helio_classes
       procedure, public :: vh2vb    => helio_coord_vh2vb_tp  !! Convert test particles from heliocentric to barycentric coordinates (velocity only)
       procedure, public :: vb2vh    => helio_coord_vb2vh_tp  !! Convert test particles from barycentric to heliocentric coordinates (velocity only)
       procedure, public :: lindrift => helio_drift_linear_tp !! Method for linear drift of massive bodies due to barycentric momentum of Sun
+      procedure, public :: drift    => helio_drift_tp        !! Method for Danby drift in Democratic Heliocentric coordinates 
       procedure, public :: accel    => helio_getacch_tp      !! Compute heliocentric accelerations of massive bodies
       procedure, public :: kick     => helio_kickvb_tp       !! Kicks the barycentric velocities
       procedure, public :: step     => helio_step_tp         !! Steps the body forward one stepsize
@@ -84,6 +85,16 @@ module subroutine helio_coord_vh2vb_tp(self, vbcb)
          class(helio_tp),        intent(inout) :: self !! Helio massive body object
          real(DP), dimension(:), intent(in)    :: vbcb !! Barycentric velocity of the central body
       end subroutine helio_coord_vh2vb_tp
+
+      module subroutine helio_drift_body(self, system, param, dt, mask)
+         use swiftest_classes, only : swiftest_body, swiftest_nbody_system, swiftest_parameters
+         implicit none
+         class(swiftest_body),         intent(inout) :: self   !! Swiftest massive body object
+         class(swiftest_nbody_system), intent(inout) :: system !! Swiftest nbody system object
+         class(swiftest_parameters),   intent(in)    :: param  !! Current run configuration parameters 
+         real(DP),                     intent(in)    :: dt     !! Stepsize
+         logical, dimension(:),        intent(in)    :: mask   !! Logical mask of size self%nbody that determines which bodies to drift
+      end subroutine helio_drift_body
 
       module subroutine helio_drift_pl(self, system, param, dt, mask)
          use swiftest_classes, only : swiftest_nbody_system, swiftest_parameters
@@ -92,9 +103,19 @@ module subroutine helio_drift_pl(self, system, param, dt, mask)
          class(swiftest_nbody_system), intent(inout) :: system !! Swiftest nbody system object
          class(swiftest_parameters),   intent(in)    :: param  !! Current run configuration parameters 
          real(DP),                     intent(in)    :: dt     !! Stepsize
-         logical, dimension(:),         intent(in)    :: mask   !! Logical mask of size self%nbody that determines which bodies to drift
+         logical, dimension(:),        intent(in)    :: mask   !! Logical mask of size self%nbody that determines which bodies to drift
       end subroutine helio_drift_pl
 
+      module subroutine helio_drift_tp(self, system, param, dt, mask)
+         use swiftest_classes, only : swiftest_nbody_system, swiftest_parameters
+         implicit none
+         class(helio_tp),              intent(inout) :: self   !! Helio massive body object
+         class(swiftest_nbody_system), intent(inout) :: system !! Swiftest nbody system object
+         class(swiftest_parameters),   intent(in)    :: param  !! Current run configuration parameters 
+         real(DP),                     intent(in)    :: dt     !! Stepsize
+         logical, dimension(:),        intent(in)    :: mask   !! Logical mask of size self%nbody that determines which bodies to drift
+      end subroutine helio_drift_tp
+
       module subroutine helio_drift_linear_pl(self, cb, dt, lbeg)
          implicit none
          class(helio_pl), intent(inout) :: self  !! Helio massive body object