Added missing 0.5 term to kinetic energy equations

MintonGroup · May 12, 2021 · 419dcc7 · 419dcc7
1 parent 9fdfa1f
commit 419dcc7
Showing 1 changed file with 6 additions and 4 deletions.
diff --git a/src/symba/symba_frag_pos.f90 b/src/symba/symba_frag_pos.f90
@@ -398,7 +398,7 @@ function symba_frag_pos_fragment_velocity(nfrag, m_frag, x_frag, v_r_unit, Lambd
       real(DP)                              :: Beta                     !! Sum of the radial kinetic energy of i>4 fragments
       real(DP), dimension(4)                :: v_r_mag_01, v_r_mag_02   !! Two initial value guesses for the radial velocity magnitude of the first four fragments
       real(DP), dimension(4)                :: f_vec_01, f_vec_02, f_vec_const   !! Equation vectors for initial value guesses 1 and 2
-      integer(I4B), parameter               :: MAXITER = 50 
+      integer(I4B), parameter               :: MAXITER = 500 
       real(DP), parameter                   :: TOL = 1e-8_DP
       real(DP)                              :: vmult
 
@@ -414,10 +414,11 @@ function symba_frag_pos_fragment_velocity(nfrag, m_frag, x_frag, v_r_unit, Lambd
          Beta = 0.0_DP
          do i = 5, nfrag
             Gam = Gam + m_frag(i) * v_r_mag(i) * v_r_unit(:,i)
-            Beta = Beta + m_frag(i) * v_r_mag(i)**2
+            Beta = Beta + 0.5_DP * m_frag(i) * v_r_mag(i)**2
          end do
          if (Beta - Lambda > 0.0_DP) then
             vmult = 0.8_DP * vmult
+            write(*,*) 'preconditioning: ',vmult
          else
             exit
          end if
@@ -437,10 +438,10 @@ function symba_frag_pos_fragment_velocity(nfrag, m_frag, x_frag, v_r_unit, Lambd
          do i = 1, 4
             if (j == 1) then
                f_vec_01(1:3) = f_vec_01(1:3) + m_frag(i) * v_r_mag_01(i) * v_r_unit(:, i)
-               f_vec_01(4)   = f_vec_01(4)   + m_frag(i) * v_r_mag_01(i)**2
+               f_vec_01(4)   = f_vec_01(4)   + 0.5_DP * m_frag(i) * v_r_mag_01(i)**2
             end if
             f_vec_02(1:3) = f_vec_02(1:3) + m_frag(i) * v_r_mag_02(i) * v_r_unit(:,i)
-            f_vec_02(4)   = f_vec_02(4)   + m_frag(i) * v_r_mag_02(i)**2
+            f_vec_02(4)   = f_vec_02(4)   + 0.5_DP * m_frag(i) * v_r_mag_02(i)**2
          end do
 
          write(*,*) 'Iteration: ',j
@@ -453,6 +454,7 @@ function symba_frag_pos_fragment_velocity(nfrag, m_frag, x_frag, v_r_unit, Lambd
          do i = 1, 4
             v_r_mag(i) = v_r_mag_02(i) - f_vec_02(i) * (v_r_mag_02(i) - v_r_mag_01(i)) / (f_vec_02(i) - f_vec_01(i))
          end do
+         v_r_mag(:) = abs(v_r_mag(:))
 
          if (norm2(f_vec_02(:)) < TOL) exit