Tweaked Fraggle so that disruption fragments act like ejecta and fly …

…away from the impact point

src/fraggle/fraggle_generate.f90

@@ -72,7 +72,7 @@ module subroutine fraggle_generate_fragments(self, colliders, system, param, lfa
          call frag%get_energy_and_momentum(colliders, system, param, lbefore=.true.)
 
          ! Start out the fragments close to the initial separation distance. This will be increased if there is any overlap or we fail to find a solution
-         r_max_start = 1 * norm2(colliders%rb(:,2) - colliders%rb(:,1))
+         r_max_start = 2 * norm2(colliders%rb(:,2) - colliders%rb(:,1))
          lfailure = .false.
          try = 1
          do while (try < MAXTRY)
@@ -177,8 +177,10 @@ subroutine fraggle_generate_pos_vec(frag, colliders, r_max_start)
       real(DP),                 intent(in)    :: r_max_start !! Initial guess for the starting maximum radial distance of fragments
       ! Internals
       real(DP)  :: dis, rad, r_max
+      real(DP), dimension(NDIM) :: runit
       logical, dimension(:), allocatable :: loverlap
       integer(I4B) :: i, j
+      logical :: lfixdir
 
       associate(nfrag => frag%nbody)
          allocate(loverlap(nfrag))
@@ -188,9 +190,13 @@ subroutine fraggle_generate_pos_vec(frag, colliders, r_max_start)
          r_max = r_max_start
          rad = sum(colliders%radius(:))
 
+         lfixdir = (frag%regime /= COLLRESOLVE_REGIME_SUPERCATASTROPHIC) ! For this style of impact, make the fragments act like eject and point away from the impact point
+
+         runit(:) = colliders%rb(:,2) - colliders%rb(:,1)
+         runit(:) = runit(:) / (.mag. runit(:))
+
          ! We will treat the first two fragments of the list as special cases. They get initialized the maximum distances apart along the original impactor distance vector.
          ! This is done because in a regular disruption, the first body is the largest, the second the second largest, and the rest are smaller equal-mass fragments.
-
          call random_number(frag%x_coll(:,3:nfrag))
          loverlap(:) = .true.
          do while (any(loverlap(3:nfrag)))
@@ -201,6 +207,8 @@ subroutine fraggle_generate_pos_vec(frag, colliders, r_max_start)
                if (loverlap(i)) then
                   call random_number(frag%x_coll(:,i))
                   frag%x_coll(:, i) = 2 * (frag%x_coll(:, i) - 0.5_DP) * r_max 
+                  frag%x_coll(:, i) = frag%x_coll(:, i) + (frag%rbimp(:) - frag%rbcom(:)) ! Shift the center of the fragment cloud to the impact point rather than the CoM
+                  if (lfixdir .and. dot_product(frag%x_coll(:,i), runit(:)) < 0.0_DP) frag%x_coll(:, i) = -frag%x_coll(:, i) ! Make sure the fragment cloud points away from the impact point
                end if
             end do
             loverlap(:) = .false.
@@ -214,12 +222,12 @@ subroutine fraggle_generate_pos_vec(frag, colliders, r_max_start)
          call fraggle_util_shift_vector_to_origin(frag%mass, frag%x_coll)
          call frag%set_coordinate_system(colliders)
 
-         do i = 1, nfrag
+         do concurrent(i = 1:nfrag)
             frag%rb(:,i) = frag%x_coll(:,i) + frag%rbcom(:)
          end do
 
          frag%rbcom(:) = 0.0_DP
-         do i = 1, nfrag
+         do concurrent(i = 1:nfrag)
             frag%rbcom(:) = frag%rbcom(:) + frag%mass(i) * frag%rb(:,i) 
          end do
          frag%rbcom(:) = frag%rbcom(:) / frag%mtot

src/fraggle/fraggle_regime.f90

@@ -27,7 +27,7 @@ module subroutine fraggle_regime_colliders(self, frag, system, param)
       integer(I4B) :: jtarg, jproj
       real(DP), dimension(2) :: radius_si, mass_si, density_si
       real(DP) :: min_mfrag_si, Mcb_si
-      real(DP), dimension(NDIM)  :: x1_si, v1_si, x2_si, v2_si
+      real(DP), dimension(NDIM)  :: x1_si, v1_si, x2_si, v2_si, runit
       real(DP) :: mlr, mslr, mtot, dentot
 
       associate(colliders => self)
@@ -71,6 +71,11 @@ module subroutine fraggle_regime_colliders(self, frag, system, param)
          frag%rbcom(:) = (colliders%mass(1) * colliders%rb(:,1) + colliders%mass(2) * colliders%rb(:,2)) / frag%mtot 
          frag%vbcom(:) = (colliders%mass(1) * colliders%vb(:,1) + colliders%mass(2) * colliders%vb(:,2)) / frag%mtot
 
+         ! Find the point of impact between the two bodies
+         runit(:) = colliders%rb(:,2) - colliders%rb(:,1)
+         runit(:) = runit(:) / (.mag. runit(:))
+         frag%rbimp(:) = colliders%rb(:,1) + colliders%radius(1) * runit(:)
+
          ! Convert quantities back to the system units and save them into the fragment system
          frag%mass_dist(:) = (frag%mass_dist(:) / param%MU2KG) 
          frag%Qloss = frag%Qloss * (param%TU2S / param%DU2M)**2 / param%MU2KG

src/fraggle/fraggle_set.f90

@@ -185,10 +185,11 @@ module subroutine fraggle_set_coordinate_system(self, colliders)
          Ltot = colliders%L_orbit(:,1) + colliders%L_orbit(:,2) + colliders%L_spin(:,1) + colliders%L_spin(:,2)
          frag%y_coll_unit(:) = delta_r(:) / r_col_norm 
          L_mag = .mag.Ltot(:)
-         if (L_mag > tiny(L_mag)) then
+         if (L_mag > sqrt(tiny(L_mag))) then
             frag%z_coll_unit(:) = Ltot(:) / L_mag
          else
-            frag%z_coll_unit(:) = 0.0_DP
+            call random_number(frag%z_coll_unit(:))
+            frag%z_coll_unit(:) = frag%z_coll_unit(:) / (.mag.frag%z_coll_unit(:))
          end if
          ! The cross product of the y- by z-axis will give us the x-axis
          frag%x_coll_unit(:) = frag%y_coll_unit(:) .cross. frag%z_coll_unit(:)
@@ -198,13 +199,15 @@ module subroutine fraggle_set_coordinate_system(self, colliders)
 
          call random_number(L_sigma(:,:)) ! Randomize the tangential velocity direction. This helps to ensure that the tangential velocity doesn't completely line up with the angular momentum vector,
                                           ! otherwise we can get an ill-conditioned system
+
          do concurrent(i = 1:nfrag, frag%rmag(i) > 0.0_DP)
             frag%v_r_unit(:, i) = frag%x_coll(:, i) / frag%rmag(i)
             frag%v_n_unit(:, i) = frag%z_coll_unit(:) + 2e-1_DP * (L_sigma(:,i) - 0.5_DP)
             frag%v_n_unit(:, i) = frag%v_n_unit(:, i) / (.mag. frag%v_n_unit(:, i))
             frag%v_t_unit(:, i) = frag%v_n_unit(:, i) .cross. frag%v_r_unit(:, i)
             frag%v_t_unit(:, i) = frag%v_t_unit(:, i) / (.mag. frag%v_t_unit(:, i))
          end do
+
       end associate
 
       return
@@ -236,6 +239,7 @@ module subroutine fraggle_set_natural_scale_factors(self, colliders)
          ! Scale all dimensioned quantities of colliders and fragments
          frag%rbcom(:) = frag%rbcom(:) / frag%dscale
          frag%vbcom(:) = frag%vbcom(:) / frag%vscale
+         frag%rbimp(:) = frag%rbimp(:) / frag%dscale
          colliders%rb(:,:) = colliders%rb(:,:) / frag%dscale
          colliders%vb(:,:) = colliders%vb(:,:) / frag%vscale
          colliders%mass(:) = colliders%mass(:) / frag%mscale
@@ -278,6 +282,7 @@ module subroutine fraggle_set_original_scale_factors(self, colliders)
          ! Restore scale factors
          frag%rbcom(:) = frag%rbcom(:) * frag%dscale
          frag%vbcom(:) = frag%vbcom(:) * frag%vscale
+         frag%rbimp(:) = frag%rbimp(:) * frag%dscale
 
          colliders%mass = colliders%mass * frag%mscale
          colliders%radius = colliders%radius * frag%dscale

src/fraggle/fraggle_util.f90

@@ -241,15 +241,15 @@ module subroutine fraggle_util_get_energy_momentum(self, colliders, system, para
             frag%ke_orbit_before = tmpsys%ke_orbit
             frag%ke_spin_before = tmpsys%ke_spin
             frag%pe_before = tmpsys%pe
-            frag%Etot_before = tmpsys%te
+            frag%Etot_before = tmpsys%te 
          else
             frag%Lorbit_after(:) = tmpsys%Lorbit(:)
             frag%Lspin_after(:) = tmpsys%Lspin(:)
             frag%Ltot_after(:) = tmpsys%Ltot(:)
             frag%ke_orbit_after = tmpsys%ke_orbit
             frag%ke_spin_after = tmpsys%ke_spin
             frag%pe_after = tmpsys%pe
-            frag%Etot_after = tmpsys%te
+            frag%Etot_after = tmpsys%te - (frag%pe_after - frag%pe_before) ! Gotta be careful with PE when number of bodies changes.
          end if
       end associate
 

src/modules/fraggle_classes.f90

@@ -54,6 +54,7 @@ module fraggle_classes
       ! Values in a coordinate frame centered on the collider barycenter and collisional system unit vectors (these are used internally by the fragment generation subroutine)
       real(DP), dimension(NDIM)              :: rbcom       !! Center of mass position vector of the collider system in system barycentric coordinates
       real(DP), dimension(NDIM)              :: vbcom       !! Velocity vector of the center of mass of the collider system in system barycentric coordinates
+      real(DP), dimension(NDIM)              :: rbimp       !! Impact point position vector of the collider system in system barycentric coordinates
       real(DP), dimension(NDIM)              :: x_coll_unit !! x-direction unit vector of collisional system
       real(DP), dimension(NDIM)              :: y_coll_unit !! y-direction unit vector of collisional system
       real(DP), dimension(NDIM)              :: z_coll_unit !! z-direction unit vector of collisional system

src/symba/symba_collision.f90

@@ -28,6 +28,7 @@ module function symba_collision_casedisruption(system, param, t)  result(status)
       integer(I4B)          :: i, ibiggest, nfrag
       logical               :: lfailure
       character(len=STRMAX) :: message 
+      real(DP) :: dpe
 
       associate(colliders => system%colliders, fragments => system%fragments)
 
@@ -46,6 +47,10 @@ module function symba_collision_casedisruption(system, param, t)  result(status)
          ! Generate the position and velocity distributions of the fragments
          call fragments%generate_fragments(colliders, system, param, lfailure)
 
+         dpe = fragments%pe_after - fragments%pe_before 
+         system%Ecollisions = system%Ecollisions - dpe 
+         system%Euntracked  = system%Euntracked + dpe 
+
          if (lfailure) then
             call io_log_one_message(FRAGGLE_LOG_OUT, "No fragment solution found, so treat as a pure hit-and-run")
             status = ACTIVE 
@@ -98,6 +103,7 @@ module function symba_collision_casehitandrun(system, param, t)  result(status)
       integer(I4B)                            :: i, ibiggest, nfrag, jtarg, jproj
       logical                                 :: lpure 
       character(len=STRMAX) :: message
+      real(DP) :: dpe
 
       associate(colliders => system%colliders, fragments => system%fragments)
          message = "Hit and run between"
@@ -123,6 +129,10 @@ module function symba_collision_casehitandrun(system, param, t)  result(status)
             ! Generate the position and velocity distributions of the fragments
             call fragments%generate_fragments(colliders, system, param, lpure)
 
+            dpe = fragments%pe_after - fragments%pe_before 
+            system%Ecollisions = system%Ecollisions - dpe 
+            system%Euntracked  = system%Euntracked + dpe 
+
             if (lpure) then
                call io_log_one_message(FRAGGLE_LOG_OUT, "Should have been a pure hit and run instead")
                nfrag = 0
@@ -172,8 +182,8 @@ module function symba_collision_casemerge(system, param, t)  result(status)
       integer(I4B)                             :: status    !! Status flag assigned to this outcome
       ! Internals
       integer(I4B)                              :: i, j, k, ibiggest
-      real(DP)                                  :: pe
       real(DP), dimension(NDIM)                 :: L_spin_new
+      real(DP)                                  :: dpe
       character(len=STRMAX) :: message
 
       associate(colliders => system%colliders, fragments => system%fragments)
@@ -207,9 +217,9 @@ module function symba_collision_casemerge(system, param, t)  result(status)
             ! Keep track of the component of potential energy due to the pre-impact colliders%idx for book-keeping
             ! Get the energy of the system after the collision
             call fragments%get_energy_and_momentum(colliders, system, param, lbefore=.false.)
-            pe = fragments%pe_after - fragments%pe_before 
-            system%Ecollisions = system%Ecollisions - pe 
-            system%Euntracked  = system%Euntracked + pe 
+            dpe = fragments%pe_after - fragments%pe_before 
+            system%Ecollisions = system%Ecollisions - dpe 
+            system%Euntracked  = system%Euntracked + dpe 
 
 
             ! Update any encounter lists that have the removed bodies in them so that they instead point to the new