function guess_target_nref() result(target_nref)
! try to extrapolate how many SIs are required for a given target population
use FciMCData, only: AllNoatHF
implicit none
integer :: target_nref
real(dp) :: alpha, beta, cAllNoatHF
target_nref = nRefs
cAllNoatHF = sum(AllNoatHF) / inum_runs
! if the reference population did not change, nothing to do
if (lastAllNoatHF.isclose.cAllNoatHF) return
! if we did not change the SI threshold, nothing to do
if (lastNRefs == nRefs) return
beta = log(lastAllNoatHF / cAllNoatHF) / (lastNRefs - nRefs)
alpha = lastAllNoatHF * exp(log(lastAllNoatHF / cAllNoatHF) * (lastNRefs / (lastNRefs - nRefs)))
! here, we really assume exponential behaviour and extrapolate the nRefs that would then
! give the reference population closest to the target
target_nref = int(log(targetRefPop / alpha) / beta)
end function guess_target_nref