From ed89188aad7447bcdffa58e1a714e23062b37449 Mon Sep 17 00:00:00 2001 From: David Minton Date: Wed, 12 May 2021 13:52:45 -0400 Subject: [PATCH] Fixed problem with the f_vectors accumulating each step of the secant method rather than being recomputed. --- src/symba/symba_frag_pos.f90 | 20 ++++++++++---------- 1 file changed, 10 insertions(+), 10 deletions(-) diff --git a/src/symba/symba_frag_pos.f90 b/src/symba/symba_frag_pos.f90 index 92d177456..3fc24a667 100644 --- a/src/symba/symba_frag_pos.f90 +++ b/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)