Merge branch 'debug'

MintonGroup · Feb 13, 2023 · 7b220e3 · 7b220e3
2 parents a4fabd0 + 704eb59
commit 7b220e3
Showing 1 changed file with 6 additions and 3 deletions.
diff --git a/src/fraggle/fraggle_generate.f90 b/src/fraggle/fraggle_generate.f90
@@ -536,7 +536,7 @@ module subroutine fraggle_generate_vel_vec(collider, nbody_system, param, lfailu
       real(DP)  :: MOMENTUM_SUCCESS_METRIC = 10*epsilon(1.0_DP) !! Relative angular momentum error to accept as a success (should be *much* stricter than energy)
       integer(I4B) :: i, j, loop, try, istart, nfrag, nsteps, nsteps_best
       logical :: lhitandrun, lsupercat
-      real(DP), dimension(NDIM) :: vimp_unit, rimp, vrot, L_residual, L_residual_unit, L_residual_best, dL, drot, rot_new, dL_metric
+      real(DP), dimension(NDIM) :: vimp_unit, rimp, vrot, vdisp, L_residual, L_residual_unit, L_residual_best, dL, drot, rot_new, dL_metric
       real(DP) :: vimp, vmag, vesc, dE, E_residual, E_residual_best, E_residual_last, ke_avail, ke_remove, dE_best, fscale, dE_metric, mfrag, rn, dL1_mag, dE_conv
       integer(I4B), dimension(:), allocatable :: vsign
       real(DP), dimension(:), allocatable :: vscale
@@ -569,6 +569,7 @@ module subroutine fraggle_generate_vel_vec(collider, nbody_system, param, lfailu
          vimp = .mag. (impactors%vc(:,2) - impactors%vc(:,1))
 
          E_residual_best = huge(1.0_DP)
+         L_residual_best(:) = 0.0_DP
          lfailure = .false.
          dE_metric = huge(1.0_DP)
          dE_best = huge(1.0_DP)
@@ -581,7 +582,7 @@ module subroutine fraggle_generate_vel_vec(collider, nbody_system, param, lfailu
                if (allocated(dLi_mag)) deallocate(dLi_mag); allocate(dLi_mag(fragments%nbody))
 
                if (lhitandrun) then
-                  vesc = sqrt(2 * sum(fragments%Gmass(istart:fragments%nbody)) / sum(fragments%radius(istart:fragments%nbody)))
+                  vesc = sqrt(2 * sum(fragments%Gmass(istart:fragments%nbody)) / impactors%radius(2))
                   vmin_guess = .mag.impactors%vc(:,2) * (1.0_DP - hitandrun_vscale)
                   vmax_guess = .mag.impactors%vc(:,2) * (1.0_DP + hitandrun_vscale)
                else
@@ -615,7 +616,8 @@ module subroutine fraggle_generate_vel_vec(collider, nbody_system, param, lfailu
                   j = fragments%origin_body(i)
                   vrot(:) = impactors%rot(:,j) .cross. (fragments%rc(:,i) - impactors%rc(:,j))
                   if (lhitandrun) then
-                     fragments%vc(:,i) = vsign(i) * impactors%bounce_unit(:) * vscale(i) + vrot(:)
+                     vdisp(:) = .unit.(fragments%rc(:,i) - impactors%rc(:,2)) * vesc
+                     fragments%vc(:,i) = vsign(i) * impactors%bounce_unit(:) * vscale(i) + vrot(:) + vdisp(:)
                   else
                      vmag = vscale(i) 
                      rimp(:) = fragments%rc(:,i) - impactors%rcimp(:)
@@ -638,6 +640,7 @@ module subroutine fraggle_generate_vel_vec(collider, nbody_system, param, lfailu
                   call collider_local%get_energy_and_momentum(nbody_system, param, phase="after")
                   L_residual(:) = (collider_local%L_total(:,2) - collider_local%L_total(:,1))
                   L_residual_unit(:) = .unit. L_residual(:)
+                  if (nsteps == 1) L_residual_best(:) = L_residual(:)
 
                   ! Use equipartition of spin kinetic energy to distribution spin angular momentum
                   do concurrent(i = istart:fragments%nbody)