Shrink Δβ if there is a NaN

1 files changed, 7 insertions, 2 deletions

diff --git a/log-fourier_integrator.cpp b/log-fourier_integrator.cpp
index 2f80f4c..8a2fe5c 100644
--- a/log-fourier_integrator.cpp
+++ b/log-fourier_integrator.cpp
@@ -94,6 +94,7 @@ int main(int argc, char* argv[]) {
   std::vector<Real> Rₜ₋₁(N);
   std::vector<Complex> Ĉₜ₋₁(N);
   std::vector<Complex> Ȓₜ₋₁(N);
+  std::vector<Real> Γ(N);
 
   if (!loadData) {
     /* Start from the exact solution for β = 0 */
@@ -111,6 +112,7 @@ int main(int argc, char* argv[]) {
 
       Ĉₜ₋₁[n] = 2 * Γ₀ / (pow(μ₀, 2) + pow(fft.ν(n), 2)) / (1 + pow(τ₀ * fft.ν(n), 2));
       Ȓₜ₋₁[n] = (Real)1.0 / (μ₀ + II * fft.ν(n));
+      Γ[n] = Γ₀ / (1 + std::pow(τ₀ * fft.ν(n), 2));
     }
   } else {
     logFourierLoad(Cₜ₋₁, Rₜ₋₁, Ĉₜ₋₁, Ȓₜ₋₁, p, s, λ, τ₀, β₀, log2n, Δτ, logShift);
@@ -138,7 +140,7 @@ int main(int argc, char* argv[]) {
 
       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();
+          Ĉₜ₊₁[n] = ((2 * Γ[n] * std::conj(Ȓₜ[n]) + std::pow(β, 2) * (RddfCt[n] * Ĉₜ[n] + dfCt[n] * std::conj(Ȓₜ[n]))) / (μₜ + II * fft.ν(n))).real();
         }
         C₀ = C0(fft, Ĉₜ₊₁);
         if (C₀ > 1) {
@@ -155,7 +157,7 @@ int main(int argc, char* argv[]) {
 
       ΔCₜ = 0;
       for (unsigned n = 0; n < N; n++) {
-        ΔCₜ += std::norm((2 * Γ₀ * std::conj(Ȓₜ[n]) / (1 + std::pow(τ₀ * fft.ν(n), 2)) + std::pow(β, 2) * (RddfCt[n] * Ĉₜ[n] + dfCt[n] * std::conj(Ȓₜ[n]))) - Ĉₜ[n] * (μₜ + II * fft.ν(n)));
+        ΔCₜ += std::norm((2 * Γ[n] * std::conj(Ȓₜ[n]) + std::pow(β, 2) * (RddfCt[n] * Ĉₜ[n] + dfCt[n] * std::conj(Ȓₜ[n]))) - Ĉₜ[n] * (μₜ + II * fft.ν(n)));
         ΔCₜ += std::norm(((Real)1.0 + std::pow(β, 2) * RddfCt[n] * Ȓₜ[n]) - Ȓₜ[n] * (μₜ + II * fft.ν(n)));
       }
       ΔCₜ = sqrt(ΔCₜ) / (2*N);
@@ -179,6 +181,9 @@ int main(int argc, char* argv[]) {
     }
 
     if (std::isnan(Cₜ[0])) {
+      β -= Δβ;
+      Δβ *= 0.1;
+      β += Δβ;
       Cₜ = Cₜ₋₁;
       Rₜ = Rₜ₋₁;
       Ĉₜ = Ĉₜ₋₁;