From e9381edef51f38c60bba8e1d940a03f5d54bfbb4 Mon Sep 17 00:00:00 2001 From: David A Minton Date: Sat, 18 Feb 2023 16:08:39 -0500 Subject: [PATCH] Fixed major issue that was causing the impactor positions to be shifted waaaay too far away from the impact point when trying to keep them from overlapping. The procedure is now greatly simplified so that bodies are only shifted just barely past their overlapping distance along their center of mass. --- src/collision/collision_resolve.f90 | 23 ++++++++++++----------- 1 file changed, 12 insertions(+), 11 deletions(-) diff --git a/src/collision/collision_resolve.f90 b/src/collision/collision_resolve.f90 index b9cf5663a..e5c4a86b6 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, b, r_dot_vunit, rrel_mag, vrel_mag, rlim, dt - real(DP), dimension(NDIM) :: xc, vc, xcom, vcom, xchild, vchild, xcrossv, rrel, vrel, vrel_unit + 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), dimension(NDIM,2) :: mxc, vcc select type(nbody_system) @@ -141,22 +141,23 @@ module subroutine collision_resolve_consolidate_impactors(self, nbody_system, pa end do lflag = .true. + mtot = sum(impactors%mass(:)) + xcom(:) = (impactors%mass(1) * impactors%rb(:, 1) + impactors%mass(2) * impactors%rb(:, 2)) / mtot + vcom(:) = (impactors%mass(1) * impactors%vb(:, 1) + impactors%mass(2) * impactors%vb(:, 2)) / mtot + ! Shift the impactors so that they are not overlapping rlim = sum(impactors%radius(1:2)) rrel = impactors%rb(:,2) - impactors%rb(:,1) rrel_mag = .mag. rrel if (rrel_mag < rlim) then - vrel = impactors%vb(:,2) - impactors%vb(:,1) - vrel_mag = .mag. vrel - vrel_unit = .unit. vrel - r_dot_vunit = dot_product(rrel,vrel_unit) - b = sqrt(rrel_mag**2 - r_dot_vunit**2) - dt = (sqrt(rlim**2 - b**2) + r_dot_vunit)/vrel_mag - impactors%rb(:,1:2) = impactors%rb(:,1:2) - dt * impactors%vb(:,1:2) + rrel_unit = .unit.rrel + dr(:) = (1.0_DP + 2*epsilon(1.0_DP)) * (rlim - rrel_mag) * rrel_unit(:) + 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) + rrel_mag = .mag. rrel end if - xcom(:) = (impactors%mass(1) * impactors%rb(:, 1) + impactors%mass(2) * impactors%rb(:, 2)) / sum(impactors%mass(:)) - vcom(:) = (impactors%mass(1) * impactors%vb(:, 1) + impactors%mass(2) * impactors%vb(:, 2)) / sum(impactors%mass(:)) mxc(:, 1) = impactors%mass(1) * (impactors%rb(:, 1) - xcom(:)) mxc(:, 2) = impactors%mass(2) * (impactors%rb(:, 2) - xcom(:)) vcc(:, 1) = impactors%vb(:, 1) - vcom(:)