Fixed bugs that were preventing disruptions from working in the simpl…

…e model
MintonGroup · Dec 24, 2022 · e28c01e · e28c01e
1 parent 09a508d
commit e28c01e
Show file tree

Hide file tree

Showing 2 changed files with 41 additions and 38 deletions.
diff --git a/src/collision/collision_generate.f90 b/src/collision/collision_generate.f90
@@ -209,7 +209,7 @@ module subroutine collision_generate_simple(self, nbody_system, param, t)
       class is (swiftest_nbody_system)
       select type(pl => nbody_system%pl)
       class is (swiftest_pl)
-         associate(impactors => self%impactors, fragments => self%fragments, status => self%status)
+         associate(impactors => self%impactors, status => self%status)
 
             select case (impactors%regime) 
             case (COLLRESOLVE_REGIME_HIT_AND_RUN)
@@ -226,30 +226,34 @@ module subroutine collision_generate_simple(self, nbody_system, param, t)
                write(*,*) "Error in swiftest_collision, unrecognized collision regime"
                call util_exit(FAILURE)
             end select
+            call self%set_mass_dist(param) 
+            call self%disrupt(nbody_system, param, t)
 
             dpe = self%pe(2) - self%pe(1) 
             nbody_system%Ecollisions = nbody_system%Ecollisions - dpe 
             nbody_system%Euntracked  = nbody_system%Euntracked + dpe 
 
-            ! Populate the list of new bodies
-            nfrag = fragments%nbody
-            write(message, *) nfrag
-            call swiftest_io_log_one_message(COLLISION_LOG_OUT, "Generating " // trim(adjustl(message)) // " fragments")
-            select case(impactors%regime)
-            case(COLLRESOLVE_REGIME_DISRUPTION)
-               status = DISRUPTED
-               ibiggest = impactors%id(maxloc(pl%Gmass(impactors%id(:)), dim=1))
-               fragments%id(1) = pl%id(ibiggest)
-               fragments%id(2:nfrag) = [(i, i = param%maxid + 1, param%maxid + nfrag - 1)]
-               param%maxid = fragments%id(nfrag)
-            case(COLLRESOLVE_REGIME_SUPERCATASTROPHIC)
-               status = SUPERCATASTROPHIC
-               fragments%id(1:nfrag) = [(i, i = param%maxid + 1, param%maxid + nfrag)]
-               param%maxid = fragments%id(nfrag)
-            end select
+            associate (fragments => self%fragments)
+               ! Populate the list of new bodies
+               nfrag = fragments%nbody
+               write(message, *) nfrag
+               call swiftest_io_log_one_message(COLLISION_LOG_OUT, "Generating " // trim(adjustl(message)) // " fragments")
+               select case(impactors%regime)
+               case(COLLRESOLVE_REGIME_DISRUPTION)
+                  status = DISRUPTED
+                  ibiggest = impactors%id(maxloc(pl%Gmass(impactors%id(:)), dim=1))
+                  fragments%id(1) = pl%id(ibiggest)
+                  fragments%id(2:nfrag) = [(i, i = param%maxid + 1, param%maxid + nfrag - 1)]
+                  param%maxid = fragments%id(nfrag)
+               case(COLLRESOLVE_REGIME_SUPERCATASTROPHIC)
+                  status = SUPERCATASTROPHIC
+                  fragments%id(1:nfrag) = [(i, i = param%maxid + 1, param%maxid + nfrag)]
+                  param%maxid = fragments%id(nfrag)
+               end select
 
-            call collision_resolve_mergeaddsub(nbody_system, param, t, status)
-      end associate
+               call collision_resolve_mergeaddsub(nbody_system, param, t, status)
+            end associate
+         end associate
       end select
       end select
       return

diff --git a/src/collision/collision_resolve.f90 b/src/collision/collision_resolve.f90
@@ -507,22 +507,21 @@ module subroutine collision_resolve_plpl(self, nbody_system, param, t, dt, irec)
          associate(plpl_encounter => self, plpl_collision => nbody_system%plpl_collision, &
                    collision_history => nbody_system%collision_history, pl => nbody_system%pl, cb => nbody_system%cb, &
                    collider => nbody_system%collider, fragments => nbody_system%collider%fragments, impactors => nbody_system%collider%impactors)
-            associate( idx1 => plpl_collision%index1, idx2 => plpl_collision%index2)
+            if (plpl_collision%nenc == 0) return ! No collisions to resolve
 
-               if (plpl_collision%nenc == 0) return ! No collisions to resolve
 
+            ! Make sure that the heliocentric and barycentric coordinates are consistent with each other
+            call pl%vb2vh(nbody_system%cb) 
+            call pl%rh2rb(nbody_system%cb)
 
-               ! Make sure that the heliocentric and barycentric coordinates are consistent with each other
-               call pl%vb2vh(nbody_system%cb) 
-               call pl%rh2rb(nbody_system%cb)
-
-               ! Get the energy before the collision is resolved
-               if (param%lenergy) then
-                  call nbody_system%get_energy_and_momentum(param)
-                  Eorbit_before = nbody_system%te
-               end if
+            ! Get the energy before the collision is resolved
+            if (param%lenergy) then
+               call nbody_system%get_energy_and_momentum(param)
+               Eorbit_before = nbody_system%te
+            end if
 
-               do loop = 1, MAXCASCADE
+            do loop = 1, MAXCASCADE
+               associate( idx1 => plpl_collision%index1, idx2 => plpl_collision%index2)
                   ncollisions = plpl_collision%nenc
                   write(timestr,*) t
                   call swiftest_io_log_one_message(COLLISION_LOG_OUT, "")
@@ -572,15 +571,15 @@ module subroutine collision_resolve_plpl(self, nbody_system, param, t, dt, irec)
                      write(*,*) "Consider reducing the step size or changing the parameters in the collisional model to reduce the number of fragments."
                      call util_exit(FAILURE)
                   end if
-               end do
+               end associate
+            end do
 
-               if (param%lenergy) then
-                  call nbody_system%get_energy_and_momentum(param)
-                  Eorbit_after = nbody_system%te
-                  nbody_system%Ecollisions = nbody_system%Ecollisions + (Eorbit_after - Eorbit_before)
-               end if
+            if (param%lenergy) then
+               call nbody_system%get_energy_and_momentum(param)
+               Eorbit_after = nbody_system%te
+               nbody_system%Ecollisions = nbody_system%Ecollisions + (Eorbit_after - Eorbit_before)
+            end if
 
-            end associate
          end associate
       end select
       end select