diff --git a/src/modules/swiftest_classes.f90 b/src/modules/swiftest_classes.f90 index 686959d87..4cf180e84 100644 --- a/src/modules/swiftest_classes.f90 +++ b/src/modules/swiftest_classes.f90 @@ -226,6 +226,7 @@ module swiftest_classes procedure :: append => util_append_pl !! Appends elements from one structure to another procedure :: h2b => util_coord_h2b_pl !! Convert massive bodies from heliocentric to barycentric coordinates (position and velocity) procedure :: b2h => util_coord_b2h_pl !! Convert massive bodies from barycentric to heliocentric coordinates (position and velocity) + procedure :: xh2xb => util_coord_xh2xb_pl !! Convert position vectors of massive bodies from heliocentric to barycentric coordinates (position and velocity) procedure :: fill => util_fill_pl !! "Fills" bodies from one object into another depending on the results of a mask (uses the UNPACK intrinsic) procedure :: resize => util_resize_pl !! Checks the current size of a Swiftest body against the requested size and resizes it if it is too small. procedure :: set_beg_end => util_set_beg_end_pl !! Sets the beginning and ending positions and velocities of planets. @@ -936,6 +937,12 @@ module subroutine util_coord_h2b_tp(self, cb) class(swiftest_cb), intent(in) :: cb !! Swiftest central body object end subroutine util_coord_h2b_tp + module subroutine util_coord_xh2xb_pl(self, cb) + implicit none + class(swiftest_pl), intent(inout) :: self !! Swiftest massive body object + class(swiftest_cb), intent(inout) :: cb !! Swiftest central body object + end subroutine util_coord_xh2xb_pl + module subroutine util_copy_encounter(self, source) implicit none class(swiftest_encounter), intent(inout) :: self !! Encounter list diff --git a/src/symba/symba_collision.f90 b/src/symba/symba_collision.f90 index 0aec3a1be..5ccbcf993 100644 --- a/src/symba/symba_collision.f90 +++ b/src/symba/symba_collision.f90 @@ -1076,7 +1076,7 @@ module subroutine symba_collision_resolve_plplenc(self, system, param, t) write(*, *) "Collision between massive bodies detected at time t = ", t call pl%vb2vh(system%cb) - call pl%h2b(system%cb) + call pl%xh2xb(system%cb) if (param%lfragmentation) then call plplcollision_list%resolve_fragmentations(system, param) else diff --git a/src/util/util_coord.f90 b/src/util/util_coord.f90 index 2a970d0dc..4daad0e53 100644 --- a/src/util/util_coord.f90 +++ b/src/util/util_coord.f90 @@ -42,6 +42,42 @@ module subroutine util_coord_h2b_pl(self, cb) end subroutine util_coord_h2b_pl + module subroutine util_coord_xh2xb_pl(self, cb) + !! author: David A. Minton + !! + !! Convert position vectors of massive bodies from heliocentric to barycentric coordinates (position and velocity) + !! + !! Adapted from David E. Kaufmann's Swifter routine coord_h2b.f90 + !! Adapted from Hal Levison's Swift routine coord_h2b.f + implicit none + ! Arguments + class(swiftest_pl), intent(inout) :: self !! Swiftest massive body object + class(swiftest_cb), intent(inout) :: cb !! Swiftest central body object + ! Internals + integer(I4B) :: i + real(DP) :: Gmtot + real(DP), dimension(NDIM) :: xtmp + + if (self%nbody == 0) return + associate(pl => self, npl => self%nbody) + Gmtot = cb%Gmass + xtmp(:) = 0.0_DP + do i = 1, npl + if (pl%status(i) == INACTIVE) cycle + Gmtot = Gmtot + pl%Gmass(i) + xtmp(:) = xtmp(:) + pl%Gmass(i) * pl%xh(:,i) + end do + cb%xb(:) = -xtmp(:) / Gmtot + do i = 1, npl + if (pl%status(i) == INACTIVE) cycle + pl%xb(:,i) = pl%xh(:,i) + cb%xb(:) + end do + end associate + + return + end subroutine util_coord_xh2xb_pl + + module subroutine util_coord_h2b_tp(self, cb) !! author: David A. Minton !!