From cf90ba82e35c4d8578825c72876ee739b2e8d95d Mon Sep 17 00:00:00 2001
From: Jaron Kent-Dobias <jaron@kent-dobias.com>
Date: Sat, 19 Apr 2025 17:27:57 -0300
Subject: Many tweaks

---
 log-fourier_integrator.cpp | 50 +++++++++++++++++++++++++++++++++++-----------
 1 file changed, 38 insertions(+), 12 deletions(-)

diff --git a/log-fourier_integrator.cpp b/log-fourier_integrator.cpp
index 7b689f1..ca378bb 100644
--- a/log-fourier_integrator.cpp
+++ b/log-fourier_integrator.cpp
@@ -67,24 +67,29 @@ int main(int argc, char* argv[]) {
     μ = (sqrt(1+4*Γ₀*τ₀)-1)/(2*τ₀);
   }
 
-  std::vector<Real> Cₜ(N);
-  std::vector<Real> Rₜ(N);
-  std::vector<Complex> Ĉₜ(N);
-  std::vector<Complex> Ȓₜ(N);
+  std::vector<Real> Cₜ₋₁(N);
+  std::vector<Real> Rₜ₋₁(N);
+  std::vector<Complex> Ĉₜ₋₁(N);
+  std::vector<Complex> Ȓₜ₋₁(N);
 
   /* Start from the exact solution for β = 0 */
   for (unsigned n = 0; n < N; n++) {
     if (τ₀ > 0) {
-      Cₜ[n] = Γ₀ * (exp(-μ * fft.t(n)) - μ * τ₀ * exp(-fft.t(n) / τ₀)) / (μ - pow(μ, 3) * pow(τ₀, 2));
+      Cₜ₋₁[n] = Γ₀ * (exp(-μ * fft.t(n)) - μ * τ₀ * exp(-fft.t(n) / τ₀)) / (μ - pow(μ, 3) * pow(τ₀, 2));
     } else {
-      Cₜ[n] = Γ₀ * exp(-μ * fft.t(n)) / μ;
+      Cₜ₋₁[n] = Γ₀ * exp(-μ * fft.t(n)) / μ;
     }
-    Rₜ[n] = exp(-μ * fft.t(n));
+    Rₜ₋₁[n] = exp(-μ * fft.t(n));
 
-    Ĉₜ[n] = 2 * Γ₀ / (pow(μ, 2) + pow(fft.ν(n), 2)) / (1 + pow(τ₀ * fft.ν(n), 2));
-    Ȓₜ[n] = 1.0 / (μ + 1i * fft.ν(n));
+    Ĉₜ₋₁[n] = 2 * Γ₀ / (pow(μ, 2) + pow(fft.ν(n), 2)) / (1 + pow(τ₀ * fft.ν(n), 2));
+    Ȓₜ₋₁[n] = 1.0 / (μ + 1i * fft.ν(n));
   }
 
+  std::vector<Real> Cₜ = Cₜ₋₁;
+  std::vector<Real> Rₜ = Rₜ₋₁;
+  std::vector<Complex> Ĉₜ = Ĉₜ₋₁;
+  std::vector<Complex> Ȓₜ = Ȓₜ₋₁;
+
   Real β = 0;
   while (β < βₘₐₓ) {
     Real μ₁ = 0;
@@ -92,6 +97,7 @@ int main(int argc, char* argv[]) {
     while (true) {
       Real ΔC = 100;
       Real ΔC₀ = 100;
+      unsigned it = 0;
       while (ΔC > ε) {
         std::vector<Real> dfC(N);
         std::vector<Real> RddfC(N);
@@ -107,8 +113,8 @@ int main(int argc, char* argv[]) {
 
         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));
+          Ĉₜ₊₁[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(Ȓₜ₊₁);
@@ -130,14 +136,29 @@ int main(int argc, char* argv[]) {
 
         if (ΔC < ΔC₀) {
           ΔC₀ = ΔC;
-          γ = std::min(1.01 * γ, 1.0);
+          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";
       }
+      if (std::isnan(Cₜ[0])) {
+        Cₜ = Cₜ₋₁;
+        Rₜ = Rₜ₋₁;
+        Ĉₜ = Ĉₜ₋₁;
+        Ȓₜ = Ȓₜ₋₁;
+        Δβ /= 2;
+        β -= Δβ;
+      } else {
       if (pow(Cₜ[0] - 1, 2) < ε) {
         break;
       }
@@ -160,6 +181,7 @@ int main(int argc, char* argv[]) {
         }
         μ = (μ₁ + μ₂) / 2;
       }
+      }
     }
 
     /* Integrate the energy using Simpson's rule */
@@ -180,6 +202,10 @@ int main(int argc, char* argv[]) {
 
     std::cerr << β << " " << μ << " " << Ĉₜ[0].real() << " " << E << " " << γ << std::endl;
     β += Δβ;
+    Cₜ₋₁ = Cₜ;
+    Rₜ₋₁ = Rₜ;
+    Ĉₜ₋₁ = Ĉₜ;
+    Ȓₜ₋₁ = Ȓₜ;
   }
 
   for (unsigned i = 0; i < N; i++) {
-- 
cgit v1.2.3-70-g09d2