Better setting of μ bounds

1 files changed, 10 insertions, 9 deletions

diff --git a/log-fourier_integrator.cpp b/log-fourier_integrator.cpp
index 7736be8..110b9e4 100644
--- a/log-fourier_integrator.cpp
+++ b/log-fourier_integrator.cpp
@@ -88,9 +88,8 @@ int main(int argc, char* argv[]) {
   Real β = 0;
   while (β < βₘₐₓ) {
     β += Δβ;
-    Real μ₁ = μ * 0.99;
+    Real μ₁ = 0;
     Real μ₂ = 0;
-    μ *= 1.01;
     while (true) {
       Real ΔC = 100;
       Real ΔC₀ = 100;
@@ -147,16 +146,18 @@ int main(int argc, char* argv[]) {
       if (std::abs(Cₜ[0] - 1) < ε) {
         break;
       }
-      if (μ₂ < μ₁) {
-        /* Cₜ[0] will generally start too big. Until we find a value too small, we
-         * scale down proportional too Cₜ[0] */
-        if (Cₜ[0] > 1) {
-          μ *= sqrt(sqrt(std::abs(Cₜ[0])));
-        } else { /* Otherwise, we've found an upper bound for μ */
+      if (μ₁ == 0 || μ₂ == 0) {
+        if (μ₁ == 0 && Cₜ[0] > 1) {
+          /* We found a lower bound */
+          μ₁ = μ;
+        }
+        if (μ₂ == 0 && Cₜ[0] < 1) {
+          /* We found an upper bound */
           μ₂ = μ;
         }
+        μ *= sqrt(sqrt(std::abs(Cₜ[0])));
       } else {
-        /* Once μ₂ is set, we can use bisection */
+        /* Once the bounds are set, we can use bisection */
         if (Cₜ[0] > 1) {
           μ₁ = μ;
           μ = μ₁ + (μ₂ - μ₁) / 2;