diff --git a/src/collision/collision_generate.f90 b/src/collision/collision_generate.f90 index 4ef6ec215..f0fc6d9e7 100644 --- a/src/collision/collision_generate.f90 +++ b/src/collision/collision_generate.f90 @@ -172,7 +172,7 @@ module subroutine collision_generate_merge(self, nbody_system, param, t) real(DP), intent(in) :: t !! The time of the collision ! Internals integer(I4B) :: i, j, k, ibiggest - real(DP), dimension(NDIM) :: L_spin_new + real(DP), dimension(NDIM) :: L_spin_new, L_residual real(DP) :: volume character(len=STRMAX) :: message @@ -228,6 +228,9 @@ module subroutine collision_generate_merge(self, nbody_system, param, t) ! Get the energy of the system after the collision call self%get_energy_and_momentum(nbody_system, param, phase="after") + L_residual(:) = (self%L_total(:,2) - self%L_total(:,1)) + call collision_util_velocity_torque(-L_residual(:), fragments%mass(1), fragments%rb(:,1), fragments%vb(:,1)) + ! Update any encounter lists that have the removed bodies in them so that they instead point to the new body do k = 1, nbody_system%plpl_encounter%nenc do j = 1, impactors%ncoll diff --git a/src/collision/collision_module.f90 b/src/collision/collision_module.f90 index d6323f522..2957dffaf 100644 --- a/src/collision/collision_module.f90 +++ b/src/collision/collision_module.f90 @@ -521,6 +521,13 @@ module subroutine collision_util_setup_fragments(self, n) integer(I4B), intent(in) :: n !! Number of particles to allocate space for end subroutine collision_util_setup_fragments + module subroutine collision_util_velocity_torque(dL, mass, r, v) + implicit none + real(DP), dimension(:), intent(in) :: dL !! Change in angular momentum to apply + real(DP), intent(in) :: mass !! Mass of body + real(DP), dimension(:), intent(in) :: r !! Position of body wrt system center of mass + real(DP), dimension(:), intent(inout) :: v !! Velocity of body wrt system center of mass + end subroutine collision_util_velocity_torque end interface contains diff --git a/src/collision/collision_util.f90 b/src/collision/collision_util.f90 index 63d623de4..7adbe2ec5 100644 --- a/src/collision/collision_util.f90 +++ b/src/collision/collision_util.f90 @@ -1028,4 +1028,32 @@ module subroutine collision_util_set_original_scale_factors(self) end subroutine collision_util_set_original_scale_factors + module subroutine collision_util_velocity_torque(dL, mass, r, v) + !! author: David A. Minton + !! + !! Applies a torque to a body's center of mass velocity given a change in angular momentum + implicit none + ! Arguments + real(DP), dimension(:), intent(in) :: dL !! Change in angular momentum to apply + real(DP), intent(in) :: mass !! Mass of body + real(DP), dimension(:), intent(in) :: r !! Position of body wrt system center of mass + real(DP), dimension(:), intent(inout) :: v !! Velocity of body wrt system center of mass + ! Internals + real(DP), dimension(NDIM) :: dL_unit, r_unit, r_lever, vapply + real(DP) :: rmag, r_lever_mag + + dL_unit(:) = .unit. dL + r_unit(:) = .unit.r(:) + rmag = .mag.r(:) + ! Project the position vector onto the plane defined by the angular momentum vector and the origin to get the "lever arm" distance + r_lever(:) = dL_unit(:) .cross. (r(:) .cross. dL_unit(:)) + r_lever_mag = .mag.r_lever(:) + if ((r_lever_mag > epsilon(1.0_DP)) .and. (rmag > epsilon(1.0_DP))) then + vapply(:) = (dL(:) .cross. r(:)) / (mass * rmag**2) + v(:) = v(:) + vapply(:) + end if + + return + end subroutine collision_util_velocity_torque + end submodule s_collision_util \ No newline at end of file diff --git a/src/fraggle/fraggle_generate.f90 b/src/fraggle/fraggle_generate.f90 index 315c8d8b5..8b805cacb 100644 --- a/src/fraggle/fraggle_generate.f90 +++ b/src/fraggle/fraggle_generate.f90 @@ -602,7 +602,7 @@ module subroutine fraggle_generate_vel_vec(collider, nbody_system, param, lfailu end if dL(:) = -L_residual(:) * fragments%mass(i) / fragments%mtot - drot(:) * fragments%Ip(3,i) * fragments%mass(i) * fragments%radius(i)**2 - call fraggle_generate_velocity_torque(dL, fragments%mass(i), fragments%rc(:,i), fragments%vc(:,i)) + call collision_util_velocity_torque(dL, fragments%mass(i), fragments%rc(:,i), fragments%vc(:,i)) call collision_util_shift_vector_to_origin(fragments%mass, fragments%vc) fragments%vmag(i) = .mag.fragments%vc(:,i) @@ -670,7 +670,7 @@ module subroutine fraggle_generate_vel_vec(collider, nbody_system, param, lfailu end do outer lfailure = (dE_best > 0.0_DP) - call fraggle_generate_velocity_torque(-L_residual_best(:), fragments%mtot, impactors%rbcom, impactors%vbcom) + call collision_util_velocity_torque(-L_residual_best(:), fragments%mtot, impactors%rbcom, impactors%vbcom) do concurrent(i = 1:collider%fragments%nbody) collider%fragments%vb(:,i) = collider%fragments%vc(:,i) + impactors%vbcom(:) @@ -695,33 +695,4 @@ module subroutine fraggle_generate_vel_vec(collider, nbody_system, param, lfailu end subroutine fraggle_generate_vel_vec - subroutine fraggle_generate_velocity_torque(dL, mass, r, v) - !! author: David A. Minton - !! - !! Applies a torque to a body's center of mass velocity given a change in angular momentum - implicit none - ! Arguments - real(DP), dimension(:), intent(in) :: dL !! Change in angular momentum to apply - real(DP), intent(in) :: mass !! Mass of body - real(DP), dimension(:), intent(in) :: r !! Position of body wrt system center of mass - real(DP), dimension(:), intent(inout) :: v !! Velocity of body wrt system center of mass - ! Internals - real(DP), dimension(NDIM) :: dL_unit, r_unit, r_lever, vapply - real(DP) :: rmag, r_lever_mag - - dL_unit(:) = .unit. dL - r_unit(:) = .unit.r(:) - rmag = .mag.r(:) - ! Project the position vector onto the plane defined by the angular momentum vector and the origin to get the "lever arm" distance - r_lever(:) = dL_unit(:) .cross. (r(:) .cross. dL_unit(:)) - r_lever_mag = .mag.r_lever(:) - if ((r_lever_mag > epsilon(1.0_DP)) .and. (rmag > epsilon(1.0_DP))) then - vapply(:) = (dL(:) .cross. r(:)) / (mass * rmag**2) - v(:) = v(:) + vapply(:) - end if - - return - end subroutine fraggle_generate_velocity_torque - - end submodule s_fraggle_generate