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

Commit

Permalink
Fixed issue that occurs when collisions are head-on and there is no i…
Browse files Browse the repository at this point in the history
…nitial orbital angular momentum. This caused a divide-by-zero error when normalizing the momentum vector
  • Loading branch information
daminton committed Oct 21, 2022
1 parent 40232d2 commit b4dabe6
Showing 1 changed file with 8 additions and 2 deletions.
10 changes: 8 additions & 2 deletions src/fraggle/fraggle_set.f90
Original file line number Diff line number Diff line change
Expand Up @@ -163,7 +163,7 @@ module subroutine fraggle_set_coordinate_system(self, colliders)
! Internals
integer(I4B) :: i
real(DP), dimension(NDIM) :: x_cross_v, delta_r, delta_v, Ltot
real(DP) :: r_col_norm, v_col_norm
real(DP) :: r_col_norm, v_col_norm, L_mag
real(DP), dimension(NDIM, self%nbody) :: L_sigma

associate(frag => self, nfrag => self%nbody)
Expand All @@ -176,7 +176,12 @@ module subroutine fraggle_set_coordinate_system(self, colliders)
! and the y-axis aligned with the pre-impact distance vector.
Ltot = colliders%L_orbit(:,1) + colliders%L_orbit(:,2) + colliders%L_spin(:,1) + colliders%L_spin(:,2)
frag%y_coll_unit(:) = delta_r(:) / r_col_norm
frag%z_coll_unit(:) = Ltot(:) / (.mag. Ltot(:))
L_mag = .mag.Ltot(:)
if (L_mag > tiny(L_mag)) then
frag%z_coll_unit(:) = Ltot(:) / L_mag
else
frag%z_coll_unit(:) = 0.0_DP
end if
! The cross product of the y- by z-axis will give us the x-axis
frag%x_coll_unit(:) = frag%y_coll_unit(:) .cross. frag%z_coll_unit(:)

Expand Down Expand Up @@ -228,6 +233,7 @@ module subroutine fraggle_set_natural_scale_factors(self, colliders)
colliders%mass(:) = colliders%mass(:) / frag%mscale
colliders%radius(:) = colliders%radius(:) / frag%dscale
colliders%L_spin(:,:) = colliders%L_spin(:,:) / frag%Lscale
colliders%L_orbit(:,:) = colliders%L_orbit(:,:) / frag%Lscale

do i = 1, 2
colliders%rot(:,i) = colliders%L_spin(:,i) / (colliders%mass(i) * colliders%radius(i)**2 * colliders%Ip(3, i))
Expand Down

0 comments on commit b4dabe6

Please sign in to comment.