diff --git a/src/fraggle/fraggle_generate.f90 b/src/fraggle/fraggle_generate.f90 index 75ab1521b..62e182fbb 100644 --- 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)