From 609ee2ed90e9072ef78b856b5280f251f35e0d79 Mon Sep 17 00:00:00 2001 From: David Minton Date: Mon, 10 May 2021 15:56:53 -0400 Subject: [PATCH] Fixed angular momentum constraint --- src/symba/symba_frag_pos.f90 | 40 +++++++++++++++++------------------- 1 file changed, 19 insertions(+), 21 deletions(-) diff --git a/src/symba/symba_frag_pos.f90 b/src/symba/symba_frag_pos.f90 index a51813ba7..e554f3e58 100644 --- a/src/symba/symba_frag_pos.f90 +++ b/src/symba/symba_frag_pos.f90 @@ -125,6 +125,13 @@ subroutine symba_frag_pos (param, symba_plA, family, x, v, L_spin, Ip, mass, rad vb_frag(:,i) = v_frag(:, i) + vcom(:) end do + ke_target = 0._DP + do i = 1, nfrag + ke_target = ke_target + 0.5_DP * m_frag(i) * dot_product(vb_frag(:, i), vb_frag(:, i)) + end do + write(*, "(' ---------------------------------------------------------------------------')") + write(*,fmtlabel) ' T_frag calc |',ke_target / abs(Etot_before) + ! Calculate the new energy of the system of fragments call symba_frag_pos_energy_calc(npl, symba_plA, lexclude, ke_after, ke_spin_after, pe_after, Ltot_after,& nfrag=nfrag, Ip_frag=Ip_frag, m_frag=m_frag, rad_frag=rad_frag, xb_frag=xb_frag, vb_frag=vb_frag, rot_frag=rot_frag) @@ -250,7 +257,7 @@ subroutine symba_frag_pos_ang_mtm(nfrag, xcom, vcom, L_orb, L_spin, m_frag, rad_ real(DP) :: rmag, L_orb_frag_mag real(DP), dimension(NDIM) :: xc, vc, v_t_unit, h_unit, v_r_unit real(DP), dimension(NDIM) :: L_orb_old, L_spin_frag, L_spin_new - real(DP), parameter :: f_spin = 0.25_DP !! Fraction of pre-impact orbital angular momentum that is converted to fragment spin + real(DP), parameter :: f_spin = 0.00_DP !! Fraction of pre-impact orbital angular momentum that is converted to fragment spin L_orb_old(:) = L_orb(:, 1) + L_orb(:, 2) @@ -295,7 +302,7 @@ subroutine symba_frag_pos_kinetic_energy(xcom, vcom, L_orb, L_spin, m_frag, x_fr logical, intent(out) :: lmerge ! Internals - real(DP) :: v_corrected, mtot, A, C, rterm + real(DP) :: v_r_mag, mtot, ke_frag 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 @@ -307,37 +314,28 @@ subroutine symba_frag_pos_kinetic_energy(xcom, vcom, L_orb, L_spin, m_frag, x_fr allocate(v_t, mold=v_frag) call symba_frag_pos_com_adjust(xcom, m_frag, x_frag) - h_unit(:) = L_orb(:, 1) + L_orb(:, 2) + L_spin(:, 1) + L_spin(:, 2) - h_unit(:) = h_unit(:) / norm2(h_unit(:)) - do i = 1, nfrag v_r_unit(:, i) = x_frag(:,i) / norm2(x_frag(:, i)) - call utiL_crossproduct(h_unit(:), v_r_unit(:, i), v_t_unit(:)) - v_t(:,i) = dot_product(v_frag(:,i), v_t_unit(:)) * v_t_unit(:) + v_t(:,i) = v_frag(:, i) - dot_product(v_frag(:,i), v_r_unit(:, i)) * v_r_unit(:, i) end do - C = 2 * ke_target - A = mtot + ke_frag = ke_target - 0.5_DP * mtot * dot_product(vcom(:), vcom(:)) do i = 1, nfrag - C = C - m_frag(i) * (dot_product(vcom(:), vcom(:)) + dot_product(v_t(:,i), v_t(:,i)) + dot_product(vcom(:), v_t(:, i))) + ke_frag = ke_frag - 0.5_DP * m_frag(i) * dot_product(v_t(:, i), v_t(:, i)) end do - if (C > 0.0_DP) then - v_corrected = sqrt(C / A) + if (ke_frag > 0.0_DP) then + ke_frag = ke_frag / nfrag ! Equipartition the energy between all fragments + do i = 1, nfrag + v_r_mag = sqrt(2 * ke_frag / m_frag(i)) + v_frag(:, i) = v_r_mag * v_r_unit(:, i) + v_t(:, i) + end do lmerge = .false. else - v_corrected = 0.0_DP + v_frag(:, :) = 0.0_DP lmerge = .true. end if - ! Alter the fragment velocity distribution to hve the new radial component required to constrain kinetic energy - do i = 1, nfrag - v_frag(:,i) = v_corrected * v_r_unit(:, i) + v_t(:, i) - end do - - ! Shift the fragments into the system barycenter frame - call symba_frag_pos_com_adjust(vcom, m_frag, v_frag) - return end subroutine symba_frag_pos_kinetic_energy