Fix logic on CalcTempAtFalseVacFraction() (#149)

* Fix logic on "CalcTempAtFalseVacFraction()"

* Bump version

* Fix "magic number"

* Increase version one more since there was a fix before the PR being approved.

* Another version bump

CMakeLists.txt

@@ -6,7 +6,7 @@
 cmake_minimum_required(VERSION 3.23)
 project(
   BSMPT
-  VERSION 3.0.4
+  VERSION 3.0.5
   LANGUAGES C CXX
   DESCRIPTION
     "BSMPT - Beyond the Standard Model Phase Transitions : A C++ package for the computation of the EWPT in BSM models"

README.md

@@ -4,7 +4,7 @@ SPDX-FileCopyrightText: 2021 Philipp Basler, Margarete Mühlleitner and Jonas M
 SPDX-License-Identifier: GPL-3.0-or-later
 -->
 
-Program: BSMPT version 3.0.4
+Program: BSMPT version 3.0.5
 
 Released by: Philipp Basler, Lisa Biermann, Margarete Mühlleitner, Jonas Müller, Rui Santos and João Viana
 

include/BSMPT/bounce_solution/bounce_solution.h

@@ -225,6 +225,25 @@ class BounceSolution
    */
   const double RelErr = 1e-6;
 
+  /**
+   * @brief Maximum relative difference in temperature on the fraction of false
+   * vacuum to be accepted.
+   *
+   */
+  const double RelativeTemperatureInCalcTempAtFalseVacFraction = 1e-8;
+
+  /**
+   * @brief Maximum relative error on the fraction of vacuum tunneled to be
+   * accepted.
+   */
+  const double RelativeErrorInCalcTempAtFalseVacFraction = 1e-3;
+
+  /**
+   * @brief Additional margin of error in the while loop without admitting
+   * failure.
+   */
+  const double MarginOfCalcTempAtFalseVacFractionBeforeFailure = 1e-4;
+
   /**
    * @brief pair of coexisiting phases
    */
@@ -547,4 +566,4 @@ struct resultErrorPair Nintegrate_Outer(BounceSolution &obj);
  */
 struct resultErrorPair Nderive_BounceRatio(BounceSolution &obj);
 
-} // namespace BSMPT
+} // namespace BSMPT

src/bounce_solution/bounce_solution.cpp

@@ -754,9 +754,14 @@ double BounceSolution::CalcTempAtFalseVacFraction(const double &false_vac_frac)
     this->SetStoredTemp(T_middle); // update temp storage for inner integral
     double IatT = prefac * Nintegrate_Outer(*this).result;
 
-    bool numerically_unstable = false;
-
-    while (std::abs(T_up / T_down - 1) > 1e-8)
+    while ((std::abs(T_up / T_down - 1) >
+            RelativeTemperatureInCalcTempAtFalseVacFraction *
+                MarginOfCalcTempAtFalseVacFractionBeforeFailure) or
+           (not almost_the_same(
+               int_at_false_vac_frac,
+               IatT,
+               RelativeErrorInCalcTempAtFalseVacFraction *
+                   MarginOfCalcTempAtFalseVacFractionBeforeFailure)))
     {
       T_middle = (T_up + T_down) / 2.;
       this->SetStoredTemp(T_middle); // update temp storage for inner integral
@@ -775,18 +780,19 @@ double BounceSolution::CalcTempAtFalseVacFraction(const double &false_vac_frac)
         T_down = T_middle;
       }
 
-      if (std::abs(T_up / T_down - 1) < 1e-8 and
-          not almost_the_same(int_at_false_vac_frac, IatT))
+      // Condition for success
+      if (std::abs(T_up / T_down - 1) <
+              RelativeTemperatureInCalcTempAtFalseVacFraction and
+          almost_the_same(int_at_false_vac_frac,
+                          IatT,
+                          RelativeErrorInCalcTempAtFalseVacFraction))
       {
-        numerically_unstable = true;
+        res_Temp = T_middle;
         break;
       }
     }
-    if (not numerically_unstable)
-    {
-      res_Temp = T_middle;
-    }
   }
+  // Not numerically stable
   return res_Temp;
 }
 
@@ -1092,4 +1098,4 @@ double BounceSolution::GetInvTimeScale()
   return this->betaH;
 }
 
-} // namespace BSMPT
+} // namespace BSMPT