diff --git a/src/kick/kick.f90 b/src/kick/kick.f90 index 56e74dfcd..58ebfd1aa 100644 --- a/src/kick/kick.f90 +++ b/src/kick/kick.f90 @@ -49,18 +49,16 @@ module pure subroutine kick_getacch_int_tp(self, GMpl, xhp, npl) ! Internals integer(I4B) :: i, j real(DP) :: rji2, irij3, fac, r2 - real(DP), dimension(NDIM) :: dx, acc + real(DP), dimension(NDIM) :: dx 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) r2 = dot_product(dx(:), dx(:)) fac = GMpl(j) / (r2 * sqrt(r2)) - acc(:) = acc(:) - fac * dx(:) + tp%ah(:, i) = tp%ah(:, i) - fac * dx(:) end do - tp%ah(:, i) = tp%ah(:, i) + acc(:) end do end associate return diff --git a/src/modules/symba_classes.f90 b/src/modules/symba_classes.f90 index 2a9602205..93e9184e2 100644 --- a/src/modules/symba_classes.f90 +++ b/src/modules/symba_classes.f90 @@ -212,6 +212,24 @@ module function symba_encounter_check_tp(self, system, dt, irec) result(lany_enc logical :: lany_encounter !! Returns true if there is at least one close encounter end function symba_encounter_check_tp + module subroutine symba_getacch_pl(self, system, param, t, lbeg) + implicit none + class(symba_pl), intent(inout) :: self !! SyMBA massive body particle data structure + 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) :: t !! Current simulation time + logical, optional, intent(in) :: lbeg !! Optional argument that determines whether or not this is the beginning or end of the step + end subroutine symba_getacch_pl + + module subroutine symba_getacch_tp(self, system, param, t, lbeg) + implicit none + class(symba_tp), intent(inout) :: self !! SyMBA test particle data structure + 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) :: t !! Current time + logical, optional, intent(in) :: lbeg !! Optional argument that determines whether or not this is the beginning or end of the step + end subroutine symba_getacch_tp + module subroutine symba_kick_pltpenc(self, system, dt, irec, sgn) implicit none class(symba_pltpenc), intent(in) :: self !! SyMBA pl-tp encounter list object diff --git a/src/symba/symba_getacch.f90 b/src/symba/symba_getacch.f90 new file mode 100644 index 000000000..b2410d99a --- /dev/null +++ b/src/symba/symba_getacch.f90 @@ -0,0 +1,90 @@ +submodule (symba_classes) s_symba_getacch + use swiftest +contains + module subroutine symba_getacch_pl(self, system, param, t, lbeg) + !! author: David A. Minton + !! + !! Compute heliocentric accelerations of massive bodies + !! + !! Adapted from David E. Kaufmann's Swifter routine symba_getacch.f90 + !! Adapted from Hal Levison's Swift routine symba5_getacch.f + implicit none + ! Arguments + class(symba_pl), intent(inout) :: self !! SyMBA massive body particle data structure + 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) :: t !! Current simulation time + logical, optional, intent(in) :: lbeg !! Optional argument that determines whether or not this is the beginning or end of the step + ! Internals + integer(I4B) :: k + real(DP) :: rji2, rlim2, faci, facj + real(DP), dimension(NDIM) :: dx + + select type(system) + class is (symba_nbody_system) + associate(pl => self, cb => system%cb, plplenc_list => system%plplenc_list, nplplenc => system%plplenc_list%nenc) + call helio_getacch_pl(pl, system, param, t, lbeg) + ! Remove accelerations from encountering pairs + do k = 1, nplplenc + associate(i => plplenc_list%index1(k), j => plplenc_list%index2(k)) + dx(:) = pl%xh(:, j) - pl%xh(:, i) + rji2 = dot_product(dx(:), dx(:)) + rlim2 = (pl%radius(i) + pl%radius(j))**2 + if (rji2 > rlim2) then + irij3 = 1.0_DP / (rji2 * sqrt(rji2)) + faci = pl%Gmass(i) * irij3 + facj = pl%Gmass(j) * irij3 + pl%ah(:, i) = pl%ah(:, i) - facj * dx(:) + pl%ah(:, j) = pl%ah(:, j) + faci * dx(:) + end if + end associate + end do + end associate + end select + + return + end subroutine symba_getacch_pl + + module subroutine symba_getacch_tp(self, system, param, t, lbeg) + !! author: David A. Minton + !! + !! Compute heliocentric accelerations of test particles + !! + !! Adapted from David E. Kaufmann's Swifter routine symba_getacch_tp.f90 + !! Adapted from Hal Levison's Swift routine symba5_getacch.f + implicit none + ! Arguments + class(symba_tp), intent(inout) :: self !! SyMBA test particle data structure + 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) :: t !! Current time + logical, optional, intent(in) :: lbeg !! Optional argument that determines whether or not this is the beginning or end of the step + ! Internals + integer(I4B) :: k + real(DP) :: rji2, fac, rlim2 + real(DP), dimension(NDIM) :: dx + + select type(system) + class is (symba_nbody_system) + associate(tp => self, cb => system%cb, pl => system%pl, pltpenc_list => system%pltpenc_list, npltpenc => system%pltpenc_list%nenc) + call helio_getacch_tp(tp, system, param, t, lbeg) + ! Remove accelerations from encountering pairs + do k = 1, npltpenc + associate(i => pltpenc_list%index1(k), j => pltpenc_list%index2(k)) + if (tp%status(j) == ACTIVE) THEN + dx(:) = tp%xh(:,j) - pl%xh(:,i) + rji2 = dot_product(dx(:), dx(:)) + rlim2 = (pl%radius(i))**2 + if (rji2 > rlim2) then + fac = pl%Gmass(i) / (rji2 * sqrt(rji2)) + tp%ah(:,j) = tp%ah(:,j) + fac * dx(:) + end if + end IF + end associate + end DO + end associate + end select + return + end subroutine symba_getacch_tp + +end submodule s_symba_getacch