Changed how angular momentum residual is handled in hit and run with …

…disruption. Also fixed an issue where hit and run "bounce_unit" vector was not being set correctly because it was being set before the regime was determined.

src/collision/collision_resolve.f90

@@ -152,9 +152,6 @@ module subroutine collision_resolve_consolidate_impactors(self, nbody_system, pa
             ! Destroy the kinship relationships for all members of this impactors%id
             call pl%reset_kinship(impactors%id(:))
 
-            ! Set the coordinate system of the impactors
-            call impactors%set_coordinate_system()
-
          end associate
       end select
       return

src/fraggle/fraggle_generate.f90

@@ -35,6 +35,9 @@ module subroutine fraggle_generate(self, nbody_system, param, t)
       select type(pl => nbody_system%pl)
       class is (swiftest_pl)
          associate(impactors => self%impactors, status => self%status, maxid => nbody_system%maxid)
+            ! Set the coordinate system of the impactors
+            call impactors%set_coordinate_system()
+
             select case (impactors%regime) 
             case (COLLRESOLVE_REGIME_HIT_AND_RUN)
                call self%hitandrun(nbody_system, param, t)
@@ -407,7 +410,7 @@ module subroutine fraggle_generate_rot_vec(collider, nbody_system, param)
       ! Internals
       real(DP), dimension(NDIM) :: Lbefore, Lafter, L_spin, rotdir
       real(DP) :: v_init, v_final, mass_init, mass_final, rotmag, dKE, KE_init, KE_final
-      real(DP), parameter :: random_scale_factor = 0.001_DP !! The relative scale factor to apply to the random component of the rotation vector
+      real(DP), parameter :: random_scale_factor = 0.01_DP !! The relative scale factor to apply to the random component of the rotation vector
       integer(I4B) :: i
 
       associate(fragments => collider%fragments, impactors => collider%impactors, nfrag => collider%fragments%nbody)
@@ -438,16 +441,14 @@ module subroutine fraggle_generate_rot_vec(collider, nbody_system, param)
 
          Lafter(:) = mass_final * (impactors%rb(:,2) - impactors%rb(:,1)) .cross. (v_final * impactors%bounce_unit(:))
          L_spin(:) = impactors%L_spin(:,1) + random_scale_factor * (Lbefore(:) - Lafter(:))
-         !fragments%rot(:,1) = L_spin(:) / (fragments%mass(1) * fragments%radius(1)**2 * fragments%Ip(3,1))
-         !fragments%rotmag(1) = .mag.fragments%rot(:,1)
          ! Add in some random spin noise. The magnitude will be scaled by the before-after amount and the direction will be random
          do i = 2,nfrag
             call random_number(rotdir)
             call random_number(rotmag)
             rotdir = rotdir - 0.5_DP
             rotdir = .unit. rotdir
             fragments%rotmag(i) = random_scale_factor * rotmag * .mag.L_spin(:) / ((nfrag - 1) * fragments%mass(i) * fragments%radius(i)**2 * fragments%Ip(3,i))
-            fragments%rot(:,i) = fragments%rotmag(i) * rotdir
+            fragments%rot(:,i) = fragments%rot(:,i) + fragments%rotmag(i) * rotdir
          end do
       end associate
 
@@ -470,16 +471,16 @@ module subroutine fraggle_generate_vel_vec(collider, nbody_system, param, lfailu
       ! Internals
       integer(I4B) :: i, j, loop, try, istart, nfrag, nsteps
       logical :: lhitandrun, lsupercat
-      real(DP), dimension(NDIM) :: vimp_unit, rimp, vrot, L_residual
-      real(DP) :: vmag, vesc, dE, E_residual, ke_min, ke_avail, ke_remove, dE_best, E_residual_best, fscale, f_spin, f_orbit, dE_metric, dM
+      real(DP), dimension(NDIM) :: vimp_unit, rimp, vrot, L_residual, vdir, L_residual_unit, Li, Lrat, r_lever
+      real(DP) :: rmag, vmag, vesc, dE, E_residual, ke_min, ke_avail, ke_remove, dE_best, E_residual_best, fscale, f_spin, f_orbit, dE_metric, dM, mfrag
       integer(I4B), dimension(collider%fragments%nbody) :: vsign
       real(DP), dimension(collider%fragments%nbody) :: vscale, ke_rot_remove, volume
       ! For the initial "guess" of fragment velocities, this is the minimum and maximum velocity relative to escape velocity that the fragments will have
       real(DP)                :: vmin_guess = 1.01_DP 
       real(DP)                :: vmax_guess = 8.0_DP
       real(DP)                :: delta_v, GC
-      integer(I4B), parameter :: MAXLOOP = 100
-      integer(I4B), parameter :: MAXTRY = 20
+      integer(I4B), parameter :: MAXLOOP = 50
+      integer(I4B), parameter :: MAXTRY = 50
       real(DP), parameter     :: mass_reduction_ratio = 0.1_DP ! Ratio of difference between first and second fragment mass to remove from the largest fragment in case of a failure
       real(DP), parameter :: SUCCESS_METRIC = 1.0e-2_DP
       class(collision_fraggle), allocatable :: collider_local
@@ -506,7 +507,7 @@ module subroutine fraggle_generate_vel_vec(collider, nbody_system, param, lfailu
                vsign(:) = 1
             end where
 
-            ! Hit and run collisions should only affect the runner
+            ! Hit and run collisions should only affect the runners, not the target
             if (lhitandrun) then
                istart = 2
             else 
@@ -546,7 +547,7 @@ module subroutine fraggle_generate_vel_vec(collider, nbody_system, param, lfailu
                      if (i ==1) then
                         fragments%vc(:,i) = impactors%vc(:,1)
                      else
-                        fragments%vc(:,i) = impactors%vc(:,2) + vsign(i) * impactors%bounce_unit(:) * (vscale(i) - 0.5_DP *(maxval(vscale(:)) - minval(vscale(:))))
+                        fragments%vc(:,i) = impactors%vc(:,2) + vsign(i) * impactors%bounce_unit(:) * vscale(i)
                      end if
                   else
                      vmag = vesc * vscale(i) 
@@ -569,24 +570,52 @@ module subroutine fraggle_generate_vel_vec(collider, nbody_system, param, lfailu
                      ke_avail = ke_avail + 0.5_DP * fragments%mass(i) * max(fragments%vmag(i) - vesc,0.0_DP)**2
                   end do
 
-                  ! Check for any residual angular momentum, and if there is any, put it into spin of the largest body
+                  ! Check for any residual angular momentum, and if there is any, put it into spin/shear depending on the type of collision (hit and runs prioritize velocity shear over spin)
                   L_residual(:) = collider_local%L_total(:,2) - collider_local%L_total(:,1)
-                  if (ke_avail < epsilon(1.0_DP)) then
+                  if (lhitandrun .or. (ke_avail < epsilon(1.0_DP))) then
+                     ! Start by putting residual angular momentum into velocity shear
+                     mfrag = sum(fragments%mass(istart:fragments%nbody))
+                     L_residual_unit(:) = .unit. L_residual(:)
+                     fragments%r_unit(:,:) = .unit.fragments%rc(:,:)
                      do i = 1, fragments%nbody
+                        r_lever(:) = (L_residual_unit(:) .cross. fragments%r_unit(:,i))
+                        rmag = .mag.r_lever(:)
+                        if (rmag > epsilon(1.0_DP)) then
+                           vdir(:) = -.unit. r_lever(:)
+                           vmag = .mag.L_residual(:) / (fragments%mtot * .mag.r_lever(:) * fragments%rmag(i))
+                           Li(:) = fragments%mass(i) * fragments%rc(:,i) .cross. (vmag * vdir(:))
+                           Lrat(:) = L_residual(:) / Li(:)
+                           fragments%vc(:,i) = fragments%vc(:,i) + vmag * vdir(:)
+                        end if
+                     end do
+                     fragments%vmag(:) = .mag.fragments%vc(:,:)
+
+                     ! Update the coordinate system now that the velocities have changed
+                     call collision_util_shift_vector_to_origin(fragments%mass, fragments%vc)            
+                     call fragments%set_coordinate_system()
+
+                     ! Try to put as much of the residual angular momentum into the spin of the fragments before the target body
+                     call collider_local%get_energy_and_momentum(nbody_system, param, phase="after")
+                     L_residual(:) = (collider_local%L_total(:,2) - collider_local%L_total(:,1)) 
+                     do i = istart, fragments%nbody
                         fragments%L_spin(:,i) = fragments%L_spin(:,i) - L_residual(:) * fragments%mass(i) / fragments%mtot
                         fragments%rot(:,i) = fragments%L_spin(:,i) / (fragments%mass(i) * fragments%radius(i)**2 * fragments%Ip(:,i)) 
                      end do
-                  else
-                     fragments%L_spin(:,1) = fragments%L_spin(:,1) - L_residual(:) 
-                     fragments%rot(:,1) = fragments%L_spin(:,1) / (fragments%mass(1) * fragments%radius(1)**2 * fragments%Ip(:,1)) 
+                     call collider_local%get_energy_and_momentum(nbody_system, param, phase="after")
+                     L_residual(:) = (collider_local%L_total(:,2) - collider_local%L_total(:,1)) 
                   end if
+
+                  ! Put any remaining residual angular momentum into the spin of the target body
+                  fragments%L_spin(:,1) = fragments%L_spin(:,1) - L_residual(:) 
+                  fragments%rot(:,1) = fragments%L_spin(:,1) / (fragments%mass(1) * fragments%radius(1)**2 * fragments%Ip(:,1)) 
                   fragments%rotmag(:) = .mag.fragments%rot(:,:)
 
                   call collider_local%get_energy_and_momentum(nbody_system, param, phase="after")
                   L_residual(:) = (collider_local%L_total(:,2) - collider_local%L_total(:,1)) 
                   dE = collider_local%te(2) - collider_local%te(1) 
                   E_residual = dE + impactors%Qloss
 
+                  ! Check if we've converged on our constraints
                   if ((abs(E_residual) < abs(E_residual_best)) .or. ((dE < 0.0_DP) .and. (dE_best >= 0.0_DP))) then ! This is our best case so far. Save it for posterity
                      E_residual_best = E_residual
                      dE_best = dE
@@ -601,6 +630,8 @@ module subroutine fraggle_generate_vel_vec(collider, nbody_system, param, lfailu
                   end if
                   if ((dE_best < 0.0_DP) .and. (dE_metric <= SUCCESS_METRIC * try)) exit outer ! As the tries increase, we relax the success metric. What was once a failure might become a success
 
+                  ! If we have an offset of energy from our target, we'll remove it from both orbit and spin, keeping the proportion of each source roughly constant
+
                   ! Remove a constant amount of velocity from the bodies so we don't shift the center of mass and screw up the momentum 
                   f_spin = (fragments%ke_spin_tot )/ (fragments%ke_spin_tot + fragments%ke_orbit_tot)
                   f_orbit = 1.0_DP - f_spin
@@ -614,7 +645,7 @@ module subroutine fraggle_generate_vel_vec(collider, nbody_system, param, lfailu
                   ke_remove = min(f_spin * E_residual, 0.9_DP*fragments%ke_spin_tot)
                   ke_rot_remove(:) = ke_remove * (fragments%ke_spin(:) / fragments%ke_spin_tot)
                   where(ke_rot_remove(:) > fragments%ke_spin(:)) ke_rot_remove(:) = fragments%ke_spin(:) 
-                  do concurrent(i = 1:fragments%nbody, fragments%ke_spin(i) > 10*sqrt(tiny(1.0_DP)))
+                  do concurrent(i = istart:fragments%nbody, fragments%ke_spin(i) > 10*sqrt(tiny(1.0_DP)))
                      fscale = sqrt((fragments%ke_spin(i) - ke_rot_remove(i))/fragments%ke_spin(i))
                      fragments%rotmag(i) = fscale * fragments%rotmag(i)
                      fragments%rot(:,i) = fscale * fragments%rot(:,i)
@@ -648,7 +679,7 @@ module subroutine fraggle_generate_vel_vec(collider, nbody_system, param, lfailu
                call fraggle_generate_rot_vec(collider_local, nbody_system, param)
 
                ! Increase the spatial size factor to get a less dense cloud
-               collider_local%fail_scale = collider_local%fail_scale*1.01_DP
+               collider_local%fail_scale = collider_local%fail_scale*1.001_DP
 
                ! Bring the minimum and maximum velocities closer together 
                delta_v = 0.125_DP * (vmax_guess - vmin_guess)