1 files changed, 114 insertions, 0 deletions
diff --git a/log-fourier_integrator.cpp b/log-fourier_integrator.cpp
new file mode 100644
index 0000000..247270b
--- /dev/null
+++ b/log-fourier_integrator.cpp
@@ -0,0 +1,114 @@
+#include "fourier.hpp"
+#include <getopt.h>
+#include <iostream>
+#include <fstream>
+
+int main(int argc, char* argv[]) {
+  unsigned p = 2;
+  unsigned s = 2;
+  Real λ = 0.5;
+  Real τ₀ = 0;
+  Real β₀ = 0;
+  Real yₘₐₓ = 0.5;
+  Real Δy = 0.05;
+
+  unsigned log2n = 8;
+  Real τₘₐₓ = 20;
+
+  unsigned maxIterations = 1000;
+  Real ε = 1e-14;
+  Real γ = 1;
+
+  int opt;
+
+  while ((opt = getopt(argc, argv, "p:s:2:T:t:0:y:d:I:g:l")) != -1) {
+    switch (opt) {
+    case 'p':
+      p = atoi(optarg);
+      break;
+    case 's':
+      s = atoi(optarg);
+      break;
+    case '2':
+      log2n = atoi(optarg);
+      break;
+    case 'T':
+      τₘₐₓ = atof(optarg);
+      break;
+    case 't':
+      τ₀ = atof(optarg);
+      break;
+    case '0':
+      β₀ = atof(optarg);
+      break;
+    case 'y':
+      yₘₐₓ = atof(optarg);
+      break;
+    case 'd':
+      Δy = atof(optarg);
+      break;
+    case 'I':
+      maxIterations = (unsigned)atof(optarg);
+      break;
+    case 'g':
+      γ = atof(optarg);
+      break;
+    default:
+      exit(1);
+    }
+  }
+
+  unsigned N = pow(2, log2n);
+
+  Real Δτ = 1e-2;
+  Real Δω = 1e-2;
+  Real Δs = 1e-2;
+  Real k = 0.1;
+
+  Real Γ₀ = 1.0;
+  Real μ = Γ₀;
+  if (τ₀ > 0) {
+    μ = (sqrt(1+4*Γ₀*τ₀)-1)/(2*τ₀);
+  }
+
+  std::vector<Real> C(N);
+  std::vector<Real> R(N);
+
+  LogarithmicFourierTransform fft(N, k, Δω, Δτ, Δs);
+
+  for (unsigned n = 0; n < N; n++) {
+    if (τ₀ > 0) {
+      C[n] = Γ₀ * (exp(-μ * fft.t(n)) - μ * τ₀ * exp(-fft.t(n) / τ₀)) / (μ - pow(μ, 3) * pow(τ₀, 2));
+    } else {
+      C[n] = Γ₀ * exp(-μ * fft.t(n)) / μ;
+    }
+    R[n] = exp(-μ * fft.t(n));
+  }
+
+  std::vector<Complex> Ct = fft.fourier(C, true);
+  std::vector<Complex> Rt = fft.fourier(R, false);
+
+  /*
+  for (unsigned n = 0; n < N; n++) {
+    std::cout << fft.t(n) << " " << C[n] << std::endl;
+  }
+  */
+
+  for (unsigned n = 0; n < N; n++) {
+    std::cout << fft.ν(n) << " " << Ct[n].real() << std::endl;
+  }
+
+  /*
+  // start from the exact solution for τ₀ = 0
+  for (unsigned i = 0; i < N + 1; i++) {
+    Real ω = i * Δω;
+    Ct[i] = 2 * Γ₀ / (pow(μ, 2) + pow(ω, 2)) / (1 + pow(τ₀ * ω, 2));
+    Rt[i] = 1.0 / (μ + 1i * ω);
+    Γt[i] = Γ₀ / (1 + pow(τ₀ * ω, 2));
+  }
+  C = fft.inverse(Ct);
+  R = fft.inverse(Rt);
+  */
+
+  return 0;
+}