1 files changed, 17 insertions, 37 deletions

diff --git a/log_integrator.cpp b/log_integrator.cpp
index 426af40..a8d9778 100644
--- a/log_integrator.cpp
+++ b/log_integrator.cpp
@@ -92,7 +92,6 @@ 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);
 
@@ -109,7 +108,6 @@ 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));
@@ -144,47 +142,23 @@ int main(int argc, char* argv[]) {
 
       for (unsigned i = 0; i < N; i++) {
         dC[i] += -μₜ * Cₜ[i];
-        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));
+        Real ΓR;
+        for (unsigned j = 0; j < N; j++) {
+          
         }
-        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);
 
-      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<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));
       }
+      std::vector<Real> Rₜ₊₁ = fft.inverse(Ȓₜ₊₁);
+      std::vector<Real> Cₜ₊₁ = fft.inverse(Ĉₜ₊₁);
 
-      μₜ /= 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];
-      }
+      μₜ *= pow(tanh(Cₜ₊₁[0]-1)+1, x);
 
       ΔCₜ = 0;
       for (unsigned i = 0; i < N; i++) {
@@ -209,6 +183,8 @@ 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";
@@ -219,6 +195,8 @@ int main(int argc, char* argv[]) {
       γ₀ /= 2;
       Cₜ = Cₜ₋₁;
       Rₜ = Rₜ₋₁;
+      Ĉₜ = Ĉₜ₋₁;
+      Ȓₜ = Ȓₜ₋₁;
       μₜ = μₜ₋₁;
     } else {
       Real E = energy(fft, Cₜ, Rₜ, p, s, λ, β);
@@ -231,6 +209,8 @@ int main(int argc, char* argv[]) {
       β += Δβ;
       Cₜ₋₁ = Cₜ;
       Rₜ₋₁ = Rₜ;
+      Ĉₜ₋₁ = Ĉₜ;
+      Ȓₜ₋₁ = Ȓₜ;
       μₜ₋₁ = μₜ;
     }
   }