Merge branch 'debug'

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

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

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

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