diff --git a/src/collision/collision_util.f90 b/src/collision/collision_util.f90 index 2e63029ce..0c435cfc6 100644 --- a/src/collision/collision_util.f90 +++ b/src/collision/collision_util.f90 @@ -347,34 +347,34 @@ module subroutine collision_util_set_coordinate_system(self) ! and the y-axis aligned with the pre-impact distance vector. ! y-axis is the separation distance - fragments%y_coll_unit(:) = .unit.delta_r(:) + impactors%y_unit(:) = .unit.delta_r(:) Ltot = impactors%Lorbit(:,1) + impactors%Lorbit(:,2) + impactors%Lspin(:,1) + impactors%Lspin(:,2) L_mag = .mag.Ltot(:) if (L_mag > sqrt(tiny(L_mag))) then - fragments%z_coll_unit(:) = .unit.Ltot(:) + impactors%z_unit(:) = .unit.Ltot(:) else ! Not enough angular momentum to determine a z-axis direction. We'll just pick a random direction - call random_number(fragments%z_coll_unit(:)) - fragments%z_coll_unit(:) = .unit.fragments%z_coll_unit(:) + call random_number(impactors%z_unit(:)) + impactors%z_unit(:) = .unit.impactors%z_unit(:) end if ! The cross product of the y- by z-axis will give us the x-axis - fragments%x_coll_unit(:) = fragments%y_coll_unit(:) .cross. fragments%z_coll_unit(:) + impactors%x_unit(:) = impactors%y_unit(:) .cross. impactors%z_unit(:) - fragments%v_coll_unit(:) = .unit.delta_v(:) + impactors%v_unit(:) = .unit.delta_v(:) - if (.not.any(fragments%r_coll(:,:) > 0.0_DP)) return - fragments%rmag(:) = .mag. fragments%r_coll(:,:) + if (.not.any(fragments%rc(:,:) > 0.0_DP)) return + fragments%rmag(:) = .mag. fragments%rc(:,:) ! 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 call random_number(L_sigma(:,:)) do concurrent(i = 1:nfrag, fragments%rmag(i) > 0.0_DP) - fragments%v_n_unit(:, i) = fragments%z_coll_unit(:) + 2e-1_DP * (L_sigma(:,i) - 0.5_DP) + fragments%v_n_unit(:, i) = impactors%z_unit(:) + 2e-1_DP * (L_sigma(:,i) - 0.5_DP) end do ! Define the radial, normal, and tangential unit vectors for each individual fragment - fragments%v_r_unit(:,:) = .unit. fragments%r_coll(:,:) + fragments%v_r_unit(:,:) = .unit. fragments%rc(:,:) fragments%v_n_unit(:,:) = .unit. fragments%v_n_unit(:,:) fragments%v_t_unit(:,:) = .unit. (fragments%v_n_unit(:,:) .cross. fragments%v_r_unit(:,:)) diff --git a/src/fraggle/fraggle_generate.f90 b/src/fraggle/fraggle_generate.f90 index e603b5b9d..ffc7acf7a 100644 --- a/src/fraggle/fraggle_generate.f90 +++ b/src/fraggle/fraggle_generate.f90 @@ -213,35 +213,35 @@ subroutine fraggle_generate_pos_vec(fragments, impactors, r_max_start) ! 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(fragments%r_coll(:,3:nfrag)) + call random_number(fragments%rc(:,3:nfrag)) loverlap(:) = .true. do while (any(loverlap(3:nfrag))) - fragments%r_coll(:, 1) = impactors%rb(:, 1) - fragments%rbcom(:) - fragments%r_coll(:, 2) = impactors%rb(:, 2) - fragments%rbcom(:) + fragments%rc(:, 1) = impactors%rb(:, 1) - fragments%rbcom(:) + fragments%rc(:, 2) = impactors%rb(:, 2) - fragments%rbcom(:) r_max = r_max + 0.1_DP * rad do i = 3, nfrag if (loverlap(i)) then - call random_number(fragments%r_coll(:,i)) - fragments%r_coll(:,i) = 2 * (fragments%r_coll(:, i) - 0.5_DP) - fragments%r_coll(:, i) = fragments%r_coll(:,i) + fdistort * vunit(:) - fragments%r_coll(:, i) = r_max * fragments%r_coll(:, i) - fragments%r_coll(:, i) = fragments%r_coll(:, i) + (fragments%rbimp(:) - fragments%rbcom(:)) ! Shift the center of the fragment cloud to the impact point rather than the CoM - !if (lnoncat .and. dot_product(fragments%r_coll(:,i), runit(:)) < 0.0_DP) fragments%r_coll(:, i) = -fragments%r_coll(:, i) ! Make sure the fragment cloud points away from the impact point + 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) + (fragments%rbimp(:) - fragments%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 end if end do loverlap(:) = .false. do j = 1, nfrag do i = j + 1, nfrag - dis = .mag.(fragments%r_coll(:,j) - fragments%r_coll(:,i)) + dis = .mag.(fragments%rc(:,j) - fragments%rc(:,i)) loverlap(i) = loverlap(i) .or. (dis <= (fragments%radius(i) + fragments%radius(j))) end do end do end do - call fraggle_util_shift_vector_to_origin(fragments%mass, fragments%r_coll) + call fraggle_util_shift_vector_to_origin(fragments%mass, fragments%rc) call fragments%set_coordinate_system(impactors) do concurrent(i = 1:nfrag) - fragments%rb(:,i) = fragments%r_coll(:,i) + fragments%rbcom(:) + fragments%rb(:,i) = fragments%rc(:,i) + fragments%rbcom(:) end do fragments%rbcom(:) = 0.0_DP @@ -408,7 +408,7 @@ subroutine fraggle_generate_tan_vel(fragments, impactors, lfailure) fragments%vb(:,1:nfrag) = fraggle_util_vmag_to_vb(fragments%v_r_mag(1:nfrag), fragments%v_r_unit(:,1:nfrag), fragments%v_t_mag(1:nfrag), & fragments%v_t_unit(:,1:nfrag), fragments%mass(1:nfrag), fragments%vbcom(:)) do concurrent (i = 1:nfrag) - fragments%v_coll(:,i) = fragments%vb(:,i) - fragments%vbcom(:) + fragments%vc(:,i) = fragments%vb(:,i) - fragments%vbcom(:) end do ! Now do a kinetic energy budget check to make sure we are still within the budget. @@ -422,7 +422,7 @@ subroutine fraggle_generate_tan_vel(fragments, impactors, lfailure) ke_diff = fragments%ke_budget - fragments%ke_spin - fragments%ke_orbit lfailure = ke_diff < 0.0_DP if (.not.lfailure) exit - fragments%r_coll(:,:) = fragments%r_coll(:,:) * 1.1_DP + fragments%rc(:,:) = fragments%rc(:,:) * 1.1_DP end do if (lfailure) then call io_log_one_message(FRAGGLE_LOG_OUT, " ") @@ -476,7 +476,7 @@ function solve_fragment_tan_vel(lfailure, v_t_mag_input) result(v_t_mag_output) else if (present(v_t_mag_input)) then vtmp(:) = v_t_mag_input(i - 6) * fragments%v_t_unit(:, i) L_lin_others(:) = L_lin_others(:) + fragments%mass(i) * vtmp(:) - L(:) = fragments%mass(i) * (fragments%r_coll(:, i) .cross. vtmp(:)) + L(:) = fragments%mass(i) * (fragments%rc(:, i) .cross. vtmp(:)) L_orb_others(:) = L_orb_others(:) + L(:) end if end do @@ -582,7 +582,7 @@ subroutine fraggle_generate_rad_vel(fragments, impactors, lfailure) fragments%vb(:,1:nfrag) = fraggle_util_vmag_to_vb(fragments%v_r_mag(1:nfrag), fragments%v_r_unit(:,1:nfrag), & fragments%v_t_mag(1:nfrag), fragments%v_t_unit(:,1:nfrag), fragments%mass(1:nfrag), fragments%vbcom(:)) do i = 1, nfrag - fragments%v_coll(:, i) = fragments%vb(:, i) - fragments%vbcom(:) + fragments%vc(:, i) = fragments%vb(:, i) - fragments%vbcom(:) fragments%ke_orbit = fragments%ke_orbit + fragments%mass(i) * norm2(fragments%vb(:, i)) fragments%ke_spin = fragments%ke_spin + fragments%mass(i) * fragments%radius(i)**2 * fragments%Ip(3,i) * norm2(fragments%rot(:,i)) end do diff --git a/src/fraggle/fraggle_set.f90 b/src/fraggle/fraggle_set.f90 index a69a69260..0da950206 100644 --- a/src/fraggle/fraggle_set.f90 +++ b/src/fraggle/fraggle_set.f90 @@ -244,12 +244,12 @@ module subroutine fraggle_set_original_scale_factors(self) fragments%mass = fragments%mass * fragments%mscale fragments%radius = fragments%radius * fragments%dscale fragments%rot = fragments%rot / fragments%tscale - fragments%r_coll = fragments%r_coll * fragments%dscale - fragments%v_coll = fragments%v_coll * fragments%vscale + fragments%rc = fragments%rc * fragments%dscale + fragments%vc = fragments%vc * fragments%vscale do i = 1, fragments%nbody - fragments%rb(:, i) = fragments%r_coll(:, i) + fragments%rbcom(:) - fragments%vb(:, i) = fragments%v_coll(:, i) + fragments%vbcom(:) + fragments%rb(:, i) = fragments%rc(:, i) + fragments%rbcom(:) + fragments%vb(:, i) = fragments%vc(:, i) + fragments%vbcom(:) end do impactors%Qloss = impactors%Qloss * fragments%Escale