From d2f426c6588b8c1ebd87c03fe740aab4134024b7 Mon Sep 17 00:00:00 2001 From: David Minton Date: Tue, 11 May 2021 14:34:46 -0400 Subject: [PATCH] Refactored to make variables names consistent with analytical equations --- src/symba/symba_frag_pos.f90 | 13 ++++++++----- 1 file changed, 8 insertions(+), 5 deletions(-) diff --git a/src/symba/symba_frag_pos.f90 b/src/symba/symba_frag_pos.f90 index 44e44063e..9a2664a7f 100644 --- a/src/symba/symba_frag_pos.f90 +++ b/src/symba/symba_frag_pos.f90 @@ -332,7 +332,10 @@ subroutine symba_frag_pos_kinetic_energy(xcom, vcom, L_orb, L_spin, m_frag, x_fr logical, intent(out) :: lmerge ! Internals - real(DP) :: mtot, ke_frag + real(DP) :: mtot !! Total mass of fragments + real(DP) :: Lambda !! Sum of the radial kinetic energy of all fragments + real(DP) :: Gam !! Sum of the radial momentum vector of i>4 fragments + real(DP), dimension(NDIM) :: tau !! Sum of the tangential momentum vector of all fragments integer(I4B) :: i, nfrag real(DP), dimension(:,:), allocatable :: v_r_unit, v_t real(DP), dimension(NDIM) :: v_t_unit, h_unit, L_orb_frag @@ -353,16 +356,16 @@ subroutine symba_frag_pos_kinetic_energy(xcom, vcom, L_orb, L_spin, m_frag, x_fr v_t(:,i) = v_frag(:, i) - dot_product(v_frag(:,i), v_r_unit(:, i)) * v_r_unit(:, i) end do - ke_frag = ke_target - 0.5_DP * mtot * dot_product(vcom(:), vcom(:)) + Lambda = ke_target - 0.5_DP * mtot * dot_product(vcom(:), vcom(:)) do i = 1, nfrag - ke_frag = ke_frag - 0.5_DP * m_frag(i) * dot_product(v_t(:, i), v_t(:, i)) + Lambda = Lambda - 0.5_DP * m_frag(i) * dot_product(v_t(:, i), v_t(:, i)) end do - if (ke_frag > 0.0_DP) then + if (Lambda > 0.0_DP) then lmerge = .false. call random_number(v_r_mag(:)) ! Our initial guess for the first 4 fragments and the values of will be based on an equipartition of KE with some random variation ! We shift the random variate to the range 0.5, 1.5 to prevent any zero values for radial velocities - v_r_mag(:) = sqrt(2 * ke_frag / nfrag / m_frag(:)) * (v_r_mag(:) + 0.5_DP) + v_r_mag(:) = sqrt(2 * Lambda / nfrag / m_frag(:)) * (v_r_mag(:) + 0.5_DP) do i = 1, nfrag v_frag(:, i) = v_r_mag(i) * v_r_unit(:, i) + v_t(:, i) end do