1 files changed, 55 insertions, 50 deletions
diff --git a/log-fourier_integrator.cpp b/log-fourier_integrator.cpp
index 7569245..94987f3 100644
--- a/log-fourier_integrator.cpp
+++ b/log-fourier_integrator.cpp
@@ -87,59 +87,64 @@ int main(int argc, char* argv[]) {
 
   Real β = 0;
   while (β < βₘₐₓ) {
-    Real ΔC = 100;
-    Real ΔC₀ = 100;
-    unsigned it = 0;
-    while (ΔC > ε) {
-      std::vector<Real> RddfC(N);
-      for (unsigned n = 0; n < N; n++) {
-        RddfC[n] = Rₜ[n] * ddf(λ, p, s, Cₜ[n]);
+    while (true) {
+      Real ΔC = 100;
+      Real ΔC₀ = 100;
+      unsigned it = 0;
+      while (ΔC > ε) {
+        std::vector<Real> RddfC(N);
+        for (unsigned n = 0; n < N; n++) {
+          RddfC[n] = Rₜ[n] * ddf(λ, p, s, Cₜ[n]);
+        }
+        std::vector<Complex> RddfCt = fft.fourier(RddfC, false);
+
+        std::vector<Complex> Ȓₜ₊₁(N);
+        std::vector<Complex> Ĉₜ₊₁(N);
+
+        for (unsigned n = 0; n < N; n++) {
+          Ȓₜ₊₁[n] = (1.0 + pow(β, 2) * RddfCt[n] * Ȓₜ[n]) / (μ + 1i * fft.ν(n));
+          Ĉₜ₊₁[n] = - 2 * Γ₀ * Ȓₜ₊₁[n].imag() / (1 + pow(τ₀ * fft.ν(n), 2)) / fft.ν(n);
+        }
+
+        std::vector<Real> Rₜ₊₁ = fft.inverse(Ȓₜ₊₁);
+        std::vector<Real> Cₜ₊₁ = fft.inverse(Ĉₜ₊₁);
+
+        ΔC = 0;
+        for (unsigned i = 0; i < N; i++) {
+          ΔC += std::norm(Ĉₜ[i] - Ĉₜ₊₁[i]);
+          ΔC += std::norm(Ȓₜ[i] - Ȓₜ₊₁[i]);
+        }
+        ΔC = sqrt(ΔC) / (2*N);
+
+        for (unsigned i = 0; i < N; i++) {
+          Cₜ[i] += γ * (Cₜ₊₁[i] - Cₜ[i]);
+          Rₜ[i] += γ * (Rₜ₊₁[i] - Rₜ[i]);
+          Ĉₜ[i] += γ * (Ĉₜ₊₁[i] - Ĉₜ[i]);
+          Ȓₜ[i] += γ * (Ȓₜ₊₁[i] - Ȓₜ[i]);
+        }
+
+        if (ΔC < ΔC₀) {
+          ΔC₀ = ΔC;
+          it = 0;
+          γ = std::min(1.001 * γ, 1.0);
+        } else {
+          it++;
+        }
+
+        if (it > 100) {
+          γ = std::max(0.5 * γ, 1e-3);
+          it = 0;
+          ΔC₀ = ΔC;
+        }
+
+        std::cerr << β << " " << μ << " " << ΔC << " " << γ;
+        std::cerr << "\r";
       }
-      std::vector<Complex> RddfCt = fft.fourier(RddfC, false);
-
-      std::vector<Complex> Ȓₜ₊₁(N);
-      std::vector<Complex> Ĉₜ₊₁(N);
-
-      for (unsigned n = 0; n < N; n++) {
-        Ȓₜ₊₁[n] = (1.0 + pow(β, 2) * RddfCt[n] * Ȓₜ[n]) / (μ + 1i * fft.ν(n));
-        Ĉₜ₊₁[n] = - 2 * Γ₀ * Ȓₜ₊₁[n].imag() / (1 + pow(τ₀ * fft.ν(n), 2)) / fft.ν(n);
-      }
-
-      std::vector<Real> Rₜ₊₁ = fft.inverse(Ȓₜ₊₁);
-      std::vector<Real> Cₜ₊₁ = fft.inverse(Ĉₜ₊₁);
-
-      ΔC = 0;
-      for (unsigned i = 0; i < N; i++) {
-        ΔC += std::norm(Ĉₜ[i] - Ĉₜ₊₁[i]);
-        ΔC += std::norm(Ȓₜ[i] - Ȓₜ₊₁[i]);
-      }
-      ΔC = sqrt(ΔC) / (2*N);
-
-      for (unsigned i = 0; i < N; i++) {
-        Cₜ[i] += γ * (Cₜ₊₁[i] - Cₜ[i]);
-        Rₜ[i] += γ * (Rₜ₊₁[i] - Rₜ[i]);
-        Ĉₜ[i] += γ * (Ĉₜ₊₁[i] - Ĉₜ[i]);
-        Ȓₜ[i] += γ * (Ȓₜ₊₁[i] - Ȓₜ[i]);
-      }
-
-      μ = γ * μ * Cₜ[0] + (1-γ) * μ;
-
-      if (ΔC < ΔC₀) {
-        ΔC₀ = ΔC;
-        it = 0;
-        γ = std::min(1.001 * γ, 1.0);
+      if (std::abs(Cₜ[0] - 1) < ε) {
+        break;
       } else {
-        it++;
+        μ *= Cₜ[0];
       }
-
-      if (it > 100) {
-        γ = std::max(0.5 * γ, 1e-3);
-        it = 0;
-        ΔC₀ = ΔC;
-      }
-
-      std::cerr << β << " " << μ << " " << ΔC << " " << γ;
-      std::cerr << "\r";
     }
 
     /* Integrate the energy using Simpson's rule */