From 7ea7b097f558475fa8f12cef2decd9edc5eb33be Mon Sep 17 00:00:00 2001
From: Jaron Kent-Dobias <jaron@kent-dobias.com>
Date: Tue, 13 May 2025 17:21:51 -0300
Subject: Adjust μ continuously
MIME-Version: 1.0
Content-Type: text/plain; charset=UTF-8
Content-Transfer-Encoding: 8bit

---
 log-fourier_integrator.cpp | 27 +++++++++++++++++++++++----
 1 file changed, 23 insertions(+), 4 deletions(-)

(limited to 'log-fourier_integrator.cpp')

diff --git a/log-fourier_integrator.cpp b/log-fourier_integrator.cpp
index 0e05366..5ff27a4 100644
--- a/log-fourier_integrator.cpp
+++ b/log-fourier_integrator.cpp
@@ -138,14 +138,33 @@ int main(int argc, char* argv[]) {
 
       std::vector<Complex> Ĉₜ₊₁(N);
       std::vector<Complex> Ȓₜ₊₁(N);
+
+      Real C₀ = 0;
+      Real μ₊ = 0;
+      Real μ₋ = 0;
+
+      while (std::abs(C₀ - 1) > ε) {
+        for (unsigned n = 0; n < N; n++) {
+          Ĉₜ₊₁[n] = ((2 * Γ₀ * std::conj(Ȓₜ[n]) / (1 + std::pow(τ₀ * fft.ν(n), 2)) + std::pow(β, 2) * (RddfCt[n] * Ĉₜ[n] + dfCt[n] * std::conj(Ȓₜ[n]))) / (μₜ + II * fft.ν(n))).real();
+        }
+        C₀ = C0(fft, Ĉₜ₊₁);
+        if (C₀ > 1) {
+          μ₋ = μₜ;
+        } else {
+          μ₊ = μₜ;
+        }
+        if (μ₋ > 0 && μ₊ > 0) {
+          μₜ = (μ₊ + μ₋) / 2;
+        } else {
+          μₜ *= C₀;
+        }
+      }
+
       for (unsigned n = 0; n < N; n++) {
         Ȓₜ₊₁[n] = ((Real)1.0 + std::pow(β, 2) * RddfCt[n] * Ȓₜ[n]) / (μₜ + II * fft.ν(n));
-        Ĉₜ₊₁[n] = ((2 * Γ₀ * std::conj(Ȓₜ[n]) / (1 + std::pow(τ₀ * fft.ν(n), 2)) + std::pow(β, 2) * (RddfCt[n] * Ĉₜ[n] + dfCt[n] * std::conj(Ȓₜ[n]))) / (μₜ + II * fft.ν(n))).real();
       }
-      std::vector<Real> Rₜ₊₁ = fft.inverse(Ȓₜ₊₁);
       std::vector<Real> Cₜ₊₁ = fft.inverse(Ĉₜ₊₁);
-
-      Real C₀ = Cₜ₊₁[0];
+      std::vector<Real> Rₜ₊₁ = fft.inverse(Ȓₜ₊₁);
 
       if (!std::isnan(Cₜ₊₁[0])) {
 
-- 
cgit v1.2.3-70-g09d2