Fraggle sort of works again. The constraints are not met as well, but…

… it does *something* at least

src/collision/collision_generate.f90

@@ -328,6 +328,8 @@ module subroutine collision_generate_merge(self, nbody_system, param, t)
                   ! The fragment trajectory will be the barycentric trajectory
                   fragments%rb(:,1) = impactors%rbcom(:)
                   fragments%vb(:,1) = impactors%vbcom(:)
+                  fragments%rc(:,1) = 0.0_DP
+                  fragments%vc(:,1) = 0.0_DP
 
                   ! Get the energy of the system after the collision
                   call self%get_energy_and_momentum(nbody_system, param, lbefore=.false.)
@@ -485,7 +487,9 @@ module subroutine collision_generate_simple_rot_vec(collider)
       class(collision_basic), intent(inout) :: collider !! Fraggle collision system object
       ! Internals
       real(DP), dimension(NDIM) :: Lresidual, Lbefore, Lafter, Lspin, rot
-      integer(I4B) :: i
+      real(DP) :: ke_residual, rotmag_old, rotmag_new, vmag_old, vmag_new
+      integer(I4B) :: i, loop
+      integer(I4B) :: MAXLOOP = 10
 
       associate(fragments => collider%fragments, impactors => collider%impactors, nfrag => collider%fragments%nbody)
 
@@ -537,7 +541,7 @@ module subroutine collision_generate_simple_vel_vec(collider)
       real(DP), dimension(2) :: vimp
       real(DP) :: vmag, vdisp
       integer(I4B), dimension(collider%fragments%nbody) :: vsign
-      real(DP), dimension(collider%fragments%nbody) :: vscale
+      real(DP), dimension(collider%fragments%nbody) :: vscale, mass_vscale
 
       associate(fragments => collider%fragments, impactors => collider%impactors, nfrag => collider%fragments%nbody)
          lhitandrun = (impactors%regime == COLLRESOLVE_REGIME_HIT_AND_RUN) 
@@ -554,7 +558,6 @@ module subroutine collision_generate_simple_vel_vec(collider)
             vdisp = 2 * sqrt(2 * impactors%Gmass(2) / impactors%radius(2))
          else
             vdisp = 2 * sqrt(2 * sum(impactors%Gmass(:)) / sum(impactors%radius(:)))
-            !vmag = abs(dot_product(impactors%vb(:,2) - impactors%vb(:,1), impactors%y_unit(:))) 
          end if
 
          vimp(:) = .mag.impactors%vc(:,:)
@@ -566,13 +569,18 @@ module subroutine collision_generate_simple_vel_vec(collider)
          end do
          vscale(:) = vscale(:)/maxval(vscale(:))
 
+         ! Give the fragment velocities a random value that is somewhat scaled with fragment mass
+         call random_number(mass_vscale)
+         mass_vscale(:) = (mass_vscale(:) + 1.0_DP) / 2
+         mass_vscale(:) = mass_vscale(:) * (fragments%mtot / fragments%mass(:))**(0.125_DP)
+         mass_vscale(:) = sqrt(2*mass_vscale(:) / maxval(mass_vscale(:)))
          fragments%vc(:,1) = .mag.impactors%vc(:,1) * impactors%bounce_unit(:) 
          do concurrent(i = 2:nfrag)
 
             j = fragments%origin_body(i)
             vrot(:) = impactors%rot(:,j) .cross. (fragments%rc(:,i) - impactors%rb(:,j) + impactors%rbcom(:))
 
-            vmag = .mag.impactors%vc(:,j) * vscale(i)
+            vmag = .mag.impactors%vc(:,j) * vscale(i) * mass_vscale(i)
 
             if (lhitandrun) then
                fragments%vc(:,i) = vmag * 0.5_DP * impactors%bounce_unit(:) * vsign(i) + vrot(:)

src/collision/collision_util.f90

@@ -216,7 +216,7 @@ module subroutine collision_util_get_kinetic_energy(self)
 
          do i = 1, nfrag
             fragments%ke_orbit = fragments%ke_orbit + fragments%mass(i) * dot_product(fragments%vc(:,i), fragments%vc(:,i))
-            fragments%ke_spin = fragments%ke_spin + fragments%mass(i) * fragments%Ip(3,i) * dot_product(fragments%rot(:,i),fragments%rot(:,i) )
+            fragments%ke_spin = fragments%ke_spin + fragments%mass(i) * fragments%radius(i)**2 * fragments%Ip(3,i) * dot_product(fragments%rot(:,i),fragments%rot(:,i) )
          end do
 
          fragments%ke_orbit = fragments%ke_orbit / 2
@@ -263,6 +263,8 @@ module subroutine collision_util_get_energy_momentum(self,  nbody_system, param,
                   ke_orbit = ke_orbit + impactors%mass(i) * dot_product(impactors%vc(:,i), impactors%vc(:,i))
                   ke_spin = ke_spin + impactors%mass(i) * impactors%radius(i)**2 * impactors%Ip(3,i) * dot_product(impactors%rot(:,i), impactors%rot(:,i))
                end do
+               ke_orbit = ke_orbit / 2
+               ke_spin = ke_spin / 2
             else
                call fragments%get_angular_momentum()
                Lorbit(:) = fragments%Lorbit(:) 

src/fraggle/fraggle_generate.f90

@@ -99,12 +99,12 @@ module subroutine fraggle_generate_disrupt(self, nbody_system, param, t, lfailur
             dEtot = self%Etot(2) - self%Etot(1)
             dLmag = .mag. (self%Ltot(:,2) - self%Ltot(:,1))
 
-            lfailure_local = ((abs(dEtot + impactors%Qloss) > FRAGGLE_ETOL) .or. (dEtot > 0.0_DP)) 
+            lfailure_local = (dEtot > 0.0_DP) 
             if (lfailure_local) then
                write(message, *) dEtot, abs(dEtot + impactors%Qloss) / FRAGGLE_ETOL
-               call swiftest_io_log_one_message(COLLISION_LOG_OUT, "Fraggle failed due to high energy error: " // &
+               call swiftest_io_log_one_message(COLLISION_LOG_OUT, "Fraggle failed due to energy gain: " // &
                                                         trim(adjustl(message)))
-               cycle
+               !cycle
                lfailure_local = .false.
                exit
             end if
@@ -160,8 +160,8 @@ subroutine fraggle_generate_minimize(collider, lfailure)
       class(collision_fraggle), intent(inout) :: collider !! Fraggle collision system object
       logical,                  intent(out)   :: lfailure  !! Logical flag indicating whether this step fails or succeeds
       ! Internals
-      real(DP), parameter :: TOL_MIN = 1.0e-6_DP
-      real(DP), parameter :: TOL_INIT = 1e-8_DP
+      real(DP), parameter :: TOL_MIN = 1.0e-5_DP
+      real(DP), parameter :: TOL_INIT = 1e-6_DP
       integer(I4B), parameter :: MAXLOOP = 50
       real(DP), dimension(collider%fragments%nbody) :: input_v
       real(DP), dimension(:), allocatable :: output_v
@@ -179,11 +179,11 @@ subroutine fraggle_generate_minimize(collider, lfailure)
             Efunc = lambda_obj(E_objective_function)
             tol = TOL_INIT
 
+            fragments%v_r_unit(:,:) = .unit. fragments%vc(:,:)
+            fragments%vmag(:) = .mag. fragments%vc(:,1:nfrag) 
+            fragments%rot(:,1:nfrag) = fragments%rot(:,1:nfrag) * 1e-12_DP
             do while(tol < TOL_MIN)
 
-               fragments%v_r_unit(:,:) = .unit. fragments%vc(:,:)
-               fragments%vmag(:) = .mag. fragments%vc(:,:)
-
                input_v(:) = fragments%vmag(1:nfrag)
                fval = E_objective_function(input_v)
                call minimize_bfgs(Efunc, nelem, input_v, tol, MAXLOOP, lfailure, output_v)
@@ -196,13 +196,15 @@ subroutine fraggle_generate_minimize(collider, lfailure)
                   fragments%vc(:,i) = abs(fragments%vmag(i)) * fragments%v_r_unit(:,i)
                end do
 
-               call collision_util_shift_vector_to_origin(fragments%mass, fragments%vc)
                ! Set spins in order to conserve angular momentum
                call fraggle_generate_set_spins(fragments)
 
                if (.not.lfailure) exit
                tol = tol * 2 ! Keep increasing the tolerance until we converge on a solution
             end do
+
+
+            call collision_util_shift_vector_to_origin(fragments%mass, fragments%vc)
          end select
       end associate
 
@@ -244,10 +246,11 @@ function E_objective_function(val_input) result(fval)
 
                   ! Use the deltaE as the basis of our objective function, with a higher penalty for having excess kinetic energy compared with having a deficit
                   if (deltaE < 0.0_DP) then
-                     fval = deltaE**4
+                     fval = deltaE**8
                   else
                      fval = deltaE**2
                   end if
+                  deallocate(tmp_frag)
 
                end select
             end associate

src/misc/minimizer_module.f90

@@ -128,7 +128,7 @@ module subroutine minimize_bfgs(f, N, x0, eps, maxloop, lerr, x1)
          real(DP), dimension(:), intent(out), allocatable :: x1
          ! Internals
          integer(I4B) ::  i, j, k, l, conv
-         real(DP), parameter     :: graddelta = 1e-4_DP !! Delta x for gradient calculations
+         real(DP), parameter     :: graddelta = 1e-2_DP !! Delta x for gradient calculations
          real(DP), dimension(N) :: S               !! Direction vectors 
          real(DP), dimension(N,N) :: H             !! Approximated inverse Hessian matrix 
          real(DP), dimension(N) :: grad1           !! gradient of f