diff --git a/src/fraggle/fraggle_generate.f90 b/src/fraggle/fraggle_generate.f90 index e555b11f2..8e52da0c5 100644 --- a/src/fraggle/fraggle_generate.f90 +++ b/src/fraggle/fraggle_generate.f90 @@ -227,7 +227,6 @@ module subroutine fraggle_generate_fragments(collider, nbody_system, param, lfai logical, dimension(size(IEEE_ALL)) :: fpe_halting_modes, fpe_quiet_modes logical, dimension(size(IEEE_USUAL)) :: fpe_flag character(len=STRMAX) :: message - real(DP), dimension(NDIM) :: runit, vunit ! The minimization and linear solvers can sometimes lead to floating point exceptions. Rather than halting the code entirely if this occurs, we ! can simply fail the attempt and try again. So we need to turn off any floating point exception halting modes temporarily @@ -267,9 +266,10 @@ module subroutine fraggle_generate_fragments(collider, nbody_system, param, lfai call collider%get_energy_and_momentum(nbody_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.2_DP * .mag.(impactors%rb(:,2) - impactors%rb(:,1)) + r_max_start = 1.1_DP * .mag.(impactors%rb(:,2) - impactors%rb(:,1)) lfailure = .false. try = 1 + f_spin = F_SPIN_FIRST do while (try < MAXTRY) write(message,*) try call swiftest_io_log_one_message(COLLISION_LOG_OUT, "Fraggle try " // trim(adjustl(message))) @@ -285,9 +285,6 @@ module subroutine fraggle_generate_fragments(collider, nbody_system, param, lfai call fraggle_generate_pos_vec(collider, r_max_start) call collider%set_coordinate_system() - ! This is a factor that will "distort" the shape of the fragment cloud in the direction of the impact velocity - f_spin= .mag. (impactors%y_unit(:) .cross. impactors%v_unit(:)) - ! Initial velocity guess will be the barycentric velocity of the colliding nbody_system so that the budgets are based on the much smaller collisional-frame velocities do concurrent (i = 1:nfrag) fragments%vb(:, i) = impactors%vbcom(:) @@ -296,12 +293,6 @@ module subroutine fraggle_generate_fragments(collider, nbody_system, param, lfai call collider%get_energy_and_momentum(nbody_system, param, lbefore=.false.) call collider%set_budgets() - call fraggle_generate_spins(collider, f_spin, lfailure) - if (lfailure) then - call swiftest_io_log_one_message(COLLISION_LOG_OUT, "Fraggle failed to find spins") - cycle - end if - call fraggle_generate_tan_vel(collider, lfailure) if (lfailure) then call swiftest_io_log_one_message(COLLISION_LOG_OUT, "Fraggle failed to find tangential velocities") @@ -314,6 +305,12 @@ module subroutine fraggle_generate_fragments(collider, nbody_system, param, lfai cycle end if + call fraggle_generate_spins(collider, f_spin, lfailure) + if (lfailure) then + call swiftest_io_log_one_message(COLLISION_LOG_OUT, "Fraggle failed to find spins") + cycle + end if + call collider%get_energy_and_momentum(nbody_system, param, lbefore=.false.) dEtot = collider%Etot(2) - collider%Etot(1) dLmag = .mag. (collider%Ltot(:,2) - collider%Ltot(:,1)) @@ -378,7 +375,6 @@ subroutine fraggle_generate_pos_vec(collider, 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, fdistort - real(DP), dimension(NDIM) :: runit, vunit logical, dimension(:), allocatable :: loverlap integer(I4B) :: i, j logical :: lnoncat, lhitandrun @@ -394,35 +390,32 @@ subroutine fraggle_generate_pos_vec(collider, r_max_start) r_max = r_max_start rad = sum(impactors%radius(:)) + ! This is a factor that will "distort" the shape of the fragment cloud in the direction of the impact velocity + fdistort = .mag. (impactors%y_unit(:) .cross. impactors%v_unit(:)) - ! Get the unit vectors for the relative position and velocity vectors. These are used to shift the fragment cloud depending on the - runit(:) = impactors%rb(:,2) - impactors%rb(:,1) - runit(:) = runit(:) / (.mag. runit(:)) - - vunit(:) = impactors%vb(:,2) - impactors%vb(:,1) - vunit(:) = vunit(:) / (.mag. vunit(:)) - - ! This is a factor that will "distort" the shape of the frgment cloud in the direction of the impact velocity - fdistort = .mag. (runit(:) .cross. vunit(:)) - - ! 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. + ! We will treat the first two fragments of the list as special cases. They get initialized at the original positions of the impactor bodies + fragments%rc(:, 1) = impactors%rb(:, 1) - impactors%rbcom(:) + fragments%rc(:, 2) = impactors%rb(:, 2) - impactors%rbcom(:) call random_number(fragments%rc(:,3:nfrag)) loverlap(:) = .true. do while (any(loverlap(3:nfrag))) - fragments%rc(:, 1) = impactors%rb(:, 1) - impactors%rbcom(:) - fragments%rc(:, 2) = impactors%rb(:, 2) - impactors%rbcom(:) - r_max = r_max + 0.1_DP * rad + if (lhitandrun) then ! For a hit-and-run with disruption, the fragment cloud size is based on the radius of the disrupted body + r_max = 2 * impactors%radius(2) + else ! For disruptions, the the fragment cloud size is based on the mutual collision system + r_max = r_max + 0.1_DP * rad + end if do i = 3, nfrag - if (loverlap(i)) then + if (loverlap(i)) then call random_number(fragments%rc(:,i)) fragments%rc(:,i) = 2 * (fragments%rc(:, i) - 0.5_DP) - fragments%rc(:, i) = fragments%rc(:,i) + fdistort * vunit(:) - fragments%rc(:, i) = r_max * fragments%rc(:, i) + fragments%rc(:, i) = fragments%rc(:,i) + fdistort * impactors%v_unit(:) + fragments%rc(:, i) = r_max * fragments%rc(:, i) fragments%rc(:, i) = fragments%rc(:, i) + (impactors%rbimp(:) - impactors%rbcom(:)) ! Shift the center of the fragment cloud to the impact point rather than the CoM - !if (lnoncat .and. dot_product(fragments%rc(:,i), runit(:)) < 0.0_DP) fragments%rc(:, i) = -fragments%rc(:, i) ! Make sure the fragment cloud points away from the impact point + if (lnoncat .and. dot_product(fragments%rc(:,i), impactors%y_unit(:)) < 0.0_DP) fragments%rc(:, i) = -fragments%rc(:, i) ! Make sure the fragment cloud points away from the impact point end if end do + + ! Check for any overlapping bodies. loverlap(:) = .false. do j = 1, nfrag do i = j + 1, nfrag