Got fourier integrator working

2 files changed, 92 insertions, 24 deletions

diff --git a/Makefile b/Makefile
index f245260..21f3def 100644
--- a/Makefile
+++ b/Makefile
@@ -1,6 +1,6 @@
 all: walk correlations integrator fourier_integrator
 
-CC := clang++ -std=c++17 -Wno-mathematical-notation-identifier-extension -O3 -march=native -mtune=native -fopenmp -lfftw3
+CC := clang++ -std=c++17 -Wno-mathematical-notation-identifier-extension -O3 -march=native -mtune=native -fopenmp -lfftw3 -g
 
 walk: walk.cpp
 	$(CC) walk.cpp -o walk
diff --git a/fourier_integrator.cpp b/fourier_integrator.cpp
index 015a4b4..9c97f1b 100644
--- a/fourier_integrator.cpp
+++ b/fourier_integrator.cpp
@@ -23,7 +23,46 @@ Real ddf(Real q) {
   return 0.5 * p * (p - 1) * pow(q, p - 2);
 }
 
+class FourierTransform {
+private:
+  std::vector<Real> a;
+  std::vector<Complex> â;
+  fftw_plan plan_r2c;
+  fftw_plan plan_c2r;
+  Real Δω;
+  Real Δτ;
+public:
+  FourierTransform(unsigned n, Real Δω, Real Δτ) : a(2 * n), â(n + 1), Δω(Δω), Δτ(Δτ) {
+    plan_r2c = fftw_plan_dft_r2c_1d(2 * n, a.data(), reinterpret_cast<fftw_complex*>(â.data()), 0);
+    plan_c2r = fftw_plan_dft_c2r_1d(2 * n, reinterpret_cast<fftw_complex*>(â.data()), a.data(), 0);
+  }
+
+  ~FourierTransform() {
+    fftw_destroy_plan(plan_r2c);
+    fftw_destroy_plan(plan_c2r);
+    fftw_cleanup();
+  }
+
+  std::vector<Complex> fourier(const std::vector<Real>& c) {
+    a = c;
+    fftw_execute(plan_r2c);
+    std::vector<Complex> ĉ(â.size());
+    for (unsigned i = 0; i < â.size(); i++) {
+      ĉ[i] = â[i] * (Δτ * M_PI);
+    }
+    return ĉ;
+  }
 
+  std::vector<Real> inverse(const std::vector<Complex>& ĉ) {
+    â = ĉ;
+    fftw_execute(plan_c2r);
+    std::vector<Real> c(a.size());
+    for (unsigned i = 0; i < a.size(); i++) {
+      c[i] = a[i] * (Δω / (2 * M_PI));
+    }
+    return c;
+  }
+};
 
 int main(int argc, char* argv[]) {
   Real Δω = 1e-3;
@@ -60,44 +99,73 @@ int main(int argc, char* argv[]) {
   Real Γ₀ = 1;
 
   Real ωₘₐₓ = 1 / Δτ;
-  unsigned N = 2 * ωₘₐₓ / Δω + 1;
-  unsigned n = ωₘₐₓ / Δω + 1;
+  unsigned n = ωₘₐₓ / Δω;
 
-  std::vector<Complex> Ĉ(N / 2 + 1);
-  std::vector<Real> C(N);
+  std::vector<Real> C(2 * n);
+  std::vector<Real> R(2 * n);
 
-  std::vector<Complex> Ř(N / 2 + 1);
-  std::vector<Real> R(N);
+  FourierTransform fft(n, Δω, Δτ);
 
   for (unsigned i = 0; i < n; i++) {
     Real τ = i * Δτ * M_PI;
     C[i] = Γ₀ / 2 * (exp(-z * τ) - z * τ₀ * exp(-τ / τ₀)) / (z - pow(z, 3) * pow(τ₀, 2));
-    C[N - 1 - i] = Γ₀ / 2 * (exp(-z * τ) - z * τ₀ * exp(-τ / τ₀)) / (z - pow(z, 3) * pow(τ₀, 2));
+    if (i > 0) {
+      C[2 * n - i] = C[i];
+    }
+    R[i] = exp(-z * τ);
   }
 
-  /*
+  std::vector<Complex> Ct = fft.fourier(C);
+  std::vector<Complex> Rt = fft.fourier(R);
+
   for (unsigned it = 0; it < iterations; it++) {
-    for (unsigned i = 0; i < n; i++) {
+    std::vector<Real> RddfC(C.size());
+    for (unsigned i = 0; i < C.size(); i++) {
+      RddfC[i] = R[i] * ddf(C[i]);
+    }
+    std::vector<Complex> RddfCt = fft.fourier(RddfC);
+
+    std::vector<Real> dfC(C.size());
+    for (unsigned i = 0; i < C.size(); i++) {
+      dfC[i] = df(C[i]);
+    }
+    std::vector<Complex> dfCt = fft.fourier(dfC);
+
+    std::vector<Complex> Ctnew(Ct.size());
+    std::vector<Complex> Rtnew(Rt.size());
+
+    for (unsigned i = 0; i < Rt.size(); i++) {
       Real ω = i * Δω;
-      ĉ[i] = Γ₀ / ((pow(z, 2) + pow(ω, 2)) * (1 + pow(τ₀ * ω ,2)));
+      Rtnew[i] = (1.0 + pow(y, 2) * RddfCt[i] * Rt[i]) / (z + 1i * ω);
     }
-  }
 
-    ř[i] = (z - 1i * ω) / (pow(z, 2) + pow(ω, 2));
-  std::vector<Real> c(n);
-  */
-  fftw_plan test = fftw_plan_dft_r2c_1d(N, C.data(), reinterpret_cast<fftw_complex*>(Ĉ.data()), 0);
+    Rt = Rtnew;
 
-  for (unsigned i = 0; i < n; i++) {
-    Real τ = i * Δτ * M_PI;
-    C[i] = Γ₀ / 2 * (exp(-z * τ) - z * τ₀ * exp(-τ / τ₀)) / (z - pow(z, 3) * pow(τ₀, 2));
-    C[N - 1 - i] = Γ₀ / 2 * (exp(-z * τ) - z * τ₀ * exp(-τ / τ₀)) / (z - pow(z, 3) * pow(τ₀, 2));
-  }
+    for (unsigned i = 0; i < Rt.size(); i++) {
+      Real ω = i * Δω;
+      Ctnew[i] = (Γ₀ * std::conj(Rt[i]) / (1 + pow(τ₀ * ω, 2)) + pow(y, 2) * (RddfCt[i] * Ct[i] + dfCt[i] * std::conj(Rt[i]))) / (z + 1i * ω);
+    }
 
-  fftw_execute(test);
+    Ct = Ctnew;
 
-  for (unsigned i = 0; i < Ĉ.size(); i++) {
-    std::cout << i * Δτ * M_PI << " " << Ĉ[i].real() << std::endl;
+    std::vector<Real> Cnew = fft.inverse(Ctnew);
+    std::vector<Real> Rnew = fft.inverse(Rtnew);
+
+    Real ΔC = 0;
+    for (unsigned i = 0; i < Cnew.size(); i++) {
+      ΔC += pow(Cnew[i] - C[i], 2);
+    }
+
+    z *= Cnew[0];
+
+    std::cerr << "Iteration " << it << ": " << sqrt(ΔC / C.size()) << " " << z << std::endl;
+
+    C = Cnew;
+    R = Rnew;
+  }
+
+  for (unsigned i = 0; i < C.size(); i++) {
+    std::cout << i * Δτ * M_PI << " " << C[i] << std::endl;
   }
 
   /*