From 3323841ab165a9b652fb04239993525741d7a6ad Mon Sep 17 00:00:00 2001 From: David Minton Date: Wed, 19 May 2021 17:17:10 -0400 Subject: [PATCH] Added check for overlap when initializing positions. Also fixed problem that was preventing fragments from being distributed properly by ensuring position vector random number generator returns values from -1 to 1 not just 0 to 1. --- src/symba/symba_frag_pos.f90 | 28 ++++++++++++++++++++++++---- 1 file changed, 24 insertions(+), 4 deletions(-) diff --git a/src/symba/symba_frag_pos.f90 b/src/symba/symba_frag_pos.f90 index 1218f42de..ac973f2d5 100644 --- a/src/symba/symba_frag_pos.f90 +++ b/src/symba/symba_frag_pos.f90 @@ -494,23 +494,43 @@ subroutine set_fragment_position_vectors() !! The initial positions do not conserve energy or momentum, so these need to be adjusted later. implicit none - real(DP) :: r_max + real(DP) :: r_max, dis real(DP), dimension(NDIM) :: L, L_sigma + logical, dimension(:), allocatable :: loverlap + integer(I4B) :: i, j allocate(rmag(nfrag)) allocate(v_r_mag(nfrag)) allocate(v_t_mag(nfrag)) allocate(v_r_unit(NDIM,nfrag)) allocate(v_t_unit(NDIM,nfrag)) + allocate(loverlap(nfrag)) ! Place the fragments into a region that is big enough that we should usually not have overlapping bodies ! An overlapping bodies will collide in the next time step, so it's not a major problem if they do (it just slows the run down) - r_max = 4 * sum(rad_frag(:)) / PI + r_max = 2 * sum(rad_frag(:)) / PI ! We will treat the first fragment of the list as a special case. It gets initialized the maximum distance along the original impactor distance vector. ! This is done because in a regular disruption, the first body is the largest fragment. - x_frag(:, 1) = -y_col_unit(:) * r_max - call random_number(x_frag(:,2:nfrag)) + x_frag(:, 1) = -y_col_unit(:) * r_max + + call random_number(x_frag(:,2:nfrag)) + loverlap(:) = .true. + do while (any(loverlap(2:nfrag))) + do i = 2, nfrag + if (loverlap(i)) then + call random_number(x_frag(:,i)) + x_frag(:, i) = 2 * (x_frag(:, i) - 0.5_DP) * r_max + end if + end do + loverlap(:) = .false. + do j = 1, nfrag + do i = j+1, nfrag + dis = norm2(x_frag(:,j) - x_frag(:,i)) + loverlap(i) = loverlap(i) .or. (dis <= (rad_frag(i) + rad_frag(j))) + end do + end do + end do x_frag(:, :) = x_frag(:, :) * r_max call shift_vector_to_origin(m_frag, x_frag)