From ff898874a530930415a7120b7bb9776b7536baea Mon Sep 17 00:00:00 2001 From: David A Minton Date: Mon, 20 Feb 2023 13:23:46 -0500 Subject: [PATCH] Fixed issue of convergence on mass multiplier for SFD --- src/fraggle/fraggle_util.f90 | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/src/fraggle/fraggle_util.f90 b/src/fraggle/fraggle_util.f90 index a0dc08af5..846122e73 100644 --- a/src/fraggle/fraggle_util.f90 +++ b/src/fraggle/fraggle_util.f90 @@ -89,7 +89,7 @@ module subroutine fraggle_util_set_mass_dist(self, param) integer(I4B) :: i, j, jproj, jtarg, nfrag, istart, nfragmax, nrem real(DP), dimension(2) :: volume real(DP), dimension(NDIM) :: Ip_avg - real(DP) :: mfrag, mremaining, mtot, mcumul, G, mass_noise, Mslr, mscale, y, yp, Mrat, m0 + real(DP) :: mfrag, mremaining, mtot, mcumul, G, mass_noise, Mslr, mscale, mscale0, y, yp, Mrat real(DP), dimension(:), allocatable :: mass real(DP) :: beta integer(I4B), parameter :: MASS_NOISE_FACTOR = 5 !! The number of digits of random noise that get added to the minimum mass value to prevent identical masses from being generated in a single run @@ -177,13 +177,13 @@ module subroutine fraggle_util_set_mass_dist(self, param) ! Use Newton's method solver to get the logspace slope of the mass function Mrat = mremaining / min_mfrag nrem = nfrag - 2 - mscale = 1.0_DP - 1.0_DP/Mrat + mscale = Mrat**(1.0_DP/nrem) do j = 1, MAXLOOP y = Mrat - (1.0_DP - mscale**nrem)/(1.0_DP - mscale) yp = (mscale + mscale**nrem * ((nrem - 1) * mscale - nrem)) / (mscale * (mscale - 1.0_DP)**2) - m0 = mscale - mscale = m0 + y/yp - if (abs((mscale - m0)/mscale) < epsilon(1.0_DP)) exit + mscale0 = mscale + mscale = mscale0 + y/yp + if (abs((mscale - mscale0)/mscale) < epsilon(1.0_DP)) exit end do Mslr = impactors%mass_dist(iMslr)