Started on new iteration scheme

2 files changed, 115 insertions, 2 deletions

diff --git a/Makefile b/Makefile
index 3519e50..f245260 100644
--- a/Makefile
+++ b/Makefile
@@ -1,6 +1,6 @@
-all: walk correlations integrator
+all: walk correlations integrator fourier_integrator
 
-CC := clang++ -std=c++17 -Wno-mathematical-notation-identifier-extension -O3 -march=native -mtune=native -fopenmp -g
+CC := clang++ -std=c++17 -Wno-mathematical-notation-identifier-extension -O3 -march=native -mtune=native -fopenmp -lfftw3
 
 walk: walk.cpp
 	$(CC) walk.cpp -o walk
@@ -10,3 +10,6 @@ correlations: correlations.cpp
 
 integrator: integrator.cpp
 	$(CC) integrator.cpp -o integrator
+
+fourier_integrator: fourier_integrator.cpp
+	$(CC) fourier_integrator.cpp -o fourier_integrator
diff --git a/fourier_integrator.cpp b/fourier_integrator.cpp
new file mode 100644
index 0000000..015a4b4
--- /dev/null
+++ b/fourier_integrator.cpp
@@ -0,0 +1,110 @@
+#include <getopt.h>
+#include <vector>
+#include <cmath>
+#include <iostream>
+#include <fftw3.h>
+#include <complex>
+
+using Real = double;
+using Complex = std::complex<Real>;
+using namespace std::complex_literals;
+
+unsigned p = 2;
+
+Real f(Real q) {
+  return 0.5 * pow(q, p);
+}
+
+Real df(Real q) {
+  return 0.5 * p * pow(q, p - 1);
+}
+
+Real ddf(Real q) {
+  return 0.5 * p * (p - 1) * pow(q, p - 2);
+}
+
+
+
+int main(int argc, char* argv[]) {
+  Real Δω = 1e-3;
+  Real Δτ = 1e-3;
+  Real τ₀ = 0;
+  Real y = 0.5;
+  unsigned iterations = 10;
+
+  int opt;
+
+  while ((opt = getopt(argc, argv, "d:T:t:y:I:")) != -1) {
+    switch (opt) {
+    case 'd':
+      Δω = atof(optarg);
+      break;
+    case 'T':
+      Δτ = atof(optarg);
+      break;
+    case 't':
+      τ₀ = atof(optarg);
+      break;
+    case 'y':
+      y = atof(optarg);
+      break;
+    case 'I':
+      iterations = (unsigned)atof(optarg);
+      break;
+    default:
+      exit(1);
+    }
+  }
+
+  Real z = (-1+sqrt(1+2*τ₀)) / (2 * τ₀);
+  Real Γ₀ = 1;
+
+  Real ωₘₐₓ = 1 / Δτ;
+  unsigned N = 2 * ωₘₐₓ / Δω + 1;
+  unsigned n = ωₘₐₓ / Δω + 1;
+
+  std::vector<Complex> Ĉ(N / 2 + 1);
+  std::vector<Real> C(N);
+
+  std::vector<Complex> Ř(N / 2 + 1);
+  std::vector<Real> R(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));
+  }
+
+  /*
+  for (unsigned it = 0; it < iterations; it++) {
+    for (unsigned i = 0; i < n; i++) {
+      Real ω = i * Δω;
+      ĉ[i] = Γ₀ / ((pow(z, 2) + pow(ω, 2)) * (1 + pow(τ₀ * ω ,2)));
+    }
+  }
+
+    ř[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);
+
+  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));
+  }
+
+  fftw_execute(test);
+
+  for (unsigned i = 0; i < Ĉ.size(); i++) {
+    std::cout << i * Δτ * M_PI << " " << Ĉ[i].real() << std::endl;
+  }
+
+  /*
+  for (unsigned i = 0; i < Ĉ.size(); i++) {
+    std::cout << i * Δω << " " << Ĉ[i].real() * Δω / (2 * M_PI) << std::endl;
+  }
+
+  */
+  return 0;
+}