From 2bc03bdcf511872e974bdd729ac98f9338eacb4a Mon Sep 17 00:00:00 2001 From: Carlisle Wishard Date: Wed, 12 May 2021 13:37:24 -0400 Subject: [PATCH] added tolerance value for radial velocity magnitudes, --- src/symba/symba_frag_pos.f90 | 25 +++++++++++++------------ 1 file changed, 13 insertions(+), 12 deletions(-) diff --git a/src/symba/symba_frag_pos.f90 b/src/symba/symba_frag_pos.f90 index 8dace48ff..e32130998 100644 --- a/src/symba/symba_frag_pos.f90 +++ b/src/symba/symba_frag_pos.f90 @@ -398,8 +398,8 @@ 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 !! Equation vectors for initial value guesses 1 and 2 - real(DP) :: KE_after !! Radial kinetic energy of the four fragments given their new velocities integer(I4B), parameter :: MAXITER = 50 + real(DP), parameter :: TOL = 1e-8_DP ! Initialize the fragment radial velocities with random values. The first 4 serve as guesses that get updated with the secant method solver. ! We shift the random variate to the range 0.5, 1.5 to prevent any zero values for radial velocities @@ -424,15 +424,17 @@ function symba_frag_pos_fragment_velocity(nfrag, m_frag, x_frag, v_r_unit, Lambd f_vec_01(4) = Beta - Lambda f_vec_02(:) = f_vec_01(:) + 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 + do j = 1, MAXITER 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 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 - KE_after = 0.0_DP write(*,*) 'Iteration: ',j write(*,*) 'v_r_mag_01: ',v_r_mag_01(:) write(*,*) 'v_r_mag_02: ',v_r_mag_02(:) @@ -441,19 +443,18 @@ 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)) - KE_after = KE_after + 0.5_DP * m_frag(i) * v_r_mag(i)**2 end do - write(*,*) 'v_r_mag: ',v_r_mag(1:4) + if (norm2(f_vec_2(:)) < TOL) exit + write(*,*) 'v_r_mag: ',v_r_mag(1:4) write(*,*) 'KE_after: ',KE_after write(*,*) 'Lambda: ',Lambda - if (KE_after <= Lambda) then - exit - else - v_r_mag_01(:) = v_r_mag_02(:) - v_r_mag_02(:) = v_r_mag(1:4) - end if + + v_r_mag_01(:) = v_r_mag_02(:) + v_r_mag_02(:) = v_r_mag(1:4) + f_vec_01(:) = f_vec_02(:) + end do return