Change convention to fixed Γ₀, and changed behavior on NaN

1 files changed, 9 insertions, 5 deletions

diff --git a/log-fourier_integrator.cpp b/log-fourier_integrator.cpp
index cf9819a..412b4d9 100644
--- a/log-fourier_integrator.cpp
+++ b/log-fourier_integrator.cpp
@@ -71,8 +71,11 @@ int main(int argc, char* argv[]) {
 
   LogarithmicFourierTransform fft(N, k, Δτ, 4);
 
-  Real Γ₀ = 1.0 + τ₀;
-  Real μₜ₋₁ = 1.0;
+  Real Γ₀ = 1.0;
+  Real μₜ₋₁ = Γ₀;
+  if (τ₀ > 0) {
+    μₜ₋₁ = (sqrt(1+4*Γ₀*τ₀)-1)/(2*τ₀);
+  }
 
   std::vector<Real> Cₜ₋₁(N);
   std::vector<Real> Rₜ₋₁(N);
@@ -82,10 +85,10 @@ int main(int argc, char* argv[]) {
   if (!loadData) {
     /* Start from the exact solution for β = 0 */
     for (unsigned n = 0; n < N; n++) {
-      if (τ₀ != 1) {
+      if (τ₀ > 0) {
         Cₜ₋₁[n] = Γ₀ * (exp(-μₜ₋₁ * fft.t(n)) - μₜ₋₁ * τ₀ * exp(-fft.t(n) / τ₀)) / (μₜ₋₁ - pow(μₜ₋₁, 3) * pow(τ₀, 2));
       } else {
-        Cₜ₋₁[n] = Γ₀ * exp(-fft.t(n)) * (1 + fft.t(n));
+        Cₜ₋₁[n] = Γ₀ * exp(-μₜ₋₁ * fft.t(n)) / μₜ₋₁;
       }
       Rₜ₋₁[n] = exp(-μₜ₋₁ * fft.t(n));
 
@@ -139,12 +142,13 @@ int main(int argc, char* argv[]) {
     }
 
     if (std::isnan(Cₜ[0])) {
+      Δβ /= 2;
+      β -= Δβ;
       Cₜ = Cₜ₋₁;
       Rₜ = Rₜ₋₁;
       Ĉₜ = Ĉₜ₋₁;
       Ȓₜ = Ȓₜ₋₁;
       μₜ = μₜ₋₁;
-      γ /= 2;
     } else {
       Real E = energy(fft, Cₜ, Rₜ, p, s, λ, β);