From eba83764923994a5109f93e46a598af2aebe11f1 Mon Sep 17 00:00:00 2001 From: David A Minton Date: Sun, 19 Feb 2023 02:03:37 -0500 Subject: [PATCH] Improved overlap adjustment --- src/collision/collision_resolve.f90 | 12 ++++++++---- 1 file changed, 8 insertions(+), 4 deletions(-) diff --git a/src/collision/collision_resolve.f90 b/src/collision/collision_resolve.f90 index e5c4a86b6..dadde1dfd 100644 --- a/src/collision/collision_resolve.f90 +++ b/src/collision/collision_resolve.f90 @@ -33,8 +33,8 @@ module subroutine collision_resolve_consolidate_impactors(self, nbody_system, pa integer(I4B), dimension(2) :: nchild integer(I4B) :: i, j, nimpactors, idx_child real(DP), dimension(2) :: volume, density - real(DP) :: mchild, volchild, rrel_mag, rlim, dt, mtot - real(DP), dimension(NDIM) :: xc, vc, xcom, vcom, xchild, vchild, xcrossv, rrel, rrel_unit, dr + real(DP) :: mchild, volchild, rrel_mag, rlim, mtot, vdotr + real(DP), dimension(NDIM) :: xc, vc, xcom, vcom, xchild, vchild, xcrossv, rrel, vrel, rrel_unit, vrel_unit, dr real(DP), dimension(NDIM,2) :: mxc, vcc select type(nbody_system) @@ -151,7 +151,11 @@ module subroutine collision_resolve_consolidate_impactors(self, nbody_system, pa rrel_mag = .mag. rrel if (rrel_mag < rlim) then rrel_unit = .unit.rrel - dr(:) = (1.0_DP + 2*epsilon(1.0_DP)) * (rlim - rrel_mag) * rrel_unit(:) + vrel = impactors%vb(:,2) - impactors%vb(:,1) + vrel_unit = .unit.vrel + vdotr = dot_product(vrel_unit, rrel) + dr(:) = -(vdotr - sign(1.0_DP, vdotr) * sqrt(rlim**2 - rrel_mag**2 + vdotr**2)) * vrel_unit(:) + dr(:) = (1.0_DP + 2*epsilon(1.0_DP)) * dr(:) impactors%rb(:,1) = impactors%rb(:,1) - dr(:) * impactors%mass(2) / mtot impactors%rb(:,2) = impactors%rb(:,2) + dr(:) * impactors%mass(1) / mtot rrel = impactors%rb(:,2) - impactors%rb(:,1) @@ -684,4 +688,4 @@ module subroutine collision_resolve_pltp(self, nbody_system, param, t, dt, irec) return end subroutine collision_resolve_pltp -end submodule s_collision_resolve \ No newline at end of file +end submodule s_collision_resolve