1 files changed, 37 insertions, 17 deletions

diff --git a/log_integrator.cpp b/log_integrator.cpp
index a8d9778..426af40 100644
--- a/log_integrator.cpp
+++ b/log_integrator.cpp
@@ -92,6 +92,7 @@ int main(int argc, char* argv[]) {
 
   std::vector<Real> Cₜ₋₁(N);
   std::vector<Real> Rₜ₋₁(N);
+  std::vector<Real> Γ(N);
   std::vector<Complex> Ĉₜ₋₁(N);
   std::vector<Complex> Ȓₜ₋₁(N);
 
@@ -108,6 +109,7 @@ int main(int argc, char* argv[]) {
         Cₜ₋₁[n] = Γ₀ * exp(-μₜ₋₁ * fft.t(n)) / μₜ₋₁;
       }
       Rₜ₋₁[n] = exp(-μₜ₋₁ * fft.t(n));
+      Γ[n] = (Γ₀ / τ₀) * exp(-fft.t(n) / τ₀);
 
       Ĉₜ₋₁[n] = 2 * Γ₀ / (pow(μₜ₋₁, 2) + pow(fft.ν(n), 2)) / (1 + pow(τ₀ * fft.ν(n), 2));
       Ȓₜ₋₁[n] = (Real)1.0 / (μₜ₋₁ + II * fft.ν(n));
@@ -142,23 +144,47 @@ int main(int argc, char* argv[]) {
 
       for (unsigned i = 0; i < N; i++) {
         dC[i] += -μₜ * Cₜ[i];
-        Real ΓR;
-        for (unsigned j = 0; j < N; j++) {
-          
+        dR[i] += -μₜ * Rₜ[i];
+        for (unsigned j = 1; j < N - i; j++) {
+          Real Γᵢ₊ⱼ = (Γ[i + j] + Γ[i + j - 1]) / 2;
+          signed jR = round(log(fft.t(i+j) - fft.t(i)) / Δτ) + N / 2;
+          Real Rⱼ = jR >= N ? 0 : jR < 0 ? 1 : Rₜ[jR];
+          Real Cⱼ = jR >= N ? 0 : jR < 0 ? 1 : Cₜ[jR];
+          Real ΓR2 = Γᵢ₊ⱼ * Rⱼ;
+          Real dfCR = (dfC[i + j] * Rⱼ + dfC[i + j - 1] * Rⱼ) / 2;
+          Real RddfCC = (RddfC[i + j] * Cⱼ + RddfC[i + j - 1] * Cⱼ) / 2;
+          dC[i] += (fft.t(i+j) - fft.t(i+j - 1)) * (ΓR2 + pow(β, 2) * (dfCR + RddfCC));
         }
+        for (unsigned j = 1; j < i; j++) {
+          signed jR = round(log(fft.t(i) - fft.t(i-j)) / Δτ) + N / 2;
+          Real Rⱼ = jR >= N ? 0 : jR < 0 ? 1 : Rₜ[jR];
+          Real Cⱼ = jR >= N ? 0 : jR < 0 ? 1 : Cₜ[jR];
+          dC[i] += (fft.t(i-j) - fft.t(i-j - 1)) * pow(β, 2) * (RddfC[i - 1 - j] * Cⱼ + RddfC[i - 1 - j+1] * Cⱼ) / 2;
+          dR[i] += (fft.t(i-j) - fft.t(i-j - 1)) * pow(β, 2) * (RddfC[i - 1 - j] * Rⱼ + RddfC[i - 1 - j+1] * Rⱼ) / 2;
+        }
+
+        if (dC[i] > 0) dC[i] = 0;
+        if (dR[i] > 0) dR[i] = 0;
       }
 
+      std::vector<Real> Rₜ₊₁(N);
+      std::vector<Real> Cₜ₊₁(N);
 
-      std::vector<Complex> Ĉₜ₊₁(N);
-      std::vector<Complex> Ȓₜ₊₁(N);
-      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));
+      Rₜ₊₁[0] = 1;
+      Cₜ₊₁[N - 1] = 0;
+
+      for (unsigned i = 1; i < N; i++) {
+        Cₜ₊₁[N - 1 - i] = Cₜ₊₁[N - i] - (fft.t(N - i) - fft.t(N - i - 1)) * (dC[N - i] + dC[N - 1 - i]) / 2;
+        Rₜ₊₁[i] = Rₜ₊₁[i - 1] + (fft.t(i) - fft.t(i - 1)) * (dR[i - 1] + dR[i]) / 2;
       }
-      std::vector<Real> Rₜ₊₁ = fft.inverse(Ȓₜ₊₁);
-      std::vector<Real> Cₜ₊₁ = fft.inverse(Ĉₜ₊₁);
 
-      μₜ *= pow(tanh(Cₜ₊₁[0]-1)+1, x);
+      μₜ /= pow(tanh(Cₜ₊₁[0]-1)+1, x);
+
+      Real C₀ = Cₜ₊₁[0];
+
+      for (unsigned i = 0; i < N; i++) {
+        Rₜ₊₁[i] = Rₜ₊₁[i] < 0 ? 0 : Rₜ₊₁[i];
+      }
 
       ΔCₜ = 0;
       for (unsigned i = 0; i < N; i++) {
@@ -183,8 +209,6 @@ int main(int argc, char* argv[]) {
       for (unsigned i = 0; i < N; i++) {
         Cₜ[i] += γ * (Cₜ₊₁[i] - Cₜ[i]);
         Rₜ[i] += γ * (Rₜ₊₁[i] - Rₜ[i]);
-        Ĉₜ[i] += γ * (Ĉₜ₊₁[i] - Ĉₜ[i]);
-        Ȓₜ[i] += γ * (Ȓₜ₊₁[i] - Ȓₜ[i]);
       }
 
       std::cerr << "\x1b[2K" << "\r";
@@ -195,8 +219,6 @@ int main(int argc, char* argv[]) {
       γ₀ /= 2;
       Cₜ = Cₜ₋₁;
       Rₜ = Rₜ₋₁;
-      Ĉₜ = Ĉₜ₋₁;
-      Ȓₜ = Ȓₜ₋₁;
       μₜ = μₜ₋₁;
     } else {
       Real E = energy(fft, Cₜ, Rₜ, p, s, λ, β);
@@ -209,8 +231,6 @@ int main(int argc, char* argv[]) {
       β += Δβ;
       Cₜ₋₁ = Cₜ;
       Rₜ₋₁ = Rₜ;
-      Ĉₜ₋₁ = Ĉₜ;
-      Ȓₜ₋₁ = Ȓₜ;
       μₜ₋₁ = μₜ;
     }
   }