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

2 files changed, 10 insertions, 6 deletions

diff --git a/log-fourier.cpp b/log-fourier.cpp
index f7d6f4b..47d5c5c 100644
--- a/log-fourier.cpp
+++ b/log-fourier.cpp
@@ -173,7 +173,7 @@ std::tuple<std::vector<Complex>, std::vector<Complex>> RddfCtdfCt(LogarithmicFou
 
 Real estimateZ(LogarithmicFourierTransform& fft, const std::vector<Real>& C, const std::vector<Complex>& Ct, const std::vector<Real>& R, const std::vector<Complex>& Rt, unsigned p, unsigned s, Real λ, Real τ₀, Real β) {
   auto [RddfCt, dfCt] = RddfCtdfCt(fft, C, R, p, s, λ);
-  Real Γ₀ = 1 + τ₀;
+  Real Γ₀ = 1.0;
 
   return ((2 * Γ₀ * std::conj(Rt[0]) + pow(β, 2) * (RddfCt[0] * Ct[0] + dfCt[0] * std::conj(Rt[0]))) / Ct[0]).real();
 }
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, λ, β);