Fixed problem with the f_vectors accumulating each step of the secant method rather than being recomputed.

src/symba/symba_frag_pos.f90

@@ -397,7 +397,7 @@ function symba_frag_pos_fragment_velocity(nfrag, m_frag, x_frag, v_r_unit, Lambd
       real(DP), dimension(NDIM)             :: Gam                      !! Sum of the radial momentum vector of i>4 fragments
       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       !! Equation vectors for initial value guesses 1 and 2
+      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 
       real(DP), parameter                   :: TOL = 1e-8_DP
 
@@ -420,17 +420,18 @@ function symba_frag_pos_fragment_velocity(nfrag, m_frag, x_frag, v_r_unit, Lambd
          Beta = Beta + m_frag(i) * v_r_mag(i)**2
       end do
 
-      f_vec_01(1:3) = Gam(:) + tau(:)
-      f_vec_01(4)   = Beta - Lambda 
-      f_vec_02(:)   = f_vec_01(:)
+      f_vec_const(1:3) = Gam(:) + tau(:)
+      f_vec_const(4)   = Beta - Lambda 
 
-      do i = 1, 4
-         f_vec_01(1:3) = f_vec_01(1:3) + m_frag(i) * v_r_mag_01(i) * x_frag(:,i)
-         f_vec_01(4)   = f_vec_01(4)   + m_frag(i) * v_r_mag_01(i)**2
-      end do
+      f_vec_01(:) = f_vec_const(:)
 
       do j = 1, MAXITER
+         f_vec_02(:) = f_vec_const(:)
          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) * x_frag(:,i)
+               f_vec_01(4)   = f_vec_01(4)   + 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) * x_frag(:,i)
             f_vec_02(4)   = f_vec_02(4)   + m_frag(i) * v_r_mag_02(i)**2
          end do
@@ -440,6 +441,7 @@ function symba_frag_pos_fragment_velocity(nfrag, m_frag, x_frag, v_r_unit, Lambd
          write(*,*) 'v_r_mag_02: ',v_r_mag_02(:)
          write(*,*) 'f_vec_01:   ',f_vec_01(:)
          write(*,*) 'f_vec_02:   ',f_vec_02(:)
+         write(*,*) 'norm2(f):   ',norm2(f_vec_02(:))
 
          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))
@@ -448,8 +450,6 @@ function symba_frag_pos_fragment_velocity(nfrag, m_frag, x_frag, v_r_unit, Lambd
          if (norm2(f_vec_02(:)) < TOL) exit
 
          write(*,*) 'v_r_mag:    ',v_r_mag(1:4)
-         write(*,*) 'KE_after:   ',KE_after
-         write(*,*) 'Lambda:     ',Lambda
 
          v_r_mag_01(:) = v_r_mag_02(:)
          v_r_mag_02(:) = v_r_mag(1:4)