Skip to content
This repository was archived by the owner on Aug 28, 2024. It is now read-only.

Commit

Permalink
Fixed distance scaling to be based on mutual radius
Browse files Browse the repository at this point in the history
  • Loading branch information
daminton committed Jun 4, 2021
1 parent a5ab709 commit ee2e682
Showing 1 changed file with 6 additions and 5 deletions.
11 changes: 6 additions & 5 deletions src/symba/symba_frag_pos.f90
Original file line number Diff line number Diff line change
Expand Up @@ -439,7 +439,7 @@ 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, dis
real(DP) :: r_max, dis, rad
real(DP), dimension(NDIM) :: L_sigma
logical, dimension(:), allocatable :: loverlap
integer(I4B) :: i, j
Expand All @@ -449,19 +449,20 @@ subroutine set_fragment_position_vectors()
! 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 = r_max_start
rad = sum(radius(:))

! 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
x_frag(:, 2) = y_col_unit(:) * r_max
x_frag(:, 1) = -y_col_unit(:) * r_max * rad
x_frag(:, 2) = y_col_unit(:) * r_max * rad

call random_number(x_frag(:,3:nfrag))
loverlap(:) = .true.
do while (any(loverlap(3:nfrag)))
do i = 3, nfrag
if (loverlap(i)) then
call random_number(x_frag(:,i))
x_frag(:, i) = 2 * (x_frag(:, i) - 0.5_DP) * r_max
x_frag(:, i) = 2 * (x_frag(:, i) - 0.5_DP) * r_max * rad
end if
end do
loverlap(:) = .false.
Expand Down Expand Up @@ -613,7 +614,7 @@ function tangential_objective_function(v_t_mag_input, lerr) result(fval)
dL = norm2(L_frag(:) - L_frag_orb(:)) / norm2(L_frag_orb(:))
keo = 0.5_DP * keo
kes = 0.5_DP * kes
fval = (keo + kes) / ke_frag_budget
fval = keo + kes
lerr = .false.
return
end function tangential_objective_function
Expand Down

0 comments on commit ee2e682

Please sign in to comment.