1 files changed, 62 insertions, 97 deletions
diff --git a/log-fourier_integrator.cpp b/log-fourier_integrator.cpp
index b15808b..531565c 100644
--- a/log-fourier_integrator.cpp
+++ b/log-fourier_integrator.cpp
@@ -93,108 +93,73 @@ int main(int argc, char* argv[]) {
   Real fac = 1;
   Real β = 0;
   while (β < βₘₐₓ) {
-    Real μ₁ = 0;
-    Real μ₂ = 0;
-    while (true) {
-      Real ΔC = 100;
-      Real ΔC₀ = 100;
-      unsigned it = 0;
-      while (ΔC > ε) {
-        std::vector<Real> dfC(N);
-        std::vector<Real> RddfC(N);
-        for (unsigned n = 0; n < N; n++) {
-          RddfC[n] = Rₜ[n] * ddf(λ, p, s, Cₜ[n]);
-          dfC[n] = df(λ, p, s, Cₜ[n]);
-        }
-        std::vector<Complex> RddfCt = fft.fourier(RddfC, false);
-        std::vector<Complex> dfCt = fft.fourier(dfC, true);
-
-        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));
+    Real ΔC = 100;
+    Real ΔC₀ = 100;
+    unsigned it = 0;
+    while (ΔC > ε) {
+      std::vector<Real> dfC(N);
+      std::vector<Real> RddfC(N);
+      for (unsigned n = 0; n < N; n++) {
+        RddfC[n] = Rₜ[n] * ddf(λ, p, s, Cₜ[n]);
+        dfC[n] = df(λ, p, s, Cₜ[n]);
+      }
+      std::vector<Complex> RddfCt = fft.fourier(RddfC, false);
+      std::vector<Complex> dfCt = fft.fourier(dfC, true);
+
+      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);
-          Ĉₜ₊₁[n] = (2 * Γ₀ * std::conj(Ȓₜ₊₁[n]) / (1 + pow(τ₀ * fft.ν(n), 2)) + pow(β, 2) * (RddfCt[n] * Ĉₜ[n] + dfCt[n] * std::conj(Ȓₜ₊₁[n]))) / (μ + 1i * fft.ν(n));
-        }
-
-        std::vector<Real> Rₜ₊₁ = fft.inverse(Ȓₜ₊₁);
-
-        for (unsigned n = 0; n < N; n++) {
-          RddfC[n] = Rₜ₊₁[n] * ddf(λ, p, s, Cₜ[n]);
-        }
-        RddfCt = fft.fourier(RddfC, false);
-        for (unsigned n = 0; n < N; n++) {
-          Ĉₜ₊₁[n] = (2 * Γ₀ * std::conj(Ȓₜ₊₁[n]) / (1 + pow(τ₀ * fft.ν(n), 2)) + pow(β, 2) * (RddfCt[n] * Ĉₜ[n] + dfCt[n] * std::conj(Ȓₜ₊₁[n]))) / (μ + 1i * fft.ν(n));
-        }
-        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 << " " << γ << " " << Cₜ[0];
-        std::cerr << "\r";
+        Ĉₜ₊₁[n] = (2 * Γ₀ * std::conj(Ȓₜ₊₁[n]) / (1 + pow(τ₀ * fft.ν(n), 2)) + pow(β, 2) * (RddfCt[n] * Ĉₜ[n] + dfCt[n] * std::conj(Ȓₜ₊₁[n]))) / (μ + 1i * fft.ν(n));
       }
-      if (std::isnan(Cₜ[0])) {
-        Cₜ = Cₜ₋₁;
-        Rₜ = Rₜ₋₁;
-        Ĉₜ = Ĉₜ₋₁;
-        Ȓₜ = Ȓₜ₋₁;
-        μ *= 1.1;
-        fac /= 2;
-        μ₁ = 0;
-        μ₂ = 0;
-      } else {
-      if (pow(Cₜ[0] - 1, 2) < ε) {
-        break;
+
+      std::vector<Real> Rₜ₊₁ = fft.inverse(Ȓₜ₊₁);
+
+      for (unsigned n = 0; n < N; n++) {
+        RddfC[n] = Rₜ₊₁[n] * ddf(λ, p, s, Cₜ[n]);
       }
-      if (μ₁ == 0 || μ₂ == 0) {
-        if (Cₜ[0] > 1 && μ₁ == 0) {
-          /* We found a lower bound */
-          μ₁ = μ;
-        }
-        if (Cₜ[0] < 1 && μ₂ == 0) {
-          /* We found an upper bound */
-          μ₂ = μ;
-        }
-        μ *= fac*std::tanh(Cₜ[0]-1)+1;
+      RddfCt = fft.fourier(RddfC, false);
+      for (unsigned n = 0; n < N; n++) {
+        Ĉₜ₊₁[n] = (2 * Γ₀ * std::conj(Ȓₜ₊₁[n]) / (1 + pow(τ₀ * fft.ν(n), 2)) + pow(β, 2) * (RddfCt[n] * Ĉₜ[n] + dfCt[n] * std::conj(Ȓₜ₊₁[n]))) / (μ + 1i * fft.ν(n));
+      }
+      std::vector<Real> Cₜ₊₁ = fft.inverse(Ĉₜ₊₁);
+
+      μ *= pow(tanh(Cₜ₊₁[0]-1)+1, 0.05);
+
+      Δ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 {
-        /* Once the bounds are set, we can use bisection */
-        if (Cₜ[0] > 1) {
-          μ₁ = μ;
-        } else {
-          μ₂ = μ;
-        }
-        μ = (μ₁ + μ₂) / 2;
+        it++;
       }
+
+      if (it > 100) {
+        γ = std::max(0.5 * γ, 1e-3);
+        it = 0;
+        ΔC₀ = ΔC;
       }
+      */
+
+      std::cerr << "\x1b[2K" << "\r";
+      std::cerr << β << " " << μ << " " << ΔC << " " << γ << " " << Cₜ[0];
     }
 
     /* Integrate the energy using Simpson's rule */
@@ -213,9 +178,9 @@ int main(int argc, char* argv[]) {
     }
     E *= β;
 
+      std::cerr << "\x1b[2K" << "\r";
     std::cerr << β << " " << μ << " " << Ĉₜ[0].real() << " " << E << " " << γ << std::endl;
     β += Δβ;
-    μ *= 2;
     Cₜ₋₁ = Cₜ;
     Rₜ₋₁ = Rₜ;
     Ĉₜ₋₁ = Ĉₜ;