Set up simple disruption model

MintonGroup · Dec 23, 2022 · 3350420 · 3350420
1 parent ba80e3c
commit 3350420
Show file tree

Hide file tree

Showing 3 changed files with 20 additions and 8 deletions.
diff --git a/src/collision/collision_generate.f90 b/src/collision/collision_generate.f90
@@ -328,28 +328,38 @@ end subroutine collision_generate_simple_pos_vec
    module subroutine collision_generate_simple_vel_vec(collider)
       !! Author:  David A. Minton
       !!
-      !! Generates an initial "guess" for the velocitity vectors 
+      !! Generates an initial velocity distribution. For disruptions, the velocity magnitude is set to be
+      !! 2x the escape velocity of the colliding pair. For hit and runs the velocity magnitude is set to be
+      !! 2x the escape velocity of the smallest of the two bodies.
       implicit none
       ! Arguments
       class(collision_simple_disruption), intent(inout) :: collider !! Fraggle collision system object
       ! Internals
       integer(I4B) :: i, j
-      logical :: lcat
-      real(DP), dimension(NDIM) :: vimp_unit
+      logical :: lhr
+      real(DP), dimension(NDIM) :: vimp_unit, vcom, rnorm
       real(DP), dimension(NDIM,collider%fragments%nbody) :: vnoise
       real(DP), parameter :: VNOISE_MAG = 0.10_DP
       real(DP) :: vmag
 
       associate(fragments => collider%fragments, impactors => collider%impactors, nfrag => collider%fragments%nbody)
-         lcat = (impactors%regime == COLLRESOLVE_REGIME_SUPERCATASTROPHIC) 
+         lhr = (impactors%regime == COLLRESOLVE_REGIME_HIT_AND_RUN) 
 
          ! "Bounce" the first two bodies
          do i = 1,2
-            fragments%vc(:,i) = impactors%vb(:,i) - impactors%vbcom(:)  - 2 * impactors%vbimp(:)
+            vcom(:) = impactors%vb(:,i) - impactors%vbcom(:)
+            rnorm(:) = impactors%y_unit(:)
+            ! Do the reflection
+            vcom(:) = vcom(:) - 2 * dot_product(vcom(:),rnorm(:)) * rnorm(:)
+            fragments%vc(:,i) = vcom(:)
          end do
 
-         vmag = .mag.impactors%vbcom(:) 
-         vimp_unit(:) = .unit. impactors%vbimp(:)
+         ! Compute the escape velocity
+         if (lhr) then
+            vmag = 2 * sqrt(2 * impactors%Gmass(2) / impactors%radius(2))
+         else
+            vmag = 2 * sqrt(2 * sum(impactors%Gmass(:)) / (.mag. (impactors%rb(:,2) - impactors%rb(:,1))))
+         end if
          call random_number(vnoise)
          vnoise = (2 * vnoise - 1.0_DP) * vmag
          do i = 3, nfrag

diff --git a/src/collision/collision_module.f90 b/src/collision/collision_module.f90
@@ -58,6 +58,7 @@ module collision
       real(DP),     dimension(NDIM,2)              :: Lspin     !! Two-body equivalent spin angular momentum vectors of the collider bodies prior to collision
       real(DP),     dimension(NDIM,2)              :: Lorbit    !! Two-body equivalent orbital angular momentum vectors of the collider bodies prior to collision
       real(DP),     dimension(NDIM,2)              :: Ip        !! Two-body equivalent principal axes moments of inertia the collider bodies prior to collision
+      real(DP),     dimension(2)                   :: Gmass     !! Two-body equivalent G*mass of the collider bodies prior to the collision
       real(DP),     dimension(2)                   :: mass      !! Two-body equivalent mass of the collider bodies prior to the collision
       real(DP),     dimension(2)                   :: radius    !! Two-body equivalent radii of the collider bodies prior to the collision
       real(DP)                                     :: Qloss     !! Energy lost during the collision

diff --git a/src/collision/collision_resolve.f90 b/src/collision/collision_resolve.f90
@@ -57,7 +57,8 @@ module subroutine collision_resolve_consolidate_impactors(self, nbody_system, pa
                pl%kin(idx_parent(2))%parent = idx_parent(1)
             end if
 
-            impactors%mass(:) = pl%mass(idx_parent(:)) ! Note: This is meant to mass, not G*mass, as the collisional regime determination uses mass values that will be converted to Si
+            impactors%Gmass(:) = pl%Gmass(idx_parent(:)) 
+            impactors%mass(:)  = pl%mass(idx_parent(:)) 
             impactors%radius(:) = pl%radius(idx_parent(:))
             volume(:) =  (4.0_DP / 3.0_DP) * PI * impactors%radius(:)**3