Found a scheme that works. Computed the available kinetic using the d…

…ifference between the velocity and the escape velocity of individual fragments.

src/fraggle/fraggle_generate.f90

@@ -335,7 +335,6 @@ module subroutine fraggle_generate_pos_vec(collider)
                istart = 2
             end if
 
-
             do i = istart, nfrag
                if (loverlap(i)) then
                   call random_number(phi(i))
@@ -344,8 +343,6 @@ module subroutine fraggle_generate_pos_vec(collider)
                end if
             end do
 
-            ! Make the fragment cloud symmertic about 0
-
             do concurrent(i = istart:nfrag, loverlap(i))
                j = fragments%origin_body(i)
 
@@ -479,9 +476,9 @@ module subroutine fraggle_generate_vel_vec(collider, nbody_system, param, lfailu
       real(DP), dimension(collider%fragments%nbody) :: vscale, ke_rot_remove
       ! 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.5_DP 
-      real(DP)                :: vmax_guess = 10.0_DP
+      real(DP)                :: vmax_guess = 4.0_DP
       real(DP)                :: delta_v, volume
-      integer(I4B), parameter :: MAXLOOP = 100
+      integer(I4B), parameter :: MAXLOOP = 200
       integer(I4B), parameter :: MAXTRY = 1000
       real(DP), parameter :: SUCCESS_METRIC = 1.0e-2_DP
       class(collision_fraggle), allocatable :: collider_local
@@ -562,7 +559,11 @@ module subroutine fraggle_generate_vel_vec(collider, nbody_system, param, lfailu
                do loop = 1, MAXLOOP
                   nsteps = loop * try
                   call collider_local%get_energy_and_momentum(nbody_system, param, phase="after")
-                  ke_avail = max(fragments%ke_orbit_tot - ke_min, 0.0_DP)
+                  ke_avail = 0.0_DP
+                  do i = 1, fragments%nbody
+                     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
                   L_residual(:) = collider_local%L_total(:,2) - collider_local%L_total(:,1)
                   if (ke_avail < epsilon(1.0_DP)) then
@@ -619,7 +620,7 @@ module subroutine fraggle_generate_vel_vec(collider, nbody_system, param, lfailu
                   call fragments%set_coordinate_system()
 
                end do
-               if (dE_best < 0.0_DP) exit outer
+               !if (dE_best < 0.0_DP) exit outer
                ! We didn't converge. Reset the fragment positions and velocities and try a new configuration with some slightly different parameters
                if (fragments%nbody == 2) exit outer
                ! Reduce the number of fragments by one

src/swiftest/swiftest_util.f90

@@ -1865,21 +1865,21 @@ module subroutine swiftest_util_rearray_pl(self, nbody_system, param)
                call plplenc_old%copy(nbody_system%plpl_encounter)
             end if
 
-            ! Re-build the encounter list
-            ! Be sure to get the level info if this is a SyMBA nbody_system
-            select type(nbody_system)
-            class is (symba_nbody_system)
-            select type(pl)
-            class is (symba_pl)
-            select type(tp)
-            class is (symba_tp)
-               lencounter = pl%encounter_check(param, nbody_system, param%dt, nbody_system%irec) 
-               if (tp%nbody > 0) then
-                  lencounter = tp%encounter_check(param, nbody_system, param%dt, nbody_system%irec)
-               end if
-            end select
-            end select
-            end select
+            ! ! Re-build the encounter list
+            ! ! Be sure to get the level info if this is a SyMBA nbody_system
+            ! select type(nbody_system)
+            ! class is (symba_nbody_system)
+            ! select type(pl)
+            ! class is (symba_pl)
+            ! select type(tp)
+            ! class is (symba_tp)
+            !    lencounter = pl%encounter_check(param, nbody_system, param%dt, nbody_system%irec) 
+            !    if (tp%nbody > 0) then
+            !       lencounter = tp%encounter_check(param, nbody_system, param%dt, nbody_system%irec)
+            !    end if
+            ! end select
+            ! end select
+            ! end select
 
             ! Re-index the encounter list as the index values may have changed
             if (nenc_old > 0) then